FESOM2-JAX v1.0: a differentiable shadow of the ocean-sea-ice model FESOM2, cast onto GPUs

FESOM2-JAX v1.0: a differentiable shadow of the ocean–sea-ice model FESOM2, cast onto GPUs Nikolay V. Koldunov1 , Sergey Danilov1 , Suvarchal Cheedela1 , Dmitry Sidorenko1 , Sebastian Beyer1 , Patrick Scholz1 , Ivan Kuznetsov1 , Jan Streffing1 , Aleksei Koldunov1 , Dmitrii Pantiukhin1 , Svetlana N. Loza1 , and Thomas Jung1,2 1 2

Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Bremerhaven, Germany Department of Physics and Electrical Engineering, University of Bremen, Bremen, Germany

arXiv:2608.01546v1 [physics.ao-ph] 2 Aug 2026

Correspondence: Nikolay V. Koldunov (nikolay.koldunov@awi.de) Abstract. We present FESOM2-JAX, a Python re-implementation of the Finite-volumE Sea ice–Ocean Model (FESOM2) in JAX. The model retains the unstructured-mesh, cell-vertex finite-volume formulation of the original, runs unchanged from a laptop CPU to 256 GPUs, and is end-to-end differentiable. FESOM2-JAX is a code shadow of the Fortran model: a projection onto the Python ecosystem, translated with large language models and verified kernel by kernel against the original. It is built to lower the barrier to experimentation, from new numerics and parameterizations to gradient-based calibration and hybrid physics–machine-learning components, while remaining close enough to the original so that what is developed in the shadow can be transferred back. In a 1958–2019 hindcast at 1◦ equivalent resolution with identical physics and forcing, the mean states of the JAX and Fortran versions differ from each other by two orders of magnitude less than either differs from observations, and the two runs agree for six decades in global temperature, salinity, heat content, and sea ice. The complete 1◦ configuration fits on a single GPU, a node of four GH200 superchips integrates ∼113 simulated years per wall-clock day, and meshes of up to 7.4 million surface vertices (∼5 km) scale to 128 GPUs. What limits the model is communication rather than arithmetic. What the shadow adds to the original is the gradient: a single reverse-mode pass through the full time loop returns the sensitivity of a model diagnostic to a parameter at every mesh vertex, verified against finite differences. To our knowledge, FESOM2-JAX is the first global ocean–sea-ice model of CMIP-class complexity written natively in a differentiable framework, and the first on an unstructured mesh.

1

Introduction

The Finite-volumE Sea ice–Ocean Model (FESOM2; Danilov et al., 2017) is a global ocean–sea-ice model built on an unstructured triangular mesh. Its cell-vertex finite-volume discretization lets the horizontal resolution vary smoothly over the globe, so a single configuration can resolve narrow straits, boundary currents, or a region of interest at eddying resolution while keeping the open ocean coarse. The model has been assessed in depth against observations and against its predecessor (Scholz et al., 2019, 2022), serves as the ocean component of the AWI climate model in upcoming CMIP7 (Streffing et al., 2022; Moon et al., 2025), part of the several kilometer-scale resolution modeling initiatives, including nextGEMS (Segura et al., 2025) and

1

Destination Earth Climate Digital Twin (Doblas-Reyes et al., 2026; Wedi et al., 2025). It scales efficiently to thousands of cores in its Fortran–MPI implementation (Koldunov et al., 2019a). Models of this class face two slow-moving but compounding pressures. The first is hardware and its toolchain. The computational capacity of current high-performance systems (HPC) is concentrated in GPUs, and accelerator toolchains for Fortran have matured later and less uniformly than their C++ and Python counterparts; there are several community’s routes around this, such as compiler directives (e.g. ICON; Giorgetta et al., 2022), domain-specific languages (Dahm et al., 2023; Adams et al., 2019), C++ performance-portability frameworks, including our own translation of FESOM2 to C++/Kokkos (Koldunov et al., 2026b), or new codes written GPU-native from the start (Ramadhan et al., 2020; Silvestri et al., 2025). The second pressure is people. The generation of scientists now entering the field is trained in Python and its machine-learning ecosystem. Much of model development is exploratory: trying a new advection scheme, a new parameterization, or a new calibration strategy. For such work, Python’s interactivity and libraries provide a better medium than Fortran production code, provided that the Python model remains faithful enough to the original for the results to carry over. Our response to these pressures, begun in Koldunov et al. (2026b) with an LLM-assisted translation of FESOM2 from Fortran through C to C++/Kokkos, is to develop and maintain code shadows. They are defined as projections of one model onto different technology stacks (C++/Kokkos for performance portability, Python for experimentation), but kept faithful to the original that casts them rather than allowed to drift into independent models. The Fortran production model remains the single source of truth (the ideal object, in the Platonic sense) so a shadow is neither a fork nor a successor; it is verified by running alongside the original on identical inputs and comparing outputs. This way developments made in the shadow, such as tuned parameter sets, trained machine-learning components, or validated numerical variants, can be carried back into the Fortran production model. This paper presents the Python shadow, FESOM2-JAX: a re-implementation of FESOM2 in python JAX (Bradbury et al., 2018). It is a NumPy-style codebase that a new researcher or student can read, run on a laptop or, unchanged, on 256 GPUs, and modify. This makes it a useful vehicle for teaching and a platform for experiments that are difficult to carry out in a Fortran production code. The clearest example of this category is the adjoint. Adjoint models, reverse-mode derivatives of the discretized dynamics (Thacker and Long, 1988; Errico, 1997), they underpin, for example, the ECCO ocean state estimates based on MITgcm ocean model (Stammer et al., 2002; Wunsch and Heimbach, 2007; Forget et al., 2015), variational data assimilation for NEMO and ROMS (Vidard et al., 2015; Moore et al., 2011), and adjoint sensitivity studies through ocean–sea-ice dynamics (Heimbach et al., 2010). Obtaining an adjoint from a Fortran model has, however, been challenging. It requires source-to-source transformation tools with tool-specific coding restrictions, hand-tuned checkpointing of the reverse sweep, and regeneration and re-verification whenever the forward model evolves (Giering and Kaminski, 1998; Heimbach et al., 2005; Utke et al., 2008; Griewank and Walther, 2000). As a result, most ocean models, including FESOM2, never acquired one, and parameter estimation has relied instead mostly on gradient-free methods (Sumata et al., 2019; Iglesias et al., 2013; Souza et al., 2020). In a differentiable framework, by contrast, the adjoint is derived by the framework from the same source that defines the forward model (Baydin et al., 2018; Paszke et al., 2019; Bezanson et al., 2017). It is exact for the discretization, remains consistent with the code as the code changes, and extends to whatever is added to the model, including neural-network components. This

2

is what makes hybrid physics–machine-learning models trainable through the dynamics (Rasp et al., 2018; Kochkov et al., 2024). Differentiable programming is therefore increasingly argued to be a foundation for the next generation of Earth-system models (Gelbrecht et al., 2023; Shen et al., 2023). This can be pursued either by rewriting models in differentiable frameworks, as we do here, or by differentiating existing compiled code through the compiler, the route pursued, for example, with Enzyme by the DJ4Earth project (Moses and Churavy, 2020; Moses et al., 2026). FESOM2-JAX belongs to a growing family of geophysical models written natively in Python and JAX. Veros (Häfner et al., 2018) translated the pyOM2 primitive-equation model into Python and, with a JAX backend, runs global structuredgrid configurations on GPUs at throughput competitive with its Fortran ancestor, including multi-node runs communicating through MPI (Häfner et al., 2021; Häfner and Vicentini, 2021). Veris (Gärtner et al., 2026) has recently added sea ice to the same ecosystem by re-implementing the MITgcm sea-ice model in JAX, matching hundreds of CPU cores on a single GPU. The translation of legacy Fortran into this ecosystem is also becoming more systematic. Large-language-model-assisted conversions to JAX have been demonstrated for a CESM photosynthesis module (Zhou et al., 2024) and, more systematically, for a 19,000-line land-surface model (Lahlou et al., 2026), with both gaining GPU execution and access to gradients. This is the same class of tooling with which the code-shadow lineage behind the present work was built (Koldunov et al., 2026b). Differentiable solvers of reduced complexity, including quasi-geostrophic models (Yan et al., 2025) and engineering CFD (Kochkov et al., 2021; Bezgin et al., 2023), have become testbeds for subgrid closures trained “online”, through the solver, which outperform the same closures fitted offline (Frezat et al., 2022). In the atmosphere, NeuralGCM (Kochkov et al., 2024) couples a differentiable spectral dynamical core in JAX to a learned physics package and trains the hybrid model end to end against reanalysis. Within this family, FESOM2-JAX is, to our knowledge, the first global ocean–sea-ice model of CMIPclass complexity: the full parameterization suite of such an ocean model and its dynamic–thermodynamic sea ice in a single differentiable time loop. It is also the first built on an unstructured mesh. The body of the paper is organized around the three requirements that FESOM2-JAX must meet for a code shadow to be useful: reproducible experimentation, scientific fidelity, and computational performance. 1. Reproducible experimentation (Sect. 3). A simulation is specified by a fully declarative YAML-based run configuration. The model supports restarts and basic diagnostics and provides all the functionality needed to conduct multi-decadal or high-resolution integrations within conventional Earth-system model data workflows. A continuous-integration ensures that the default configuration remains bit-reproducible across releases. 2. Scientific fidelity (Sect. 4). In a 1958–2019 hindcast on a mesh with an equivalent resolution of 1◦ , configured identically to a Fortran run, the JAX model closely reproduces the Fortran climatology. The comparison is statistical and includes climatological bias maps, multi-decadal drift, and sea-ice metrics. 3. Computational performance (Sect. 5). The model runs on meshes ranging from ∼ 105 to ∼ 7 × 106 surface vertices and on hardware ranging from a laptop CPU to 256 GPUs. It is fast in absolute terms but remains communication-bound. Section 2 first describes the model, its numerics, and the procedure by which the port was verified; Sects. 3–5 then address the three requirements in turn. Section 6 demonstrates the differentiability of the full time loop. The applications enabled by this 3

capability, including gradient-based calibration, hybrid physics–machine-learning modelling, and the transfer of their products back into the Fortran original, are the subject of ongoing work and lie beyond the scope of this model-description paper.

2

Model and numerics

FESOM2-JAX solves the same equations as the Fortran version of FESOM2, on the same geometry, with the same discretization and the same production physics. The two codes therefore differ in implementation, not in formulation, and everything described below is common to both unless explicitly noted. The port covers the components required for the forced ocean– sea-ice configuration of FESOM2; the components that are not ported are listed in Sect. 2.1. Table 1 summarizes the ported components and the settings used in this study, and Appendix A describes them one by one. The port itself was carried out kernel by kernel, following the LLM-agentic workflow described by Koldunov et al. (2026b); in the present work, the harness was Claude Code and the language models used were Opus 4.8 and Fable. 2.1

FESOM2 and the scope of the port

FESOM2 solves the hydrostatic, Boussinesq primitive equations on an unstructured triangular mesh with a cell-vertex finitevolume discretization. Scalar fields, including temperature, salinity, pressure, and sea-surface height, are located at mesh vertices and advanced on the control volumes, while horizontal velocity is located at triangle centroids. The vertical coordinate is the arbitrary Lagrangian–Eulerian (ALE) coordinate, used throughout this study in its z ⋆ form (Adcroft and Campin, 2004). Sea ice is a single-class dynamic–thermodynamic model, based on the Finite-Element Sea Ice Model (FESIM; Danilov et al., 2015) and discretized at the same mesh vertices as the ocean scalars. The discretization and numerics are described by Danilov et al. (2017), and further model elements, together with their assessment against observations and against the predecessor model, by Scholz et al. (2019, 2022). We do not reproduce that material here. Appendix A gives the component-by-component description for readers who want it without consulting those papers. The scope of the present model is the forced ocean–sea-ice problem. FESOM2-JAX is driven by prescribed atmospheric forcing and has no tides, icebergs, ocean biogeochemistry, ice-shelf cavities, floating ice, or partial bottom cells. Among the components that are ported, several of FESOM2’s alternatives are retained and are chosen from the configuration file: the ALE coordinate is ported in its z ⋆ form, with the linear free surface as the alternative surface treatment, while FESOM2’s z-level mode is not ported; vertical mixing can be run with the turbulent-kinetic-energy closure of Gaspar et al. (1990) as implemented in the CVMix library (cvmix-TKE; Griffies et al., 2015; Van Roekel et al., 2018), with the K-profile parameterization (Large et al., 1994), or with the Richardson-number scheme of Pacanowski and Philander (1981); and the sea-ice rheology can be either mEVP or standard EVP. Atmospheric forcing is the one component in which the port is narrower than the original, in that the JRA55-do reader (Tsujino et al., 2018) is its only forcing implementation. These choices define the boundaries of a first model-description paper, not limitations of the approach.

4

2.2

The JAX implementation

The implementation is ordinary python JAX array code (Bradbury et al., 2018), using JAX version 10.1. Three choices shape it: the timestep is a pure function of the model state and compiles to a single program, differentiability is treated as a property of every kernel rather than as a later addition, and one source runs unchanged on any device count. 2.2.1

One pure step, one compiled program.

The model state, including every prognostic field, the Adams–Bashforth history, and the warm-start vector of the CG solver, is stored as a single immutable pytree, and one timestep is a pure function of that state. Fields are stored as dense vertex-by-level arrays. The number of wet layers varies from column to column with the bottom topography, and that variation is carried by masking the cells below the bottom rather than by storing only the wet ones, so every array has a fixed shape and the XLA (Accelerated Linear Algebra) compiler can fuse the full timestep into one compiled program. The time loop is a single jax.lax.scan over this function, compiled once and iterated without returning to the Python interpreter. A long campaign is then a chain of such scans, reusing the compiled executable across job boundaries. Configuration is static at compile time: an inactive scheme, for example KPP when TKE is selected, is absent from the compiled program rather than branched over at run time. An all-on run therefore incurs no cost for options it does not use. All state and arithmetic are in float64, matching the Fortran model; this is required both for multi-decadal conservation and for the verification procedure of Sect. 2.4. 2.2.2

Differentiability by construction.

Reverse-mode differentiability is treated as a correctness property of every kernel, rather than as a feature added afterwards. Three patterns recur. First, iterative solvers with data-dependent stopping are avoided in favour of fixed iteration counts: five passes of the bulk-flux stability iteration, five Newton steps in the ice thermodynamics, and 120 mEVP iterations. This makes the backward sweep a well-defined unroll. The one exception is the CG solver for the free surface, which is wrapped in JAX’s custom_linear_solve so that the forward pass reproduces the reference’s early-stopped iterate, while the backward pass differentiates the linear system implicitly instead of unrolling the iteration. Second, non-smooth primitives are given finite gradients at their switching points: guarded square roots and powers at zero, clamped denominators, and identity rows in the tridiagonal solves for dry cells so that no gradient leaks below the bottom. Third, the memory of the reverse sweep over long integrations is contained by nested jax.checkpoint blocks around groups of timesteps (Sect. 6). 2.2.3

One source, any device count.

Distributed execution uses jax.shard_map over a one-dimensional device axis. The domain decomposition is FESOM2’s own: each device holds the entities it owns, followed by a halo rim of entities it reads but does not own, padded to a common size so that the per-device arrays keep a fixed shape. The layout is the familiar overlap arrangement of a Fortran ocean model. The difference lies in the exchange itself, which is not a call into a message-passing library but an operation written in the

5

same language as the rest of the model and compiled together with it. We implemented four transports, described in Sect. 2.3. Global sums are psum reductions. The same source also defines the distributed adjoint. In a Fortran model the adjoint of the exchange is written and maintained by hand. This holds whether the surrounding adjoint is generated by source transformation, as for the MIT general circulation model (Heimbach et al., 2005), in which case the exchange is a side-effecting subroutine wrapping message-passing calls that the tool cannot traverse, or hand-coded throughout, as for ROMS and NEMO (Moore et al., 2011; Vidard et al., 2015); the handling of parallelisation is among the reasons given for the latter choice (Vidard et al., 2015). The underlying mathematics is the same in every case: an exchange is linear in the field it moves, so its adjoint is its transpose, in which the adjoint values of the halo copies are returned to the owner and summed there. In the present case the exchange is free of side effects and is built from primitives whose transposes the framework already defines, so that this transpose is obtained by composition rather than by hand. This reduces the amount that must be verified but does not by itself establish correctness. We verify the distributed adjoint by requiring the sharded gradient to match the single-device gradient, which is itself checked against finite differences (Sect. 2.4). The broadcast, padded, and coloured transports carry that adjoint; the ragged all-to-all is used forward-only, because we found its reverse-mode rule to be defective in JAX 10.1. The step body is unchanged between a laptop CPU, a single GPU, and 256 GPUs. Running the sharded code path on one device reproduces the single-device path bit for bit, and this identity is the invariant by which we verify the distributed implementation. Restarts and output are written in the mesh’s global vertex numbering rather than in any device’s local order, so that the files a run produces do not depend on the number of devices used: a restart written on four GPUs can be read on 256, and output from either can be compared directly (Sect. 3). 2.3

The halo exchange, four ways

Because the halo exchange lies on the critical path of every stencil, we implemented it four ways; Fig. 1 shows what each of them puts on the wire. They differ along two axes. The first is whether the adjoint is exact. The all-gather broadcast, the padded all-to-all, and the coloured ppermute all differentiate correctly, because the adjoint of a halo read is an additive scatter back to the owner. By contrast, lax.ragged_all_to_all, which ships the least data, has a defective reverse-mode rule in JAX 10.1 and is usable forward-only. The second axis is how the per-device wire volume grows with the device count P . The all-gather moves the whole field, however little of it is halo. The ragged and padded all-to-alls carry one message slot for every device, whether or not that device owns any of the halo, so their cost rises as the partition is spread over more of them; the padded variant also pads every slot to a common size, so the devices that own none of the halo send full-sized messages of zeros. The coloured transport instead exchanges with one partner at a time, in K rounds coloured so that no device is asked to talk to two partners at once, and its cost is bounded by the largest number of neighbours any device has, 6 to 14 on the meshes used here. That bound does not grow with P , because a spatial partition keeps a bounded number of neighbours at any scale. All but the ragged transport also run on the CPU backend. Section 5.7 measures the four against each other.

6

(b) all-gather broadcast

(c) ragged all-to-all

(a) one device's share of the mesh

11 2 12 3 6 14

11

3

6

owned

12

2

9 peers with no halo here

owned vertices

9 peers with no halo here

the 30 entries it needs and no more, but a slot for every one of the 15 peers; 9 slots are zero-length and carry nothing

1138 entries received per step

30 entries received per step

(d) padded all-to-all

(e) coloured ppermute d11

14

P = 16 devices. This one owns 62 vertices and reads 30 more, from 6 of its 15 peers. Arrow colours are the rounds of (e).

11 2 12 3 6 14

a whole subdomain from every peer, 75 entries each, 7× off this scale; 30 of them are used

11 2 12 3 6 14

9 peers with no halo here

1

d2

2

d12

3

d3

d6

d14

4

5

6

every message padded to the same 8 entries, round (one ppermute each), in sequence so those 9 peers send full messages of zeros 6 messages, one per neighbour, against the 15 slots of (b) to (d); 6 is set by the colouring 120 entries received per step 30 entries received per step halo vertices and their messages one colour per ppermute round, in (a) and (e) on the wire, not needed

Figure 1. What each halo transport sends, for one device of a schematic partition. (a) A small unstructured mesh cut into P = 16 subdomains. The highlighted device owns 62 vertices (blue) and its stencils read 30 more that it does not own (orange): its halo. Six of the other 15 devices own part of that halo and are numbered; the remaining nine own none of it. (b to d) What the highlighted device receives per step under each of the three exchanges, one bar per device, all panels on the same vertical scale. The all-gather delivers every other device’s whole subdomain, which runs off the top of the scale, and 30 of the 1138 values received are used. The ragged all-to-all delivers exactly the halo entries each device owns and nothing more, but still carries a slot for all 15; the nine drawn as dashes are empty and move no data. This is the transport with no usable adjoint. The padded all-to-all delivers the same entries but requires every slot to be the same size, so those nine devices send full-sized messages of zeros that do cross the wire. (e) The coloured ppermute is a schedule rather than a single exchange, so its bars are rounds rather than devices: the halo arrives in K rounds, and in each round every device exchanges with exactly one partner. The arrows in (a) are coloured by round, and each bar is labelled with the partner it serves. It sends six messages against fifteen slots, and K is set by the number of neighbours rather than by the device count, which is what keeps its cost flat as the partition is spread wider; but the rounds run one after another, so it pays K start-up costs where a single exchange pays one. The totals are exact for this partition; on the meshes of Sect. 2.5 a device has 6 to 14 neighbours while P reaches 256, which widens the same gaps. Table 4 gives the measured per-step times.

2.4

Verification of the port

A port of this size is only as credible as its verification, so we describe the verification procedure and the reference from which the port was made. The proximate reference for FESOM2-JAX was not the Fortran source directly, but the instrumented serial

7

C translation of FESOM2 produced as the intermediate of the Fortran-to-Kokkos port of Koldunov et al. (2026b): a kernelby-kernel mirror of the Fortran code that writes every intermediate field of every substep to disk. For the present work, that C reference was extended beyond the configuration of the Kokkos study with the components used in this paper: the z ⋆ coordinate, the cvmix-TKE closure, and the mEVP sea-ice rheology. Each extension was itself a literal translation of the corresponding Fortran routine, validated against instrumented Fortran output before being adopted as a baseline. The unit of porting was then the kernel, corresponding to one Fortran subroutine’s worth of physics, and every JAX kernel was gated against the C reference’s dumps. JAX does not reproduce the loop order of the reference, so agreement cannot be exact, and the tolerance is therefore set by what the kernel does. Kernels that work point by point, or that only read values at neighbouring points, are required to agree to about 10−15 relative, which is the level of rounding alone. Kernels that accumulate contributions, such as the assembly of edge fluxes onto the control volumes and the global sums, are added up in a different order by JAX and are required to agree to about 10−12 . Multi-step replay tests then verify the state threading that single-step tests cannot see, including the Adams–Bashforth history, the CG warm start, and the layer-thickness commit, and the sharded code path must reproduce the dense path exactly on one device and to reassociation tolerance on many. Differentiability is tested with the same rigour as the forward model: gradients are checked against finite differences kernel by kernel, at states away from the switching points at which the model is not differentiable, and through short full-model integrations. The suite comprises 78 test modules (roughly 640 tests), of which 41 modules exercise gradients; the continuousintegration gate of Sect. 3 additionally pins the shipped configuration to a stored reference trajectory bit for bit. Finally, the chain closes back on the original: the hindcast of Sect. 4 compares the assembled model against Fortran FESOM2 itself, rather than against the C reference the port was built from, so an infidelity introduced at either translation step, Fortran to C or C to JAX, would appear in this comparison. The same procedure localized the freshwater-budget defect of Sect. 4. Because every kernel had already been verified on its own, a drift that appears only over climate timescales cannot originate inside a kernel and must instead arise where two kernels are joined. Checking which of the conservation budgets fails to close then locates the defect. 2.5

Meshes

Five meshes are used in this paper, and they are the same community FESOM2 meshes the Fortran model runs on, each introduced in the study cited beside it in Table 2; the geometry tabulated there is measured from the mesh files the model reads. They span two orders of magnitude in size, from a 1◦ -equivalent configuration to a mesh with 7.4 million surface vertices, and they follow different design principles. FORCA20 is a Mercator mesh, on which resolution is set by latitude and increases gradually towards the poles. DARS follows the local Rossby radius (Sein et al., 2017) and the energetically active regions. NG5 is close to uniform at 5 km over most of the globe and is coarsened only where the flow is less energetic. fArc is CORE2 refined in a single region, the 4.5 km Arctic. One property of the hierarchy matters for the rest of the paper. The vertical extent grows with the horizontal one, from 47 layers on CORE2 to 69 on NG5, so the cost of a mesh grows faster than its surface-vertex count: NG5 has 58 times CORE2’s vertices but 86 times its vertex-levels, the product of vertices and layers that the model actually computes on. CORE2 is the 8

Table 1. Model components ported in FESOM2-JAX and the schemes used in this study. The Gent–McWilliams/Redi parameterization is active in the CORE2 hindcast and switched off in the configurations used for the performance measurements; all other components are common to every configuration. Component

Scheme / setting

Reference

Vertical coordinate

ALE, z ⋆ (linear free surface also ported)

Adcroft and Campin (2004)

Equation of state

polynomial, in-situ density

Jackett and McDougall (1995)

Pressure gradient

density-Jacobian

Free surface

Shchepetkin and McWilliams (2003) −5

semi-implicit (α = θ = 1), preconditioned CG to 10

,

warm-started

Danilov et al. (2017); Marshall et al. (1997)

Time stepping

Adams–Bashforth 2 (ϵ = 0.1); implicit vertical mixing

Tracer advection

MUSCL-type–FCT; 3rd-order horizontal, 4th-order ver-

Danilov et al. (2017); Zalesak (1979)

tical Vertical mixing

cvmix-TKE, ck = 0.1, cϵ = 0.7 (KPP, PP also ported)

Gaspar et al. (1990); Griffies et al. (2015); Van Roekel et al. (2018)

Mesoscale eddies

GM (streamfunction BVP) + Redi, coarse mesh only

Horizontal viscosity

biharmonic filter, flow-aware coefficient

Danilov et al. (2017); Juricke et al. (2020)

Sea-ice dynamics

mEVP, α = β = 250, 120 iterations (EVP also ported)

Bouillon et al. (2013); Kimmritz et al.

Gent and McWilliams (1990); Ferrari et al. (2010)

(2015) Sea-ice thermodynamics

zero-layer, 7 thickness classes

Semtner (1976); Hibler (1979)

Sea-ice advection

Taylor–Galerkin FE-FCT

Löhner et al. (1987)

Atmospheric forcing

NCAR bulk formulae, JRA55-do

Large and Yeager (2009); Tsujino et al. (2018)

Surface salinity

relaxation to monthly climatology, 10 m per 60 days

Precision

float64 throughout

Time integration

single jax.lax.scan; shard_map over devices

mesh of the hindcast in Sect. 4 and the only mesh on which the GM/Redi parameterization is active in this study; FORCA20 carries the eddying comparison of Sect. 4; and all five are used in the performance measurements of Sect. 5. Appendix A7 notes how far the nominal resolution of Table 2 departs from the area-weighted one.

3

Usability and reproducibility

For the Python JAX shadow, usability is a design requirement rather than a byproduct. This section describes the design choices intended to make FESOM2-JAX a practical tool rather than a research prototype: single-file configuration, restarts,

9

Table 2. The meshes used in this study. Layers is the number of vertical layers and vertex-levels the product of the two, which is the number of entries in the model’s dense state arrays: because the varying number of wet layers is carried as masks rather than by storing only the wet cells (Sect. 2.2), every timestep advances all of those entries, and not only the 48–80 % of them that are wet ocean, so the full product is the work a timestep performs. Resolution is the nominal resolution each mesh is designed and named for, with the region it applies to where the mesh is regionally refined; the vertex, layer and vertex-level counts are measured from the mesh files the model reads. ∆t is the production timestep, at which all throughput in Sect. 5 is reported. Measured per-step cost on each mesh is given in Table 3. surface

∆t

vertex-levels

Mesh

vertices

layers

(106 )

resolution

(s)

CORE2

126 858

47

6.0

1◦ equivalent

1800

Wang et al. (2014)

fArc

638 387

47

30.0

4.5 km, Arctic

1200

Wang et al. (2018)

FORCA20

2 127 871

69

146.8

1/5◦ equivalent

240

Koldunov et al. (2026a)

DARS

3 160 340

56

177.0

∼10 km

240

Streffing et al. (2022)

NG5

7 402 886

69

510.8

∼5 km

240

Rackow et al. (2025)

Reference

data output, and an explicitly scoped reproducibility guarantee. In the resulting workflow, a single YAML RunConfig drives one invocation: load a restart or an initial condition, integrate N steps as a jax.lax.scan, stream online per-period means and subsampled snapshots, and write restarts. A scheduler dependency then chains these invocations across batch jobs to assemble a long integration. 3.1

One run, one file

A simulation is specified by a single YAML RunConfig: the mesh and its partition, the physics, the timestep, the forcing, and the output and restart targets. Physics components are enabled by their presence and configured by their contents. Mutually exclusive choices, such as the two vertical-mixing schemes, are rejected at load time rather than at run time. Every key is optional, and an absent key falls back to a default. This property is enforced by a regression test, so that adding configurability does not silently change an existing result. The abridged listing below is the all-on CORE2 configuration used for the hindcast of Sect. 4: ale: {}

# z* vertical coordinate (absent => linear free surface)

tke: {}

# cvmix-TKE vertical mixing (exclusive with kpp)

ice: {whichEVP: 1}

# sea ice, modified-EVP rheology

gm:

# Gent-McWilliams + Redi (coarse mesh)

{}

visc: {}

# flow-aware viscosity, default coefficients

tracer: {}

# FCT tracer advection

dt: 1800.0

# timestep [s]

10

mesh: data/mesh_core2 partition: dist_4

# device count of the on-disk partition

forcing: {kind: core2, start_year: 1958}

# JRA55-do, real calendar

snapshot_every: 48

# steps (~daily); 0 = off

checkpoint_every: 1440

# rolling-restart cadence (~monthly)

restart_in:

# null => cold start from climatology

null

restart_out: runs/core2/restart restart_archive_out: runs/core2/restart_archive restart_archive_period: year duration: 10yr

# immutable archival restarts

# or an explicit step count

One invocation performs a single task: it loads an initial condition or restart, integrates for a requested number of steps or a specified duration (“10yr”, “3mo”), streams diagnostics, writes restarts, and exits. Multi-job integrations, such as the 62-year hindcast below, are assembled outside the model by an ordinary batch-scheduler dependency chain, so no bespoke driver is required. 3.2

Restarts and output

The model writes two restart streams for different purposes. A rolling restart, written every checkpoint_every steps into a single directory, supports crash recovery and the handoff between chained batch jobs, while an archival stream writes an immutable restart at calendar boundaries, yearly here and monthly for the high-resolution runs, under a unique name, fesom.<YYYY>.<DDD>.<SSSSS>, together with a restart.latest pointer. Diagnostics are computed online as percalendar-period means, with daily and monthly streams whose field lists are configurable, and are written alongside timesubsampled snapshots in a single pass and with constant memory use. Restarts and output alike are written in the mesh’s global vertex numbering rather than in any device’s local order, so that a run checkpointed on N GPUs can be resumed on any number M , and a field can be analysed without knowing the decomposition it was produced on. Output is written as plain Zarr stores with vertex coordinates and a CF-convention time axis on a real calendar, so that a store can be opened directly in xarray or in an unstructured-mesh viewer. 3.3

Reproducibility and its scope

The shipped baseline configuration reproduces a stored reference trajectory bit for bit, and a continuous-integration test fails on any deviation. This gate guarantees that a given release computes a specific, reproducible result, and it is one reason the model is kept in float64 throughout. The guarantee holds on a single device. Multi-GPU runs are reproducible only to round-off, because the reduction order of the collectives is not associative (Sect. 4 accounts for this in the comparison methodology). The model, a user guide, and the tutorial notebook are archived as a tagged release (see Code availability).

11

4

Fidelity: a CORE2 1958–2019 hindcast against Fortran FESOM2

4.1

Experimental design

We run FESOM2-JAX and Fortran FESOM2 in the same configuration: the CORE2 mesh (Table 2), 1958–2019 forced by the JRA55-do reanalysis (Tsujino et al., 2018), a cold start from the PHC3.0 winter climatology (Steele et al., 2001), the z ⋆ coordinate with cvmix-TKE vertical mixing, mEVP sea ice and GM/Redi (Table 1), and ∆t = 1800 s. The comparison is against Fortran FESOM2, not against the C reference from which the port was built (Sect. 2.4). Both translation steps are therefore under test: an error introduced in the Fortran-to-C step would be invisible in a comparison against C, but appears here. The ocean is chaotic and the two runs are not bit-identical, so their internal variability diverges however faithful the port is (Sect. 3). We therefore compare statistics rather than instantaneous fields: the time-mean state, the biases against observations, the multi-decadal drift and the sea-ice cycle. Monthly-mean climatologies are compared over 1980–2009 and integrals over the full hindcast. 4.2

Mean state

Figure 2 shows the 1980–2009 annual-mean surface bias of each model against PHC3.0. The two bias patterns are visually indistinguishable, with root-mean-square (RMS) errors against observations of 0.61 ◦ C in both runs for sea-surface temperature (SST) and 0.38 for salinity (SSS). Both models carry the familiar CORE2 mean-state errors: the warm Gulf Stream-separation bias, the cold subpolar gyres, and the Southern Ocean pattern. The JAX−Fortran difference in the third column is two orders of magnitude below either model’s bias against observations, at 0.004 ◦ C and 0.002 RMS. It is not spatially uniform but concentrated along the western boundary currents, the Southern Ocean fronts and the marginal ice zones. That structure may be worth investigating, but at this amplitude it has no practical consequence. Figure 3 repeats the comparison in the vertical, as global zonal-mean sections. The bias sections are again indistinguishable, with RMS biases of 0.379 ◦ C in both runs for temperature and 0.077 for salinity. The JAX−Fortran difference is 1.2 × 10−3 ◦ C and 1.4 × 10−4 RMS, and is largest in the upper few hundred metres at high latitudes, where the two runs diverge fastest through the mixed layer and the ice edge; the deep ocean contributes almost nothing. 4.3

Multidecadal drift

Figure 4 shows the difference between the two models’ global integrals over the hindcast. Both models cool from the PHC initial state and settle without runaway drift, the volume-mean temperature falling from 3.636 ◦ C to 3.619 ◦ C in both runs, a change of −1.7 × 10−2 ◦ C over six decades. The two codes differ by +8.8 × 10−5 ◦ C, and that difference fluctuates about zero rather than growing. The 0–700 m mean diverges faster than the deep ocean and carries the larger residual, but stays within

12

CORE2 annual-mean surface bias vs PHC3.0 (1980-2009, 360 months) JAX Fortran

0.02

RMSE = 0.606

1 2 3 4

RMSE = 0.380

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0

0.00

RMSE = 0.380

JAX Fortran [°C]

0.04

0.02 0.04 RMS = 0.004

0.02 0.01

SSS bias [psu]

RMSE = 0.605

4 3 2 1 0

0.00

JAX Fortran [psu]

Fortran PHC

SST bias [°C]

JAX PHC

0.01 0.02 RMS = 0.002

Figure 2. CORE2 1980–2009 annual-mean surface bias against PHC3.0. Rows: SST (top), SSS (bottom). Columns: JAX−PHC,

CORE2 global zonal-mean T/S bias vs PHC3.0 (1980-2009) Fortran PHC

0

depth [m]

0.02

T bias [°C]

1

0.00

1

0.02

2

4000 RMS bias = 0.379

RMS bias = 0.379

0.04

RMS = 1.2×10 3

3 1.0

0.010

0.5

0.005

0.0

2000

0.000

0.5

4000 6000

0.04

2

2000

6000 0 50 100 250 500 1000

JAX Fortran

3

RMS bias = 0.077

75

50

25

RMS bias = 0.077

0 25 latitude [°N]

50

75

75

50

25

1.0 0 25 latitude [°N]

50

75

JAX Fortran [°C]

0 50 100 250 500 1000

JAX Fortran [psu]

JAX PHC

S bias [psu]

depth [m]

Fortran−PHC, and JAX−Fortran. Note the much tighter colour bar for the difference column (±0.05 ◦ C, ±0.025). RMSE is annotated.

0.005 RMS = 1.4×10 4

75

50

0.010 25

0 25 latitude [°N]

50

75

Figure 3. CORE2 global zonal-mean temperature (top) and salinity (bottom) bias against PHC3.0 over 1980–2009, on a sqrt-warped depth axis. Columns: JAX−PHC, Fortran−PHC, and JAX−Fortran. Note the much tighter colour bar for the difference column (±0.05 ◦ C, ±0.01). The annotated RMS is that of the zonal-mean bias, so it is smaller than a full three-dimensional RMS.

13

5.4 × 10−3 ◦ C of the Fortran run even in its most extreme single month. Ocean heat content stays within 0.96 ZJ throughout and ends +0.48 ZJ apart. Volume-mean salinity changes by +5 × 10−5 on the practical salinity scale over 1958–2019, against +2 × 10−5 for Fortran. Salinity is the one integral whose JAX−Fortran residual has a consistent sign, growing smoothly to +2.9 × 10−5 by the end of the hindcast rather than fluctuating about zero as the temperature difference does. The vertical structure of the difference, averaged over 1980–2009 in panel (d), reaches 1.3 × 10−3 ◦ C at 255 m and 3.8 × 10−4 at 15 m and decays towards the abyss, with the same upper-ocean confinement as the sections of Fig. 3. A comparison of this length is what makes such residuals visible at all: an earlier version of the port leaked freshwater through the global freshwater-flux normalization of the z ⋆ coordinate, a defect that showed up only as a slow salinity drift over decades and that a shorter run would have missed (Sect. 2.4). 4.4

Sea ice

Figure 5 shows the seasonal cycle of the JAX−Fortran ice-area difference. The largest monthly difference is 0.06 ×1012 m2 in the Arctic (June) and 0.16 in the Antarctic (November), or 0.5 % and 1.1 % of the ice area in those months; the larger Antarctic value falls near the summer minimum, when the ice area itself is small. At the four seasonal extremes the two models differ by at most 0.07 ×1012 m2 , and they share the same errors against OSI-SAF (EUMETSAT Ocean and Sea Ice Satellite Application Facility, 2022): both carry too much Arctic ice through the cycle, with a seasonal RMS difference of 1.62 ×1012 m2 , and both reproduce the Southern-Hemisphere cycle to 0.78 ×1012 m2 with a low summer minimum. The agreement holds at the minima as well as the maxima. 4.5

Towards high resolution

The test above is at coarse resolution, but the unstructured mesh is what makes FESOM2 attractive at high resolution, and the JAX port runs the same physics on eddy-permitting meshes without modification. Figure 6 compares the Gulf Stream region in the second year of a hindcast on FORCA20 (Table 2). Both implementations produce an eddying ocean of the same character, with a Gulf Stream separating at Cape Hatteras, meanders and rings along the front, and comparable jet sharpness and surface speed. A quantitative eddying evaluation, including spectra, eddy statistics and boundary-current transports, is left to future work.

5

Performance and scalability

FESOM2-JAX runs the whole mesh hierarchy of Table 2 at useful speed, on device counts from one GPU to 256, and from the same source on two different accelerators. What limits it is the shape of the scaling curve rather than the absolute throughput: the model is communication-bound, and it becomes so at partition sizes at which the Fortran original is still scaling. We quantify both parts of that statement below. Figure 7 shows strong scaling and throughput across the hierarchy, Fig. 8 reduces those curves to the single quantity that governs them, and Table 3 gives representative production numbers. 14

CORE2 global drift

JAX minus matched Fortran (1958-2019, 744 months)

0 700 m

0.5 0.0 0.5 1.0

1.0 JAX Fortran [ZJ]

full depth

JAX Fortran [10 5 psu]

1.0

1960 1970 1980 1990 2000 2010 2020 year (c) ocean heat content

0.5

depth [m]

JAX Fortran [10 3 °C]

(a) global mean ocean temperature

0.0 0.5 1.0

3 2 1 0 1 2 3

0 100 250 500 1000 2000

(b) global mean ocean salinity

1960 1970 1980 1990 2000 2010 2020 year (d) mean-state profile difference, 1980 2009

4000 1960 1970 1980 1990 2000 2010 2020 year

6000

T S

1.5 0.4

1.0

0.5 0.0 0.5 T [10 3 °C] 0.2

0.0 0.2 3 S [10 psu]

1.0

1.5 0.4

Figure 4. Global drift over the CORE2 hindcast (1958–2019), shown throughout as JAX minus the matched Fortran run. (a) volume-mean temperature, full depth and 0–700 m; (b) volume-mean salinity; (c) ocean heat content; all three as annual means. (d) the JAX−Fortran temperature and salinity profiles averaged over 1980–2009, on a sqrt-warped depth axis, each with its own x-axis. Panel (c) is panel (a)’s full-depth curve rescaled, since the heat-content diagnostic uses a constant ρ0 cp and the fixed mesh volume.

5.1

What is measured, and how

All numbers are for the complete ocean–sea-ice model in each mesh’s production configuration, with the same physics as the hindcasts of Sect. 4: z ⋆ vertical coordinate, cvmix-TKE mixing, mEVP sea ice with its 120-iteration subcycle, and, on CORE2, the only mesh coarse enough to need it, GM/Redi. The model is driven by real JRA55-do forcing from the PHC

15

JAX Fortran ice area [1012 m2]

CORE2 sea-ice area

JAX minus matched Fortran (1980-2009; obs reference OSI-SAF CDR v3p0 (ICDC)) (a) Northern Hemisphere

(b) Southern Hemisphere

J F M A M J J A S O N D month

J F M A M J J A S O N D month

0.2 0.1 0.0 0.1 0.2

Figure 5. CORE2 sea ice (1980–2009): the seasonal cycle of the ice-area difference, JAX minus the matched Fortran run, for the Northern (a) and Southern (b) Hemisphere.

JAX

SST

FORCA20 Gulf Stream, daily mean 1959-01-05 (hindcast year 2) Fortran

SST 25 20 SST [°C]

15 10 5 0

JAX

surface speed

Fortran

surface speed 1.2 1.0 speed [m s 1]

0.8 0.6 0.4 0.2 0.0

Figure 6. Daily-mean SST (top) and surface speed (bottom) in the Gulf Stream region on 5 January of hindcast year 2 on the FORCA20 mesh (2.13 × 106 surface vertices, ∼14 km local resolution): JAX (left) versus the matched Fortran run (right). The Fortran fields come from a one-month re-integration of the reference run’s January 1959, restarted from its end-of-1958 state with daily output enabled (the reference run stores yearly means). The two models decorrelate during the preceding year of eddying integration, so eddy positions differ by construction.

initial state, and throughput is always reported at the production timestep of Table 2. Measurements were made on two GPU systems: the DKRZ Levante GPU partition, with four NVIDIA A100-80 GB per node, and the JUPITER booster at the Jülich 16

Supercomputing Centre, with four NVIDIA GH200 Grace–Hopper superchips per node, 96 GB of HBM3 each, and NDR200 InfiniBand. Unless a number is labelled GH200, it was measured on A100. Each point runs on its fastest halo transport (Table 4). Timing follows a compile-once protocol: the step chain is compiled, invoked once to warm it, and the same executable is then invoked a second time and timed over 150 steps with device synchronization, twice per point, so that XLA compilation never contaminates the per-step numbers. Compiling the full step chain takes roughly half a minute on CORE2 and DARS and about two minutes for NG5 at 32 GPUs; the executable is cached on disk and reused across the chunks of a long run, so it is paid once per configuration rather than once per job. In every run the temperature and velocity fields remain free of NaNs over the timed window, so no timing comes from a diverged integration. The protocol has two departures. DARS is benchmarked at ∆t = 120 s rather than at the production 240 s, because a coldstarted benchmark window is not stable at the production step on every partition of this mesh; the per-step cost is timestepindependent, so the measured step time is unaffected. The DARS curve also starts at 16 GPUs: at 8 GPUs the productionphysics step exceeds a practical XLA compilation budget, an upstream compiler pathology at very large per-device partitions that is under investigation. The same configuration compiles in about a minute from 16 GPUs upward. 5.2

Absolute throughput

The model is fast in absolute terms, and the whole hierarchy is affordable. A single A100 runs the complete CORE2 configuration at 0.229 s per step (21 simulated years per wall-clock day, SYPD) within 11.2 GiB of device memory, so that configuration fits on one accelerator; a single node of four brings it to 0.085 s (58 SYPD) at 67 % parallel efficiency, which makes the 62-year hindcast of Sect. 4 a modest computation, and a node of four GH200s runs it at 0.044 s (112.8 SYPD). At the other end of the hierarchy, NG5 sustains 1.50 SYPD on 64 A100s and 2.79 SYPD on 64 GH200s. Table 3 gives the same numbers for every mesh. 5.3

One quantity governs the scaling

We measure the work a device performs in vertex-levels: one layer of the water column beneath one surface vertex, so a mesh of N surface vertices and L layers has N × L of them, from 6 million on CORE2 to 511 million on NG5 (Table 2). The model holds its state in dense arrays of that shape and advances every entry each timestep, including the masked ones below the sea floor, so vertex-levels per second is a rate of work comparable across meshes of different size and depth. Counted this way, the rate one GPU sustains depends on how many surface vertices it owns and on little else (Fig. 8a). Above 120 000 surface vertices per device the rate is 24 million vertex-levels per second on an A100 and 57 million on a GH200, and the five meshes agree to ±11 % (21–27 million) although they differ by a factor of 58 in size, run from a 1◦ -equivalent configuration to ∼5 km, and carry between 47 and 69 layers. Below that threshold the rate falls away steadily: half as many vertices gives 79 % of the plateau, and fewer than 30 000 vertices only 30 %. The halo is a rim around a two-dimensional subdomain, so the communication-to-work ratio grows as the inverse square root of the vertices per device and is independent of the layer count, which is why meshes of 47 and 69 layers fall on one curve.

17

(a) strong scaling (dotted = ideal 1/N)

(b) throughput, CORE2 & farc

time per step [s]

0.4 0.2 0.1 0.05 1

CORE2 (1°)

4

16

64

number of GPUs

256

fArc (4.5 km, Arctic)

(c) throughput, multi-million-vertex meshes

100

SYPD (simulated yr / wall day)

SYPD (simulated yr / wall day)

0.8

80 60 40 20 0

1

DARS (10 km)

2

4

8

16

number of GPUs

32

FORCA20 (1/5°)

5 4 3 2 1 0

64

16

32

64

128

number of GPUs

NG5 (5 km)

A100 (Levante)

256

GH200 (JUPITER)

Figure 7. Strong scaling of FESOM2-JAX for the full ocean+sea-ice model, using the same source code on two GPU systems: A100 (DKRZ Levante; solid, circles) and GH200 (JUPITER booster; dash-dotted, triangles). (a) Time per step versus the number of GPUs, per mesh; dotted lines are ideal 1/N scaling anchored at each mesh’s first A100 point. Every point uses its fastest halo transport (Table 4): on A100 the padded or ragged all-to-all at 32 GPUs and below and the coloured ppermute at 64 and above, on GH200 the coloured transport at every multi-node point. (b, c) Throughput in simulated years per wall-clock day at each mesh’s stable production timestep, split by mesh size so that the small meshes’ throughput (b) does not compress the multi-million-vertex meshes’ (c). The small CORE2 mesh saturates beyond one node on either machine while the larger meshes strong-scale sub-linearly. GH200 is roughly twice as fast per GPU wherever the per-device partition is large, and each mesh’s knee correspondingly moves to about half the device count. FORCA20 was measured on A100 only. Table 3. FESOM2-JAX at the production timestep on both GPU systems, run with the same source, JAX version, physics and protocol on each. FORCA20 was measured on A100 only. Speed-up is the A100 per-step time divided by the GH200 one at the same GPU count. Peak device memory is a property of the partition rather than of the accelerator (NG5 at 32 GPUs measures 48.5 GiB on GH200 against the 48.2 GiB tabulated here). A100 (Levante)

GH200 (JUPITER)

Mesh (#GPU)

vertices

s/step

SYPD

GPU-h/yr

GiB/GPU

s/step

SYPD

GPU-h/yr

speed-up

CORE2 (4)

127 k

0.085

57.7

2

3.1

0.044

112.8

1

1.96

CORE2 (8)

127 k

0.074

66.5

3

1.7

0.046

106.9

2

1.61

DARS (16)

3.16M

0.450

1.5

263

28.1

0.183

3.6

107

2.46

FORCA20 (16)

2.13M

0.389

1.7

227

21.6

–

–

–

–

NG5 (32)

7.40M

0.755

0.9

883

48.2

0.317

2.1

370

2.38

NG5 (64)

7.40M

0.438

1.5

1024

34.2

0.236

2.8

551

1.86

18

50 20

(b) the price of the same simulation, faster cost [GPU-h per simulated year]

throughput [M vertex-levels s 1 GPU 1]

(a) per-GPU throughput vs shard size GH200 kernel-bound, 57 A100 kernel-bound, 24

10 5 2

CORE2 (1°)

10

20

50

100

surface vertices per device [thousands]

fArc (4.5 km, Arctic)

DARS (10 km)

200

1000 100 10 1 1

FORCA20 (1/5°)

4

16

number of GPUs

NG5 (5 km)

64

A100 (Levante)

256

GH200 (JUPITER)

Figure 8. The same measurements as Fig. 7, reduced to what governs them. (a) The rate at which one GPU advances vertex-levels, against the number of the mesh’s surface vertices it owns, taken as the mesh’s vertex count divided by the device count (the FESOM2 partitions used here are balanced to within 1.2 % at every point plotted); grey guides mark the kernel-bound plateau. CORE2 lies above the common curve at small shards because it is the only mesh carrying GM/Redi, so that the same halo is overlapped with more arithmetic. (b) Accelerator hours per simulated year against device count.

A mesh therefore strong-scales while its shards are large and saturates once they are not, so the device count at which it stops improving grows with its size, and a faster accelerator raises the throughput at a given shard size without changing the shard size at which scaling ends. 5.4

Strong scaling

CORE2 exhausts its parallelism within two nodes: the second node still reduces the step time, but at 58 % doubling efficiency, and at 8 GPUs each device holds about 16 000 vertices, the small-shard end of Fig. 8a. The eddy-permitting meshes strongscale sub-linearly and keep improving out to the largest A100 partition we tested: the first doubling of each retains 75–86 % parallel efficiency, and DARS still returns 63 % on its second (32→64 GPUs). At 128 GPUs a DARS step is 0.19 s and an NG5 step 0.35 s, both still falling, while FORCA20, the smallest of the three, has flattened at 0.23 s. fArc, between these two regimes, gains a factor 1.85 from 4 to 16 GPUs. Whether these meshes keep scaling at high device counts also depends on the halo transport. With a single P -way all-to-all, DARS and FORCA20 turn over past 64 GPUs, because that transport’s per-device message count grows with the device count P . The coloured ppermute, whose cost is set by the mesh’s colouring rather than by P , removes the turnover. Section 5.7 shows how their ranking changes with scale.

19

5.5

The same code on a second machine

Because the model is ordinary JAX code, running it on a different accelerator required no change to the source: the JUPITER installation is a virtual environment built from the same jax[cuda12] wheels, which are available for the machine’s aarch64 hosts, with mesh exports regenerated on site and the same YAML run configurations, and needed no container image or vendorspecific build. We repeated the campaign there from a single GPU to 256 (dash-dotted curves in Fig. 7; the GH200 columns of Table 3). The gain is confined to throughput per GPU, where it is a factor of 2.3 at the kernel-bound end of Fig. 8a; it decays as the partitions shrink and communication becomes the limiting cost, and at the smallest partitions it inverts, fArc on 16 GPUs being slower on GH200 than on A100 (0.90×). These kernels are gather- and bandwidth-bound rather than FLOP-bound, so the ratio follows the memory bandwidth of the two devices, a factor of about two, far more closely than their arithmetic peak, a factor of about three and a half. The range over which the model scales is unchanged: every mesh turns over at roughly half the device count it does on A100, NG5 for instance saturating at 64 GPUs (0.236 s), unchanged at 128 (0.236 s) and slightly slower at 256 (0.248 s). This follows from the per-GPU gain: a device twice as fast completes a shard’s work in half the time while the communication per step is unchanged, so the compute–communication crossover arrives at half the shard count, and a faster machine is best used with fewer, larger shards. The multi-node points were measured on a busy shared system, where repeats within one allocation agree to 0.05 % but the spread between allocations is ∼ 5 %, so the turnover locations carry that uncertainty. 5.5.1

Two lessons about portability.

A source that transfers unchanged does not carry its tuning with it, in two respects. The first is that the ranking of the halo transports depends on the machine as well as on the mesh, and the Levante ranking does not transfer. On GH200, the padded all-to-all is slower than the coloured ppermute at every multi-node point we measured, by factors between 1.1 and 5.5, and its timings do not reproduce between allocations; the coloured transport reproduces and is fastest at every point beyond one node, so every multi-node GH200 point uses it. Contention from our own concurrent jobs, a degraded fabric (a direct bandwidth probe through the same communication library found it healthy), and partition imbalance were each tested and rejected as explanations, and the cause remains open. Because the four transports are interchangeable behind a run-time flag, adapting to the new machine amounted to selecting a different one. The second is that a code whose inner loops are written by a compiler also inherits the compiler’s judgement about how to write them, including its mistakes. At the NG5 partition into 128 devices, and at no other partition, XLA merged the global surface-flux balance sum into the sea-ice thermodynamics, which it compiles as a single block of 1425 operations, and switched that block to a slower code path, producing one 383 ms kernel and a 2.7× increase in step time. The effect was bracketed by unaffected points at 64 and 256 GPUs, reproduced on separate allocations, and present under all three transports. Forcing the compiler to keep the sum separate, which leaves the outputs bit-identical, removes it and reduces the step time from 0.637 to 0.236 s; that barrier is now the default and the point plotted in Fig. 7 is the corrected one. This kind of failure belongs to compiler-generated code, and is part of the price of the portability demonstrated above.

20

5.6

Where the time goes

Switching components off decomposes the CORE2 4-GPU production step. The bare ocean core, comprising dynamics, the free surface, tracers, and forcing, accounts for ∼ 58 ms per step. The mEVP sea ice adds ∼ 14 ms, almost entirely from the 120 subcycled rheology iterations, each of which performs its own halo exchange. The GM/Redi eddy parameterization adds ∼ 9 ms, the z ⋆ moving vertical coordinate ∼ 5 ms, and the TKE mixing closure only ∼ 3 ms. The single-GPU baseline prices the communication directly: relative to ideal scaling from 0.229 s, communication and duplicated halo work account for about 33 % of the measured 4-GPU step. Both point at the halo exchange, which lies on the critical path of every stencil and every sea-ice subcycle. 5.7

Which transport is fastest

The ranking of the four transports of Sect. 2.3 changes with scale (Table 4). On the small CORE2 mesh, the minimal-volume ragged and padded exchanges are fastest and the coloured transport is slowest: the SSH conjugate-gradient solve fires a few hundred small two-dimensional exchanges per step, these are latency-bound, and K sequential ppermute rounds pay K launch latencies where a single all-to-all pays one. On the largest configuration, NG5 at 64 GPUs, the order reverses. The three-dimensional tracer exchanges dominate, these are bandwidth-bound, and the coloured transport’s flat-in-P volume beats both the ragged all-to-all (by 8 %) and the padded one (by 37 %). DARS at 32 GPUs sits just before the crossover: ragged and padded tie at 0.35 s and the coloured transport remains slightly slower at 0.39 s. An exact adjoint therefore costs nothing: at both ends of the mesh range the fastest transport with an exact adjoint is faster than the forward-only ragged all-to-all (Table 4), and the large end is where gradients are scientifically useful. Choosing the transport per point recovers part of the communication cost but not all of it. The unstructured partition makes the exchanges irregular, since each rank communicates with a data-dependent set of neighbours, and their relative cost grows as the subdomains shrink, which is the effect Fig. 8a measures. Two improvements are already reflected in Fig. 7: the padded and coloured transports move each stencil’s halo in one statically shaped collective, and fusing the sea-ice subcycle’s velocity pair and the solver’s paired reductions removed a further ∼ 9 % of the CORE2 8-GPU step. 5.8

Optional accelerations

Two further optimizations are implemented as opt-in switches and are OFF in every number reported here, because each changes results within solver or floating-point tolerance rather than bit for bit. The first computes the per-step atmospheric forcing on the device, removing host-side interpolation and transfer. The second replaces the preconditioner of the sea-surfaceheight solve with a degree-3 Chebyshev polynomial, which reaches the same residual in roughly one third of the iterations (127 to 42 on the CORE2 operator); because the adjoint of the solve is obtained by implicit differentiation, gradients are unchanged to 10−13 . With both on, annual means of a one-year integration deviate from the production configuration by no more than repeated baseline runs deviate from one another, while the end-to-end production loop accelerates by 21–53 %, the gain growing with mesh size because the host-forcing and solver shares both grow. On NG5 the on-device forcing path

21

Table 4. The four halo-exchange transports, and measured per-step time (full model, uniform benchmark configuration of the transport study, so that the ranking rather than the absolute times is the object of comparison) at the two ends of the mesh range. “Adjoint” is whether reversemode is exact; P is the device count. The all-gather is impractical at NG5 (it gathers the full ∼4 GB field onto every device). The ragged all-to-all is forward-only because its transpose is defective in JAX 10.1. Transport

Adjoint

Per-device volume

CORE2/4

NG5/64

all-gather broadcast

exact

whole field

0.090

—

ragged all-to-all

forward-only

halo only, P -way

0.086

0.589

padded all-to-all

exact

halo × P slots

0.080

0.742

coloured ppermute

exact

halo × K rounds

0.093

0.542

currently costs more device time than it saves on the host, a compiler-optimization pathology under investigation, so only the Chebyshev preconditioner is worth enabling there (−4 to −5 %). 5.9

The cost of a gradient

Reverse-mode differentiation has a measurable cost: modest in time and larger in memory. On one A100, with per-step checkpointing, a reverse-mode step of the complete model, defined here as the value_and_grad of a scalar diagnostic over an N -step window including the checkpointed recomputation of the forward pass, costs 0.92 s, compared with 0.196 s for the forward step. This is a factor of 4.7 and is stable between N = 12 and N = 24 windows. Compiling the reverse program is a one-time cost of about six minutes. Memory, rather than time, is the binding constraint on the window: the backward pass stores one model state per step, measured at 0.83 GiB per step on CORE2 on top of a ∼19 GiB working set (29 GiB at N = 12, 39 GiB at N = 24). Day-scale adjoint windows therefore fit on one 80-GB device, while longer windows use the two-level checkpointing of Sect. 2.2, trading recomputation for storage. 5.10

CPU execution and its scope

The same source also runs on CPUs. This is a convenience rather than a performance claim: it supports development, continuous integration, teaching, and laptop-scale experimentation. On one 128-core node, the full CORE2 model executes at 5.76 s per step (0.9 SYPD) as a single process parallelized by XLA’s thread pool, with no partition files. Launched as sixteen processes of eight cores each under the shard_map layer, with the coloured halo transport, it runs at 1.10 s per step (4.5 SYPD), which also exercises the sharded code path, including gradients, without a GPU. It scales across nodes: DARS takes 16.4 s per step on eight nodes and 8.6 s on sixteen, a speed-up of 1.91×. NG5 does not run on CPU at all, because every process builds the global model state before taking its share and so needs close to the whole-model footprint whatever the process count. The gap to the Fortran original has two parts. First, the individual kernels run about 2.7× slower, the price an array framework pays for generating code from whole-array operations rather than from hand-written loops. Veros reports 1.0–1.4× against

22

its own Fortran ancestor (Häfner et al., 2021). Second, parallel efficiency is about 2.5× worse, and the cause is not the halo exchange but the global sums: the CPU all-reduce grows linearly with process count where MPI’s grows logarithmically, and the surface-elevation solve needs two per iteration, so adding processes stops reducing the time per step beyond about 32. The Fortran original runs 512 ranks and reaches 87.8 SYPD on four CPU nodes, faster than any of our GPU configurations. Many-CPU production therefore belongs to the Fortran original, while the JAX shadow targets GPU execution, gradient-based applications, and the rapid implementation of new model capabilities.

6

Outlook: differentiability

The capability that most clearly distinguishes the shadow from the Fortran model it shadows is differentiability. Because the entire time loop, including ocean dynamics, tracer transport, vertical mixing, and sea-ice rheology, is written as pure JAX functions, the gradient of any scalar diagnostic with respect to any input, boundary condition, or model parameter is available by reverse-mode automatic differentiation through the full integration, at a cost independent of the number of parameters and measured in Sect. 5. Figure 9 illustrates this with two parameter-sensitivity maps on the CORE2 mesh. In each, a scalar parameter is promoted to a field with one value per surface vertex, so that a single reverse-mode pass returns the complete map ∂J/∂θ(x); the same map by finite differences would need one forward integration per vertex, here 1.3×105 runs. The top panel differentiates the globalmean mixed-layer depth with respect to the TKE mixing parameter ck , the bottom the mean temperature of the upper ten layers with respect to the Gent–McWilliams/Redi coefficient kgm . Both are ocean-only, with sea ice off and only the differentiated parameterization active. Differentiation through the sea ice, including the 120-iteration mEVP subcycle, is available, but the iterated explicit rheology amplifies reverse-mode signals over the window, so our gradient applications to date hold the ice fixed; sea ice has required comparable care in the adjoint of the MIT general circulation model (Heimbach et al., 2010). The sum of each map matches a central finite difference of the forward model to relative 6 × 10−7 , a pointwise finite difference reproduces the local value at the most sensitive vertex, and the kgm gradient agrees to within 7 % with an independent twelve-member forward-ensemble estimate. Component-by-component gradient checks run in the continuous-integration suite (Sect. 2.4). Forward-mode differentiation of the same source yields the tangent-linear model, which needs no stored trajectory: its memory does not grow with the window length. The two modes are complementary: the adjoint returns the sensitivity of one scalar diagnostic to a field of parameters at a cost independent of the field’s size, while the tangent-linear model returns the response of every model field, at every point, to one scalar parameter, at a cost independent of the number of outputs. Because the two are transposes of the same linearization, the same scalar derivative can be computed both ways. In half-day windows on the CORE2 mesh, with the full physics active and the ice held fixed in the derivative, the two routes agree to better than 1 %. Two caveats bound the capability. A gradient through chaotic dynamics is informative only over finite horizons, beyond which sensitivities grow exponentially (Lea et al., 2000; Köhl and Willebrand, 2002), so the six-hour maps of Fig. 9 are instan-

23

Figure 9. Adjoint sensitivity maps computed by reverse-mode automatic differentiation through FESOM2-JAX on the CORE2 mesh. Top: gradient of the global-mean mixed-layer depth with respect to a vertex-wise field of the TKE mixing parameter ck ; bottom: gradient of upper-ocean (top ten layers) mean temperature with respect to the Gent–McWilliams/Redi coefficient kgm . Each map is one backward pass over a six-hour window from the PHC initial state under JRA55-do forcing. Circles mark the vertices at which the map value is verified against a pointwise finite difference; the inset compares the scalar kgm gradient, the sum of the map, with the forward-ensemble estimate.

taneous sensitivities rather than an equilibrium response. And the model contains switches: sea-ice concentration thresholds, the triggers in the mixing schemes, and the flux limiter in tracer advection. The gradient is defined everywhere except at those switching points, which is why the finite-difference checks of Sect. 2.4 are taken at states away from them. What the capability enables, gradient-based calibration against observations and hybrid models whose trainable components are trained through the dynamical core, is ongoing work and lies beyond the scope of this paper.

7

Conclusions

We have presented FESOM2-JAX, a re-implementation of the unstructured-mesh ocean–sea-ice model FESOM2 in Python and JAX. The Fortran original is a production code: it is the ocean component of the AWI climate model for CMIP7 and runs in kilometre-scale initiatives including the Destination Earth Climate Digital Twin. FESOM2-JAX is its code shadow, a projection of the same model onto Python and its machine-learning ecosystem, kept demonstrably faithful to the original that casts it, and we have evaluated it against the three requirements such a shadow must meet.

24

It is faithful. In a 1958–2019 hindcast at 1◦ -equivalent resolution, with matched physics and forcing, the mean states of the two codes differ from each other by two orders of magnitude less than either differs from observations, and the runs agree for six decades in temperature, heat content and the sea-ice seasonal cycle. It is usable and it is fast. The model is a NumPy-style codebase that a student can read and modify; a run is specified by a single YAML file; restarts and output do not depend on how many devices produced them; and a regression test fails on any change that alters the shipped configuration’s results. The same source runs, unchanged, from a laptop CPU to 256 GPUs. The complete model at 1◦ fits on a single GPU, and a node of four modern GPUs (GH200 superchips) integrates 112.8 simulated years per wall-clock day, so the hindcast is about half a day of computation; meshes of up to 7.4 million surface vertices, about 5 km, scale to 128 GPUs. What limits the model is communication rather than arithmetic. The small meshes stop gaining beyond one node; the large ones scale only because the halo exchange is implemented so that its cost per device does not grow with the number of devices. What the shadow adds to the original is the gradient, a capability that has been hard for ocean models to acquire and maintain, and one the Fortran original never had. Because the whole time loop is differentiable, a single backward pass returns the sensitivity of a model diagnostic to a parameter at every one of the mesh’s 1.3 × 105 surface vertices at once, where the same map by finite differences would take one integration per vertex. Differentiated forward, the same loop yields the tangentlinear model: the response of every field to one parameter, with no stored trajectory. Both derivatives come from the source that defines the forward model and remain consistent with it as the code evolves. Among models written natively in differentiable frameworks, FESOM2-JAX is to our knowledge the first global ocean–seaice model of CMIP-class complexity, and the first on an unstructured mesh. What this paper demonstrates is the capability itself, on short-window sensitivities verified against finite differences. The applications it points to are potential rather than demonstrated. The same machinery underpins ocean state estimation and variational data assimilation. Model parameters could be calibrated against observations by gradient descent, through the tangent-linear model when the tuned coefficients are few and through the adjoint when they are a field. Sensitivity maps computed with respect to the ocean state rather than to a parameter show where an observation would most constrain a chosen diagnostic, the starting point of observing-system design. Because the derivative extends to whatever is added to the code, neural-network components could be trained through the ocean physics rather than beside it. Realizing them is ongoing work. A production model ordinarily faces a choice between its validated code and the platforms where new methods are developed. Code shadows avoid that choice: the Fortran model remains the single source of truth; its shadows (a C++/Kokkos translation for performance portability, Koldunov et al., 2026b, and FESOM2-JAX for experimentation and gradients) project it onto new technology stacks and are verified by running alongside it; and what a shadow produces (tuned parameter sets, trained machine-learning components, validated numerical variants) returns to the production model. Through its shadows, a Fortran model meets modern hardware, machine learning, and a new generation of developers without a full rewrite.

25

Code availability. FESOM2-JAX is developed openly at https://github.com/koldunovn/fesom_jax; the exact version used here will be tagged and archived at Zenodo on acceptance (DOI: TODO). The reduction and plotting scripts that produce every figure in this paper from the derived data will be archived together with the code on acceptance; they depend only on standard scientific-Python libraries.

Data availability. The model is forced by JRA55-do (Tsujino et al., 2018). The observational references are the PHC3.0 hydrography (Steele et al., 2001), which is also the initial condition, and the OSI-SAF sea-ice concentration climate-data record. A runnable data package for the CORE2 configuration (the mesh, the derived PHC initial state, the domain decompositions, and one year of JRA55-do forcing) is archived at Zenodo (https://doi.org/10.5281/zenodo.21324319) and retrieved by a fetch script shipped with the code. The small derived arrays that produce every figure will be shared on acceptance, archived together with the scripts; the raw model output (the CORE2 and FORCA20 hindcasts and the matched Fortran runs) is several terabytes and is available from the authors on request.

Author contributions. TODO: author contributions to be completed.

Competing interests. The authors declare that they have no conflict of interest.

Acknowledgements. We thank the FESOM2 development team. Computations were performed on the DKRZ Levante supercomputer and, for the GH200 scaling campaign of Sect. 5, on the JUPITER supercomputer at the Jülich Supercomputing Centre.

Financial support. This work was primarily supported through the core funding of the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (AWI), within the Helmholtz Association. Additional support for individual authors was provided through the following projects and funding sources. The contribution by NK, IK, DP and SL was supported by the TerraDT (Digital Twin of Earth system for Cryosphere, Land surface and related interactions) project, which has received funding from the European Union’s Horizon Europe research and innovation programme under Grant Agreement no. 101187992. SD, PS and NK were supported by projects M5, S2 and S1 of the Collaborative Research Centre TRR181 “Energy Transfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 274762653. SC is funded by the WarmWorld Better project of the German Federal Ministry of Research, Technology and Space under the funding code 01LK2202A; the responsibility for the content of this publication lies with the authors. The work of SB has been supported by the European Union’s Destination Earth Initiative and relates to tasks entrusted by the European Union to the European Centre for Medium-Range Weather Forecasts implementing part of this Initiative with funding by the European Union. Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. AK was supported by the HClimRep project, funded by the Helmholtz Foundation Model Initiative (HFMI). TJ was supported by the EERIE project (Grant Agreement No 101081383) funded by the European Union. Views and opinions expressed are however those of the author(s) only and do

26

not necessarily reflect those of the European Union or the European Climate Infrastructure and Environment Executive Agency (CINEA). Neither the European Union nor the granting authority can be held responsible for them.

Appendix A: Component-by-component description of the ported model This appendix collects the description of the ported components, so that the model can be read without consulting Danilov et al. (2017) and Scholz et al. (2019, 2022). Numerical settings quoted here are those used in the present study and are summarized in Table 1. A1

Continuous equations and horizontal discretization

The ocean component solves the hydrostatic, Boussinesq primitive equations. The horizontal mesh is an unstructured triangulation, and the discretization is cell-vertex finite-volume. Scalar fields, including temperature, salinity, pressure, and sea-surface height, are located at mesh vertices and advanced on the control volumes, while horizontal velocity is located at triangle centroids (Danilov et al., 2017). Horizontal operators are therefore indirect: divergences, gradients, and fluxes are assembled by gathering values along precomputed vertex, edge, and element neighbor lists and scattering edge-wise antisymmetric contributions back. The same gather/scatter stencil structure, rather than structuredgrid array shifts, gives the model its smoothly variable resolution and also accounts for much of its computational cost (Sect. 5). Lateral boundaries are no-slip. A2

Vertical coordinate

The vertical coordinate is the arbitrary Lagrangian–Eulerian (ALE) coordinate of FESOM2, used throughout this study in its z ⋆ form (Adcroft and Campin, 2004). In this formulation, the free-surface displacement is distributed uniformly over the stretchable part of the water column, so layer thicknesses vary with the change in sea-surface height while layer interfaces retain their nominal stacking. Two consequences are important for what follows. First, surface freshwater exchange is represented as a real volume flux, not as a virtual salt flux, so the global freshwater budget must close exactly. The salinity-affecting bug in the python JAX port, discussed in Sect. 4, violated precisely this closure. Second, the pressure-gradient force must be evaluated on the time-varying geometry. FESOM2-JAX therefore uses the density-Jacobian formulation of Shchepetkin and McWilliams (2003), as the original does for z ⋆ . The linear free surface option is also available. A3

Equation of state

In-situ density and hydrostatic pressure are computed from the polynomial equation of state of Jackett and McDougall (1995); the thermal-expansion and haline-contraction coefficients that feed the mixing and eddy parameterizations follow McDougall (1987), and the Brunt–Väisälä frequency is smoothed once over the vertex neighbourhood, as in the original.

27

A4

Time stepping and the free surface

One model timestep advances momentum, the free surface, and the tracers in sequence, with the sea-ice component executed first (Sect. A6). Momentum uses a second-order Adams–Bashforth estimate, with FESOM2’s stabilizing offset ϵ = 0.1, for the Coriolis and advection terms; momentum advection is evaluated in flux form on the scalar control volumes. Sea-surface height is advanced by treating the external mode semi-implicitly, with two implicitness weights in the notation of Danilov et al. (2017): α on the transport divergence in the elevation equation, and θ on the surface elevation in the momentum equation. Setting both to 1/2 recovers the Crank–Nicolson scheme; we use FESOM2’s default α = θ = 1, which is fully implicit. The resulting elliptic system is re-assembled on each time step under the full free surface and is static under the linear free surface. It is solved at every step by a preconditioned conjugate-gradient (CG) iteration, with the symmetric preconditioner of Marshall et al. (1997) that FESOM2 uses, to a relative residual of 10−5 . The solve is warm-started from the previous step solution and typically converges in a few iterations. Layer thicknesses are then updated from the new surface height (z ⋆ ), and the quasi-vertical diasurface velocity is diagnosed by integrating the horizontal transport divergence upward from the no-flux bottom. Vertical viscosity and diffusion are treated implicitly with tridiagonal solves in each column, removing the vertical CFL restriction from mixing. FESOM2’s w_split, which moves the part of the vertical advective velocity that exceeds the explicit CFL limit into the implicit solve, is also ported. It is normally enabled on high-resolution meshes to improve stability, but none of the runs reported here use it. The CORE2 configuration of Sect. 4 runs at ∆t = 1800, s. A5

Interior physics and surface fluxes

Tracers are advected with the third-order upwind MUSCL-type scheme (Danilov et al., 2017) in the horizontal direction and the fourth-order centered scheme in the vertical direction, limited by the flux-corrected-transport (FCT) scheme of Zalesak (1979): a first-order upwind solution is corrected by limited antidiffusive fluxes from the high-order estimate. The advected field enters the fluxes through the same Adams–Bashforth interpolation as the Coriolis and momentum advection terms. Vertical mixing uses the prognostic turbulent-kinetic-energy closure of Gaspar et al. (1990) in the CVMix formulation adopted by FESOM2 (cvmix-TKE; Griffies et al., 2015; Van Roekel et al., 2018; Scholz et al., 2022): one TKE equation per column, integrated implicitly, with the standard constants (ck = 0.1, cϵ = 0.7) and background TKE value and viscosity of 10−5 and 10−4 m2 s−1 respectively. The K-profile parameterization (Large et al., 1994) and the Richardson-number scheme of Pacanowski and Philander (1981) are also ported and configuration-selectable; all runs in this paper use cvmix-TKE. Where the mesh does not resolve eddies, eddy effects are parameterized by Gent–McWilliams bolus advection (Gent and McWilliams, 1990) and Redi isoneutral diffusion (Redi, 1982). As in FESOM2, the bolus transport is obtained by solving a vertical boundary-value problem for the eddy streamfunction in each column (Ferrari et al., 2010) rather than from local slopes alone; the transfer coefficient (here capped at 1000 m2 s−1 ) is scaled with the local horizontal resolution and attenuated in the vertical, and neutral slopes are tapered where they steepen. The parameterization is a pure diagnostic of the current density

28

field, recomputed every step. It is active in the CORE2 hindcast and switched off in the configurations used for the performance measurements of Sect. 5. Long integrations on fArc, which is coarse outside its refined region, also use it. Momentum dissipation uses FESOM2’s biharmonic filter, assembled edge-wise in two passes of an approximate Laplacian with a coefficient set by the local mesh scale and the velocity difference across the edge, so that dissipation concentrates where the flow is rough and relaxes to a small background value where it is smooth. The operator is that of Danilov et al. (2017), iterated to biharmonic form; the flow-aware coefficient is the one introduced by Juricke et al. (2020), who note that, being based on velocity differences, it is reminiscent of the Smagorinsky viscosity (Smagorinsky, 1963). The model is driven by the JRA55-do reanalysis dataset (Tsujino et al., 2018): of three-hourly fields (10-m winds, air temperature, specific humidity, downward short- and longwave radiation, rain, and snow). They are interpolated bilinearly to the mesh and linearly in time, and turbulent fluxes are computed interactively with the NCAR bulk formulae (Large and Yeager, 2009), whose Monin–Obukhov stability iteration runs a fixed number of passes (a choice made for differentiability, Sect. 2.2). Shortwave radiation penetrates the water column with a chlorophyll-dependent attenuation profile (Sweeney et al., 2005). River runoff is a climatological field, and sea-surface salinity is restored weakly to a monthly climatology with a piston velocity of 10 m per 60 days. A6

Sea ice

Sea ice is FESOM2’s single-class dynamic–thermodynamic model, based on the Finite-Element Sea Ice Model (FESIM; Danilov et al., 2015) and discretized at the same mesh vertices as the ocean scalars. Its prognostic variables are ice concentration, ice and snow volume per unit area, and ice velocity at the vertices. The momentum balance uses the viscous–plastic rheology with an elliptical yield curve (aspect ratio 2) and the Hibler (1979) strength P = P ∗ h e−C(1−A) (here P ∗ = 3 × 104 N m−2 , C = 20), solved with the modified elastic–viscous–plastic (mEVP) iteration (Bouillon et al., 2013; Kimmritz et al., 2015), the scheme FESOM2 uses in production for its favourable convergence and cost (Koldunov et al., 2019b). Each ocean timestep runs a fixed 120 pseudotime iterations with stabilization parameters α = β = 250; each iteration evaluates element strain rates and stresses, scatters the stress divergence to the vertices, and solves the pointwise implicit momentum update including Coriolis and an ice–ocean drag (coefficient 5.5 × 10−3 ). Standard EVP is retained as a configuration option. The three ice tracers are advected with a Taylor–Galerkin finite-element FCT scheme of Löhner et al. (1987), the sea-ice finite-element counterpart of the ocean FCT. Thermodynamics is the zero-layer scheme of Semtner (1976): a fixed-iteration Newton solve for the ice skin temperature, growth rates averaged over seven ice-thickness classes, lead closing following Hibler (1979), snow-to-ice conversion by flooding, and a salinity-dependent freezing point. The ice component runs first within each timestep: it receives the ocean surface state (temperature, salinity, surface currents, and the sea-surface tilt), executes dynamics, advection, and thermodynamics, and returns to the ocean the ice-modulated surface stress, the net heat flux, the shortwave flux gated by ice cover, and the surface salinity forcing. Under z ⋆ the ice–ocean freshwater exchange is a real volume flux with an explicit salt flux for the brine content, and the global-mean freshwater flux is

29

removed each step so that the free surface conserves volume — the budget closure whose violation produced the salinity drift dissected in Sect. 4. A7

Nominal and area-weighted mesh resolution

The nominal resolution of Table 2 is a design figure, not a uniform one: a mesh built to concentrate resolution is coarse over most of its area. On fArc, which is CORE2 outside the Arctic, four-fifths of the vertices lie in the 4.5 km refined region, yet the area-weighted mean resolution over the global ocean is 85 km. Koldunov et al. (2026a) compare the latitude-based and eddyrefined mesh families directly and find that eddy-refined meshes reproduce global energetics comparable to latitude-based ones at lower cost.

30

References Adams, S. V., Ford, R. W., Hambley, M., Hobson, J. M., Kavčič, I., Maynard, C. M., et al.: LFRic: Meeting the challenges of scalability and performance portability in weather and climate models, Journal of Parallel and Distributed Computing, 132, 383–396, https://doi.org/10.1016/j.jpdc.2019.02.007, 2019. Adcroft, A. and Campin, J.-M.: Rescaled height coordinates for accurate representation of free-surface flows in ocean circulation models, Ocean Modelling, 7, 269–284, https://doi.org/10.1016/j.ocemod.2003.09.003, 2004. Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M.: Automatic differentiation in machine learning: a survey, Journal of Machine Learning Research, 18, 1–43, 2018. Bezanson, J., Edelman, A., Karpinski, S., and Shah, V. B.: Julia: A fresh approach to numerical computing, SIAM Review, 59, 65–98, https://doi.org/10.1137/141000671, 2017. Bezgin, D. A., Buhendwa, A. B., and Adams, N. A.: JAX-Fluids: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows, Computer Physics Communications, 282, 108 527, https://doi.org/10.1016/j.cpc.2022.108527, 2023. Bouillon, S., Fichefet, T., Legat, V., and Madec, G.: The elastic–viscous–plastic method revisited, Ocean Modelling, 71, 2–12, https://doi.org/10.1016/j.ocemod.2013.05.013, 2013. Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., and Zhang, Q.: JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax, 2018. Dahm, J., Davis, E., Deconinck, F., Elbert, O., George, R., McGibbon, J., others, and Fuhrer, O.: Pace v0.2: a Python-based performanceportable atmospheric model, Geoscientific Model Development, 16, 2719–2736, https://doi.org/10.5194/gmd-16-2719-2023, 2023. Danilov, S., Wang, Q., Timmermann, R., Iakovlev, N., Sidorenko, D., Kimmritz, M., Jung, T., and Schröter, J.: Finite-Element Sea Ice Model (FESIM), version 2, Geoscientific Model Development, 8, 1747–1761, https://doi.org/10.5194/gmd-8-1747-2015, 2015. Danilov, S., Sidorenko, D., Wang, Q., and Jung, T.: The Finite-volumE Sea ice–Ocean Model (FESOM2), Geoscientific Model Development, 10, 765–789, https://doi.org/10.5194/gmd-10-765-2017, 2017. Doblas-Reyes, F. J., Kontkanen, J., Sandu, I., Acosta, M., Al Turjmam, M. H., Alsina-Ferrer, I., Andrés-Martínez, M., Anerdi, C., Arriola, L., Axness, M., Batlle Martín, M., Bauer, P., Becker, T., Beltrán, D., Beyer, S., Bockelmann, H., Bretonnière, P.-A., Cabaniols, S., Caprioli, S., Castrillo, M., Chandrasekar, A., Cheedela, S., Correal, V., Danovaro, E., Davini, P., Enkovaara, J., Frauen, C., Früh, B., Gaya Àvila, A., Ghinassi, P., Ghosh, R., Ghosh, S., González, I., Grayson, K., Griffith, M., Hadade, I., Haine, C., Hartick, C., Haus, U.-U., Hearne, S., Järvinen, H., Jiménez, B., John, A., Juchem, M., Jung, T., Kegel, J., Kelbling, M., Keller, K., Kinoshita, B., Kiszler, T., Klocke, D., Kluft, L., Koldunov, N., Kölling, T., Kolstela, J., Kornblueh, L., Kosukhin, S., Lacima-Nadolnik, A., Leal Rojas, J. J., Lehtiranta, J., Lunttila, T., Luoma, A., Manninen, P., Medvedev, A., Milinski, S., Mohammed, A., Müller, S., Naryanappa, D., Nazarova, N., Niemelä, S., Niraula, B., Nortamo, H., Nummelin, A., Nurisso, M., Ortega, P., Paronuzzi, S., Pedruzo-Bagazgoitia, X., Pelletier, C., Peña, C., Polade, S., Pradhan, H. K., Quintanilla, R., Quintino, T., Rackow, T., Räisänen, J., Rajput, M. M., Redler, R., Reuter, B., Rocha Monteiro, N., Roura-Adserias, F., Ruppert, S., Sayed, S., Schnur, R., Sharma, T., Sidorenko, D., Sievi-Korte, O., Soret, A., Steger, C., Stevens, B., Streffing, J., Sunny, J., Tenorio, L., Thober, S., Tigerstedt, U., Tinto, O., Tonttila, J., Tuomenvirta, H., Tuppi, L., Van Thielen, G., Vitali, E., von Hardenberg, J., Wagner, I., Wedi, N., Wehner, J., Willner, S., Yepes-Arbós, X., Ziemen, F., and Zimmermann, J.: The Destination Earth digital twin for climate change adaptation, Geoscientific Model Development, 19, 2821–2848, https://doi.org/10.5194/gmd-19-2821-2026, 2026. Errico, R. M.: What is an adjoint model?, Bulletin of the American Meteorological Society, 78, 2577–2591, https://doi.org/10.1175/15200477(1997)078<2577:WIAAM>2.0.CO;2, 1997.

31

EUMETSAT Ocean and Sea Ice Satellite Application Facility: Global sea ice concentration climate data record 1978–2020 (OSI-450-a), version 3.0, https://doi.org/10.15770/EUM_SAF_OSI_0013, 2022. Ferrari, R., Griffies, S. M., Nurser, A. J. G., and Vallis, G. K.: A boundary-value problem for the parameterized mesoscale eddy transport, Ocean Modelling, 32, 143–156, https://doi.org/10.1016/j.ocemod.2010.01.004, 2010. Forget, G., Campin, J.-M., Heimbach, P., Hill, C. N., Ponte, R. M., and Wunsch, C.: ECCO version 4: an integrated framework for non-linear inverse modeling and global ocean state estimation, Geoscientific Model Development, 8, 3071–3104, https://doi.org/10.5194/gmd-83071-2015, 2015. Frezat, H., Le Sommer, J., Fablet, R., Balarac, G., and Lguensat, R.: A posteriori learning for quasi-geostrophic turbulence parametrization, Journal of Advances in Modeling Earth Systems, 14, e2022MS003 124, https://doi.org/10.1029/2022MS003124, 2022. Gärtner, J. P., Losch, M., Cheedela, S. K., Jochum, M., and Nuterman, R.: Veris: fast & efficient sea-ice modeling in Python with GPU acceleration, Geoscientific Model Development, 19, 5225–5236, https://doi.org/10.5194/gmd-19-5225-2026, 2026. Gaspar, P., Grégoris, Y., and Lefevre, J.-M.: A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and long-term upper ocean study site, Journal of Geophysical Research: Oceans, 95, 16 179–16 193, https://doi.org/10.1029/JC095iC09p16179, 1990. Gelbrecht, M., White, A., Bathiany, S., and Boers, N.: Differentiable programming for Earth system modeling, Geoscientific Model Development, 16, 3123–3135, https://doi.org/10.5194/gmd-16-3123-2023, 2023. Gent, P. R. and McWilliams, J. C.: Isopycnal mixing in ocean circulation models, Journal of Physical Oceanography, 20, 150–155, https://doi.org/10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2, 1990. Giering, R. and Kaminski, T.: Recipes for adjoint code construction, ACM Transactions on Mathematical Software, 24, 437–474, https://doi.org/10.1145/293686.293695, 1998. Giorgetta, M. A., Sawyer, W., Lapillonne, X., Adamidis, P., Alexeev, D., Clément, V., et al.: The ICON-A model for direct QBO simulations on GPUs (version icon-cscs:baf28a514), Geoscientific Model Development, 15, 6985–7016, https://doi.org/10.5194/gmd-15-6985-2022, 2022. Griewank, A. and Walther, A.: Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation, ACM Transactions on Mathematical Software, 26, 19–45, https://doi.org/10.1145/347837.347846, 2000. Griffies, S. M., Levy, M., Adcroft, A. J., Danabasoglu, G., Hallberg, R. W., Jacobsen, D., Large, W. G., and Ringler, T. D.: Theory and numerics of the Community Ocean Vertical Mixing (CVMix) project, Tech. rep., NOAA/Geophysical Fluid Dynamics Laboratory, Princeton, NJ, available at: https://github.com/CVMix/CVMix-description (last access: 2 August 2026), 2015. Häfner, D. and Vicentini, F.: mpi4jax: Zero-copy MPI communication of JAX arrays, Journal of Open Source Software, 6, 3419, https://doi.org/10.21105/joss.03419, 2021. Häfner, D., Jacobsen, R. L., Eden, C., Kristensen, M. R. B., Jochum, M., Nuterman, R., and Vinter, B.: Veros v0.1 – a fast and versatile ocean simulator in pure Python, Geoscientific Model Development, 11, 3299–3312, https://doi.org/10.5194/gmd-11-3299-2018, 2018. Häfner, D., Nuterman, R., and Jochum, M.: Fast, cheap, and turbulent — Global ocean modeling with GPU acceleration in Python, Journal of Advances in Modeling Earth Systems, 13, e2021MS002 717, https://doi.org/10.1029/2021MS002717, 2021. Heimbach, P., Hill, C., and Giering, R.: An efficient exact adjoint of the parallel MIT general circulation model, generated via automatic differentiation, Future Generation Computer Systems, 21, 1356–1371, https://doi.org/10.1016/j.future.2004.11.010, 2005.

32

Heimbach, P., Menemenlis, D., Losch, M., Campin, J.-M., and Hill, C.: On the formulation of sea-ice models. Part 2: Lessons from multi-year adjoint sea-ice export sensitivities through the Canadian Arctic Archipelago, Ocean Modelling, 33, 145–158, https://doi.org/10.1016/j.ocemod.2010.02.002, 2010. Hibler, W. D.: A dynamic thermodynamic sea ice model, Journal of Physical Oceanography, 9, 815–846, https://doi.org/10.1175/15200485(1979)009<0815:ADTSIM>2.0.CO;2, 1979. Iglesias, M. A., Law, K. J. H., and Stuart, A. M.: Ensemble Kalman methods for inverse problems, Inverse Problems, 29, 045 001, https://doi.org/10.1088/0266-5611/29/4/045001, 2013. Jackett, D. R. and McDougall, T. J.: Minimal adjustment of hydrographic profiles to achieve static stability, Journal of Atmospheric and Oceanic Technology, 12, 381–389, https://doi.org/10.1175/1520-0426(1995)012<0381:MAOHPT>2.0.CO;2, 1995. Juricke, S., Danilov, S., Koldunov, N., Oliver, M., Sein, D. V., Sidorenko, D., and Wang, Q.: A kinematic kinetic energy backscatter parametrization: From implementation to global ocean simulations, Journal of Advances in Modeling Earth Systems, 12, e2020MS002 175, https://doi.org/10.1029/2020MS002175, 2020. Kimmritz, M., Danilov, S., and Losch, M.: On the convergence of the modified elastic–viscous–plastic method for solving the sea ice momentum equation, Journal of Computational Physics, 296, 90–100, https://doi.org/10.1016/j.jcp.2015.04.051, 2015. Kochkov, D., Smith, J. A., Alieva, A., Wang, Q., Brenner, M. P., and Hoyer, S.: Machine learning-accelerated computational fluid dynamics, Proceedings of the National Academy of Sciences, 118, e2101784 118, https://doi.org/10.1073/pnas.2101784118, 2021. Kochkov, D., Yuval, J., Langmore, I., Norgaard, P., Smith, J., Mooers, G., Klöwer, M., others, and Hoyer, S.: Neural general circulation models for weather and climate, Nature, 632, 1060–1066, https://doi.org/10.1038/s41586-024-07744-y, 2024. Köhl, A. and Willebrand, J.: An adjoint method for the assimilation of statistical characteristics into eddy-resolving ocean models, Tellus A, 54, 406–425, https://doi.org/10.3402/tellusa.v54i4.12150, 2002. Koldunov, A., Sidorenko, D., Koldunov, N., Danilov, S., and Scholz, P.: Lorenz Energy Cycle of the Global Eddying Ocean Simulated on Meshes With Different Designs and Resolutions, Journal of Advances in Modeling Earth Systems, 18, e2025MS005 551, https://doi.org/10.1029/2025MS005551, 2026a. Koldunov, N. V., Aizinger, V., Rakowsky, N., Scholz, P., Sidorenko, D., Danilov, S., and Jung, T.: Scalability and some optimization of the Finite-volumE Sea ice–Ocean Model, Version 2.0 (FESOM2), Geoscientific Model Development, 12, 3991–4012, https://doi.org/10.5194/gmd-12-3991-2019, 2019a. Koldunov, N. V., Danilov, S., Sidorenko, D., Hutter, N., Losch, M., Goessling, H., Rakowsky, N., Scholz, P., Sein, D., Wang, Q., and Jung, T.: Fast EVP solutions in a high-resolution sea ice model, Journal of Advances in Modeling Earth Systems, 11, 1269–1284, https://doi.org/10.1029/2018MS001485, 2019b. Koldunov, N. V., Cheedela, S. K., Danilov, S., Sidorenko, D., Beyer, S., and Jung, T.: An ocean model ported by a large language model: experience and lessons from FESOM2 (Fortran to C to C++/Kokkos), arXiv preprint arXiv:2606.11356, https://doi.org/10.48550/arXiv.2606.11356, 2026b. Lahlou, A., Hawkins, L., and Gentine, P.: Systematic LLM translation of legacy scientific code to differentiable frameworks: application to a land surface model, arXiv preprint arXiv:2606.07681, https://doi.org/10.48550/arXiv.2606.07681, 2026. Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying air–sea flux data set, Climate Dynamics, 33, 341–364, https://doi.org/10.1007/s00382-008-0441-3, 2009. Large, W. G., McWilliams, J. C., and Doney, S. C.: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization, Reviews of Geophysics, 32, 363–403, https://doi.org/10.1029/94RG01872, 1994.

33

Lea, D. J., Allen, M. R., and Haine, T. W. N.: Sensitivity analysis of the climate of a chaotic system, Tellus A, 52, 523–532, https://doi.org/10.3402/tellusa.v52i5.12283, 2000. Löhner, R., Morgan, K., Peraire, J., and Vahdati, M.: Finite element flux-corrected transport (FEM–FCT) for the Euler and Navier–Stokes equations, International Journal for Numerical Methods in Fluids, 7, 1093–1109, https://doi.org/10.1002/fld.1650071007, 1987. Marshall, J., Adcroft, A., Hill, C., Perelman, L., and Heisey, C.: A finite-volume, incompressible Navier–Stokes model for studies of the ocean on parallel computers, Journal of Geophysical Research: Oceans, 102, 5753–5766, https://doi.org/10.1029/96JC02775, 1997. McDougall,

T.

J.:

Neutral

surfaces,

Journal

of

Physical

Oceanography,

17,

1950–1964,

https://doi.org/10.1175/1520-

0485(1987)017<1950:NS>2.0.CO;2, 1987. Moon, J.-Y., Streffing, J., Lee, S.-S., Semmler, T., Andrés-Martínez, M., Chen, J., Cho, E.-B., Chu, J.-E., Franzke, C. L., Gärtner, J. P., et al.: Earth’s future climate and its variability simulated at 9 km global resolution, Earth System Dynamics, 16, 1103–1134, https://doi.org/10.5194/esd-16-1103-2025, 2025. Moore, A. M., Arango, H. G., Broquet, G., Powell, B. S., Weaver, A. T., and Zavala-Garay, J.: The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems: Part I – System overview and formulation, Progress in Oceanography, 91, 34–49, https://doi.org/10.1016/j.pocean.2011.05.004, 2011. Moses, W. S. and Churavy, V.: Instead of rewriting foreign code for machine learning, automatically synthesize fast gradients, in: Advances in Neural Information Processing Systems 33, pp. 12 472–12 485, 2020. Moses, W. S., Cheng, G., Churavy, V., Gelbrecht, M., Klöwer, M., Kump, J., Morlighem, M., Williamson, S., Apte, D., Berg, P., Giordano, M., Hill, C., Loose, N., Montoison, A., Narayanan, S. H. K., Pal, A., Schanen, M., Silvestri, S., Wagner, G., and Heimbach, P.: DJ4Earth: Differentiable, and performance-portable Earth system modeling via program transformations, Journal of Advances in Modeling Earth Systems, 18, e2025MS005 615, https://doi.org/10.1029/2025MS005615, 2026. Pacanowski, R. C. and Philander, S. G. H.: Parameterization of vertical mixing in numerical models of tropical oceans, Journal of Physical Oceanography, 11, 1443–1451, https://doi.org/10.1175/1520-0485(1981)011<1443:POVMIN>2.0.CO;2, 1981. Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., et al.: PyTorch: An imperative style, high-performance deep learning library, in: Advances in Neural Information Processing Systems 32, pp. 8024–8035, 2019. Rackow, T., Pedruzo-Bagazgoitia, X., Becker, T., Milinski, S., Sandu, I., Aguridan, R., Bechtold, P., Beyer, S., Bidlot, J., Boussetta, S., Deconinck, W., Diamantakis, M., Dueben, P., Dutra, E., Forbes, R., Ghosh, R., Goessling, H. F., Hadade, I., Hegewald, J., Jung, T., Keeley, S., Kluft, L., Koldunov, N., Koldunov, A., Kölling, T., Kousal, J., Kühnlein, C., Maciel, P., Mogensen, K., Quintino, T., Polichtchouk, I., Reuter, B., Sármány, D., Scholz, P., Sidorenko, D., Streffing, J., Sützl, B., Takasuka, D., Tietsche, S., Valentini, M., Vannière, B., Wedi, N., Zampieri, L., and Ziemen, F.: Multi-year simulations at kilometre scale with the Integrated Forecasting System coupled to FESOM2.5 and NEMOv3.4, Geoscientific Model Development, 18, 33–69, https://doi.org/10.5194/gmd-18-33-2025, 2025. Ramadhan, A., Wagner, G. L., Hill, C., Campin, J.-M., Churavy, V., Besard, T., Souza, A., Edelman, A., Ferrari, R., and Marshall, J.: Oceananigans.jl: Fast and friendly geophysical fluid dynamics on GPUs, Journal of Open Source Software, 5, 2018, https://doi.org/10.21105/joss.02018, 2020. Rasp, S., Pritchard, M. S., and Gentine, P.: Deep learning to represent subgrid processes in climate models, Proceedings of the National Academy of Sciences, 115, 9684–9689, https://doi.org/10.1073/pnas.1810286115, 2018. Redi, M. H.: Oceanic isopycnal mixing by coordinate rotation, Journal of Physical Oceanography, 12, 1154–1158, https://doi.org/10.1175/1520-0485(1982)012<1154:OIMBCR>2.0.CO;2, 1982.

34

Scholz, P., Sidorenko, D., Gurses, O., Danilov, S., Koldunov, N., Wang, Q., Sein, D., Smolentseva, M., Rakowsky, N., and Jung, T.: Assessment of the Finite-volumE Sea ice–Ocean Model (FESOM2.0) – Part 1: Description of selected key model elements and comparison to its predecessor version, Geoscientific Model Development, 12, 4875–4899, https://doi.org/10.5194/gmd-12-4875-2019, 2019. Scholz, P., Sidorenko, D., Danilov, S., Wang, Q., Koldunov, N., Sein, D., and Jung, T.: Assessment of the Finite-volumE Sea ice–Ocean Model (FESOM2.0) – Part 2: Partial bottom cells, embedded sea ice and vertical mixing library CVMix, Geoscientific Model Development, 15, 335–363, https://doi.org/10.5194/gmd-15-335-2022, 2022. Segura, H., Pedruzo-Bagazgoitia, X., Weiss, P., Müller, S. K., Rackow, T., Lee, J., Dolores-Tesillos, E., Benedict, I., Aengenheyster, M., Aguridan, R., Arduini, G., Baker, A. J., Bao, J., Bastin, S., Baulenas, E., Becker, T., Beyer, S., Bockelmann, H., Brüggemann, N., Brunner, L., Cheedela, S. K., Das, S., Denissen, J., Dragaud, I., Dziekan, P., Ekblom, M., Engels, J. F., Esch, M., Forbes, R., Frauen, C., Freischem, L., García-Maroto, D., Geier, P., Gierz, P., González-Cervera, A., Grayson, K., Griffith, M., Gutjahr, O., Haak, H., Hadade, I., Haslehner, K., ul Hasson, S., Hegewald, J., Kluft, L., Koldunov, A., Koldunov, N., Kölling, T., Koseki, S., Kosukhin, S., Kousal, J., Kuma, P., Kumar, A. U., Li, R., Maury, N., Meindl, M., Milinski, S., Mogensen, K., Niraula, B., Nowak, J., Praturi, D. S., Proske, U., Putrasahan, D., Redler, R., Santuy, D., Sármány, D., Schnur, R., Scholz, P., Sidorenko, D., Spät, D., Sützl, B., Takasuka, D., Tompkins, A., Uribe, A., Valentini, M., Veerman, M., Voigt, A., Warnau, S., Wachsmann, F., Wacławczyk, M., Wedi, N., Wieners, K.-H., Wille, J., Winkler, M., Wu, Y., Ziemen, F., Zimmermann, J., Bender, F. A.-M., Bojovic, D., Bony, S., Bordoni, S., Brehmer, P., Dengler, M., Dutra, E., Faye, S., Fischer, E., van Heerwaarden, C., Hohenegger, C., Järvinen, H., Jochum, M., Jung, T., Jungclaus, J. H., Keenlyside, N. S., Klocke, D., Konow, H., Klose, M., Malinowski, S., Martius, O., Mauritsen, T., Mellado, J. P., Mieslinger, T., Mohino, E., Pawłowska, H., Peters-von Gehlen, K., Sarré, A., Sobhani, P., Stier, P., Tuppi, L., Vidale, P. L., Sandu, I., and Stevens, B.: nextGEMS: entering the era of kilometer-scale Earth system modeling, Geoscientific Model Development, 18, 7735–7761, https://doi.org/10.5194/gmd-18-7735-2025, 2025. Sein, D. V., Koldunov, N. V., Danilov, S., Wang, Q., Sidorenko, D., Fast, I., Rackow, T., Cabos, W., and Jung, T.: Ocean Modeling on a Mesh With Resolution Following the Local Rossby Radius, Journal of Advances in Modeling Earth Systems, 9, 2601–2614, https://doi.org/10.1002/2017MS001099, 2017. Semtner, A. J.: A model for the thermodynamic growth of sea ice in numerical investigations of climate, Journal of Physical Oceanography, 6, 379–389, https://doi.org/10.1175/1520-0485(1976)006<0379:AMFTTG>2.0.CO;2, 1976. Shchepetkin, A. F. and McWilliams, J. C.: A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate, Journal of Geophysical Research: Oceans, 108, 3090, https://doi.org/10.1029/2001JC001047, 2003. Shen, C., Appling, A. P., Gentine, P., Bandai, T., Gupta, H., Tartakovsky, A., et al.: Differentiable modelling to unify machine learning and physical models for geosciences, Nature Reviews Earth & Environment, 4, 552–567, https://doi.org/10.1038/s43017-023-00450-9, 2023. Silvestri, S., Wagner, G. L., Constantinou, N. C., Hill, C. N., Campin, J.-M., Souza, A. N., Bishnu, S., Churavy, V., Marshall, J., and Ferrari, R.: A GPU-based ocean dynamical core for routine mesoscale-resolving climate simulations, Journal of Advances in Modeling Earth Systems, 17, e2024MS004 465, https://doi.org/10.1029/2024MS004465, 2025. Smagorinsky, J.: General circulation experiments with the primitive equations: I. The basic experiment, Monthly Weather Review, 91, 99– 164, https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2, 1963. Souza, A. N., Wagner, G. L., Ramadhan, A., Allen, B., Churavy, V., Schloss, J., Campin, J.-M., Hill, C., Edelman, A., Marshall, J., Flierl, G., and Ferrari, R.: Uncertainty quantification of ocean parameterizations: Application to the K-Profile-Parameterization for penetrative convection, Journal of Advances in Modeling Earth Systems, 12, e2020MS002 108, https://doi.org/10.1029/2020MS002108, 2020.

35

Stammer, D., Wunsch, C., Giering, R., Eckert, C., Heimbach, P., Marotzke, J., Adcroft, A., Hill, C. N., and Marshall, J.: Global ocean circulation during 1992–1997, estimated from ocean observations and a general circulation model, Journal of Geophysical Research: Oceans, 107, 3118, https://doi.org/10.1029/2001JC000888, 2002. Steele, M., Morley, R., and Ermold, W.: PHC: A global ocean hydrography with a high-quality Arctic Ocean, Journal of Climate, 14, 2079– 2087, https://doi.org/10.1175/1520-0442(2001)014<2079:PAGOHW>2.0.CO;2, 2001. Streffing, J., Sidorenko, D., Semmler, T., Zampieri, L., Scholz, P., Andrés-Martínez, M., Koldunov, N., Rackow, T., Kjellsson, J., Goessling, H., Athanase, M., Wang, Q., Hegewald, J., Sein, D. V., Mu, L., Fladrich, U., Barbi, D., Gierz, P., Danilov, S., Juricke, S., Lohmann, G., and Jung, T.: AWI-CM3 coupled climate model: description and evaluation experiments for a prototype post-CMIP6 model, Geoscientific Model Development, 15, 6399–6427, https://doi.org/10.5194/gmd-15-6399-2022, 2022. Sumata, H., Kauker, F., Karcher, M., and Gerdes, R.: Simultaneous parameter optimization of an Arctic sea ice–ocean model by a genetic algorithm, Monthly Weather Review, 147, 1899–1926, https://doi.org/10.1175/MWR-D-18-0360.1, 2019. Sweeney, C., Gnanadesikan, A., Griffies, S. M., Harrison, M. J., Rosati, A. J., and Samuels, B. L.: Impacts of shortwave penetration depth on large-scale ocean circulation and heat transport, Journal of Physical Oceanography, 35, 1103–1119, https://doi.org/10.1175/JPO2740.1, 2005. Thacker, W. C. and Long, R. B.: Fitting dynamics to data, Journal of Geophysical Research: Oceans, 93, 1227–1240, https://doi.org/10.1029/JC093iC02p01227, 1988. Tsujino, H., Urakawa, S., Nakano, H., et al.: JRA-55 based surface dataset for driving ocean–sea-ice models (JRA55-do), Ocean Modelling, 130, 79–139, https://doi.org/10.1016/j.ocemod.2018.07.002, 2018. Utke, J., Naumann, U., Fagan, M., Tallent, N., Strout, M., Heimbach, P., Hill, C., and Wunsch, C.: OpenAD/F: A modular open-source tool for automatic differentiation of Fortran codes, ACM Transactions on Mathematical Software, 34, 1–36, https://doi.org/10.1145/1377596.1377598, article 18, 2008. Van Roekel, L., Adcroft, A. J., Danabasoglu, G., Griffies, S. M., Kauffman, B., Large, W., Levy, M., Reichl, B. G., Ringler, T., and Schmidt, M.: The KPP boundary layer scheme for the ocean: Revisiting its formulation and benchmarking one-dimensional simulations relative to LES, Journal of Advances in Modeling Earth Systems, 10, 2647–2685, https://doi.org/10.1029/2018MS001336, 2018. Vidard, A., Bouttier, P.-A., and Vigilant, F.: NEMOTAM: tangent and adjoint models for the ocean modelling platform NEMO, Geoscientific Model Development, 8, 1245–1257, https://doi.org/10.5194/gmd-8-1245-2015, 2015. Wang, Q., Danilov, S., Sidorenko, D., Timmermann, R., Wekerle, C., Wang, X., Jung, T., and Schröter, J.: The Finite Element Sea IceOcean Model (FESOM) v.1.4: formulation of an ocean general circulation model, Geoscientific Model Development, 7, 663–693, https://doi.org/10.5194/gmd-7-663-2014, 2014. Wang, Q., Wekerle, C., Danilov, S., Wang, X., and Jung, T.: A 4.5 km resolution Arctic Ocean simulation with the global multi-resolution model FESOM 1.4, Geoscientific Model Development, 11, 1229–1255, https://doi.org/10.5194/gmd-11-1229-2018, 2018. Wedi, N., Sandu, I., Bauer, P., Acosta, M., Andersen, R. C., Andrae, U., Auger, L., Balsamo, G., Baousis, V., Bennett, V., Bennett, A., Buontempo, C., Bretonnière, P.-A., Capell, R., Castrillo, M., Chantry, M., Chevallier, M., Correa, R., Davini, P., Denby, L., Doblas-Reyes, F., Dueben, P., Fischer, C., Frauen, C., Frogner, I.-L., Früh, B., Gascón, E., Gérard, E., Gorwits, O., Geenen, T., Grayson, K., GuenovaRubio, N., Hadade, I., von Hardenberg, J., Haus, U.-U., Hawkes, J., Hirtl, M., Hoffmann, J., Horvath, K., Järvinen, H., Jung, T., Kann, A., Klocke, D., Koldunov, N., Kontkanen, J., Sievi-Korte, O., Kristiansen, J., Kuwertz, E., Mäkelä, J., Maljutenko, I., Manninen, P., McKnight, U. S., Milinski, S., Mueller, A., McNally, A., Modigliani, U., Narayanappa, D., Nielsen, K. P., Nipen, T., Nortamo, H., Peuch, V.-H., Polade, S., Quintino, T., Schicker, I., Reuter, B., Smart, S., Sleigh, M., Suttie, M., Termonia, P., Thober, S., Randriamampianina,

36

R., Theeuwes, N., Thiemert, D., Vannière, B., Vannitsem, S., Wittmann, C., Yang, X., Pontaud, M., Stevens, B., and Pappenberger, F.: Implementing digital twin technology of the earth system in Destination Earth, Journal of the European Meteorological Society, 3, 100 015, https://doi.org/10.1016/j.jemets.2025.100015, 2025. Wunsch, C. and Heimbach, P.: Practical global oceanic state estimation, Physica D: Nonlinear Phenomena, 230, 197–208, https://doi.org/10.1016/j.physd.2006.09.040, 2007. Yan, F. E., Frezat, H., Le Sommer, J., Mak, J., and Otness, K.: Adjoint-based online learning of two-layer quasi-geostrophic baroclinic turbulence, Journal of Advances in Modeling Earth Systems, 17, e2024MS004 857, https://doi.org/10.1029/2024MS004857, 2025. Zalesak, S. T.: Fully multidimensional flux-corrected transport algorithms for fluids, Journal of Computational Physics, 31, 335–362, https://doi.org/10.1016/0021-9991(79)90051-2, 1979. Zhou, A., Hawkins, L., and Gentine, P.: Proof-of-concept: Using ChatGPT to translate and modernize an Earth system model from Fortran to Python/JAX, arXiv preprint arXiv:2405.00018, https://doi.org/10.48550/arXiv.2405.00018, 2024.

37

Related documents

• Sea-bed photographs (benthos) from the Marguarite Bay west of the Antarctic Peninsula along ROV profile JR37/015 (©AWI and BAS, Starmans, Clarke and Gutt 2003)#283681
• Sea-bed photographs (benthos) along profile PS56/111-2#241711
• High-resolution ocean and atmosphere pCO2 time-series measurements from mooring JKEO_147E_38N in the North Pacific Ocean (NCEI Accession 0100070)#408155
• Laboratory Evidence that Methanol-Rich Ice Mantles Lower Methyl Formate Binding Energies#1063242
• Mineral and oxygen isotopic compositions in a profile through the altered lower oceanic crust, ODP Hole 176-735B#660237