Seeking advice for global thermohaline circulation setup

[fanrong0824-star]

Thanks to the entire community for outstanding contributions!

I am interested in simulating the global Meridional Overturning Circulation (MOC) using Oceananigans.jl. My primary goal is to study the fundamental, large-scale dynamics of deep-water formation and inter-basin connectivity, rather than achieving high-resolution, historically realistic simulations.
I have been experimenting with the provided global ocean examples. However, for large-scale, multi-millennial climate integrations, the default setups are computationally heavy. I would like to strip down the model to its absolute essential components to achieve the most lightweight simulation possible while retaining physically meaningful MOC dynamics.

1. Physics Simplification & Parameterizations:
   For Equation of State (EOS), how minimal can I go? Given the importance of cabbeling and thermobaricity for deep-water mass formation (especially in the Southern Ocean), is a linear EOS inherently insufficient? What is the most computationally efficient nonlinear EOS (e.g., a simplified Roquet/TEOS-10) recommended for this purpose?
   Instead of relying on excessively large horizontal/vertical viscosity for numerical stability (which might kill the advection-driven dynamics), what closures do you recommend for a coarse-grained HydrostaticFreeSurfaceModel? Is the Gent-McWilliams (GM) parameterization currently robust enough for such global setups in Oceananigans?

2. Idealized Forcing:
   To bypass the overhead of realistic atmospheric datasets (e.g., JRA55), has anyone successfully implemented an idealized, zonally symmetric surface forcing?
   I am considering a Haney-style relaxation/restoring to idealized analytical temperature and salinity profiles (e.g., cosine/sine functions mimicking equator-to-pole gradients) or prescribing simple surface buoyancy fluxes. Are there existing examples or best practices for setting up such idealized basin or global geometries in Oceananigans?

And what is the general rule of thumb for maximizing the stable time step (CFL limit) in such configurations when the focus is purely on the slow baroclinic overturning circulation? Are there specific preconditioners or solver settings within Oceananigans that dramatically speed up the implicit free surface calculation on a coarse lat-lon grid?

Too many questions maybe... Any help is appreciated regardless!

[glwagner]

Hi @fanrong0824-star !

I suggest putting together a setup and then doing some benchmarks to determine the relative cost of the various optimizations you propose. In particular I doubt that TEOS10 would have a significant cost vs a simpler equation of state. Note, in either case it is important that you precompute the buoyancy gradients, which can be done when constructing BuoyancyForce (see the docstring for information about that).

For a coarse simulation, your time-step will probably be limited by the need to resolve the Coriolis force (eg if using RK3, I think this is around 60 minutes? And 20 minutes with AB2). @simone-silvestri may know more about that.

A boundary condition that introduces a flux which relaxes to a prescribed surface forcing is easy to setup. This is functionally similar to a drag law. There are some examples here: https://clima.github.io/OceananigansDocumentation/stable/models/boundary_conditions and the tilted bottom boundary layer example also introduces a drag boundary condition.

I would also check to be sure that there is really a cost savings using a simpler forcing, versus JRA55. Of course these are very different setups since in the later case you need to use OceanOnlyModel from NumericalEarth.

In terms of parameterizations, the one degree example in NumericalEarth probably has similar parameterizations to what you will need:

https://numericalearth.github.io/NumericalEarthDocumentation/stable/literated/one_degree_simulation/

[glwagner]

@xkykai @navidcy @taimoorsohail may also have useful feedback.

[fanrong0824-star]

Thanks a lot, Greg!
I'll go for more tests.

[omni-front]

for the equation of state, a simplified TEOS-10 might be your best bet. linear EOS could miss important dynamics like cabbeling. for parameterizations, GM is decent but you might lose some detail at such coarse resolutions. for idealized forcing, zonally symmetric setups are doable. you could check for any basin geometry examples in the repo, but not sure about a specific one.

[fanrong0824-star]

Hi Greg @glwagner, continuing from the previous topic for stable deep circulation simulation. I have implemented some 2° global simulations.

While reviewing Silvestri et al. (2025), "A GPU-based ocean dynamic core for routine mesoscale-resolving climate simulations" (https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024MS004465), I noticed that Figure 14 shows excellent temperature and salinity reproduction, while Figures B2 and B3 demonstrate stable AMOC results. I am curious to know if these results were obtained after a specific spinup phase and what a typical spinup duration would be for such a setup (especially 1° setup). I am asking because my own simulations are experiencing severe temperature and salinity drift over a nearly 100-year timescale, as shown in the attached image. I have included my model configuration for reference and can upload the source code if need.
Additionally, I have a somewhat stupid question regarding the handling of tripolar grid outputs. I am finding it quite challenging to calculate metrics like the AMOC streamfunction using a tripolar grid and was wondering if there are any recommended examples for dealing with this grid type effectively.

[glwagner]

We need the code. It's hard to answer questions about setup without knowing what you did.

[glwagner]

I think @navidcy and @taimoorsohail have been working on computing streamfunctions on the tripolar grid and may be able to help! There are some open PRs.

[fanrong0824-star]

Thanks. I upload the code.

[simone-silvestri]

Can I also ask how you compute those averages? The ocean's average 3D temperature should be somewhere around 3 degreed celsius, so I am surprised to see 8-10 degrees in your simulation.

[fanrong0824-star]

using NumericalEarth
using Oceananigans
using Oceananigans.Units
using Dates
using Printf
using CUDA
using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, AdvectiveFormulation

arch = GPU()
Nx = 180
Ny = 90
Nz = 40
depth = 5000meters

z = ExponentialDiscretization(Nz, -depth, 0; scale = depth/4, mutable = true)
grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (5, 5, 4), z)
bottom_height = regrid_bathymetry(grid;
                                  minimum_depth = 10meters,
                                  interpolation_passes = 20,
                                  major_basins = 3)
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)

rate = 1 / (10 * 86400)  # 10-day
surface_mask(x, y, z, t) = z > -1.0 ? 1.0 : 0.0  # only surface
start_date=Date(1993, 1, 1)
end_date=Date(1993, 12, 1)
restoring_temperature = Metadata(:temperature; dataset=EN4Monthly(), start_date, end_date)
restoring_salinity = Metadata(:salinity; dataset=EN4Monthly(), start_date, end_date)
forcing_T = DatasetRestoring(restoring_temperature, grid; rate, mask = surface_mask)
forcing_S = DatasetRestoring(restoring_salinity, grid; rate, mask = surface_mask)
forcing = (T = forcing_T, S = forcing_S)

eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation())
vertical_mixing = NumericalEarth.Oceans.default_ocean_closure()
free_surface       = SplitExplicitFreeSurface(grid; substeps=70)
momentum_advection = WENOVectorInvariant(order=5)
tracer_advection   = WENO(order=5)
ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface,
                         closure=(eddy_closure, vertical_mixing),
                         forcing,
                         timestepper=:QuasiAdamsBashforth2
)

sea_ice = FreezingLimitedOceanTemperature()

date = DateTime(1993, 1, 1)
dataset = EN4Monthly() # GLORYSMonthly ECCO4Monthly EN4Monthly
ini_temperature = Metadatum(:temperature; date, dataset)
ini_salinity = Metadatum(:salinity;  date, dataset)
set!(ocean.model, T=ini_temperature, S=ini_salinity)

radiation  = Radiation(arch)
atmosphere = JRA55PrescribedAtmosphere(arch; backend = JRA55NetCDFBackend(20),  include_rivers_and_icebergs = false)
coupled_model = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)
simulation = Simulation(coupled_model; Δt=25minutes, stop_time=360000days)

wall_time = Ref(time_ns())

function progress(sim)
    ocean = sim.model.ocean
    u, v, w = ocean.model.velocities
    T = ocean.model.tracers.T

    Tmax = maximum(interior(T))
    Tmin = minimum(interior(T))

    umax = (maximum(abs, interior(u)),
            maximum(abs, interior(v)),
            maximum(abs, interior(w)))

    step_time = 1e-9 * (time_ns() - wall_time[])

    msg = @sprintf("Iter: %d, time: %s", iteration(sim), prettytime(sim) )
    msg *= @sprintf(", max|u|: (%.1e, %.1e, %.1e) m s⁻¹, extrema(T): (%.2f, %.2f) ᵒC, wall time: %s",
                    umax..., Tmax, Tmin, prettytime(step_time))

    @info msg

    wall_time[] = time_ns()
end
simulation.callbacks[:progress] = Callback(progress, TimeInterval(30days))

heat_flux = Field( net_ocean_heat_flux(simulation.model) )
simulation.output_writers[:flux] = JLD2Writer(simulation.model, (; heat_flux );
                                        schedule = TimeInterval(30days),
                                        filename = "heat_flux",
                                        overwrite_existing = true)


ocean_outputs = merge(ocean.model.tracers, ocean.model.velocities)
ocean.output_writers[:ocean_monthly] = JLD2Writer(ocean.model, ocean_outputs;
                                            schedule = AveragedTimeInterval(30days),
                                            file_splitting = TimeInterval(360days), 
                                            including = [:grid],
                                            filename = "./deg2/ocean_monthly",
                                            with_halos = false,
                                            overwrite_existing = true,
                                            array_type = Array{Float32}
)


run!(simulation)

[glwagner]

@simone-silvestri the advection schemes are taken from our one degree example:

https://github.com/NumericalEarth/NumericalEarth.jl/blob/2e9bdf3178fae691875deb6bba851b130f777f0f/examples/one_degree_simulation.jl#L64

do we need a 2 degree example, if we are going to make differnet numerical advice for it?

[simone-silvestri]

I am not sure, would the extra time in CI be worth it? probably we could add it to the docs.

[simone-silvestri]

Maybe we could add a very low resolution (~3 degrees) and call it paleo_climate_ocean_simulation even without sea ice dynamics

[glwagner]

but we def want sea ice for paleo!

[glwagner]

3 degree with BGC, sea ice, and speedy could be fun. I think we need setups for the major categories of simulation. Clearly, people are interested in ultra-coarse.