Global climate simulation

This example configures a global ocean–sea ice simulation at 1.5ᵒ horizontal resolution with realistic bathymetry and a few closures including the "Gent-McWilliams" IsopycnalSkewSymmetricDiffusivity. The atmosphere is represented by a 4-layer SpeedyWeather simulation on the T63 spectral grid (this grid has approximately 1.875ᵒ resolution).

The atmosphere is initialized with the Jablonowski and Williamson (2006) initial conditions, which consist of a zonal wind centered at mid-latitudes and higher altitudes and a temperature profile that is baroclinically unstable. The surface pressure is adjusted by orography for approximately globally constant mean-sea level pressure. The initial specific humidity is calculated from temperature for a constant relative humidity everywhere. The ocean and sea ice are initialized by ocean temperature, salinity, sea ice concentration, and sea ice thickness from the ECCO state estimate.

For this example, we need Oceananigans.HydrostaticFreeSurfaceModel (the ocean), ClimaSeaIce.SeaIceModel (the sea ice) and SpeedyWeather.PrimitiveWetModel (the atmosphere). All these are coupled and orchestrated by the NumericalEarth.EarthSystemModel (the coupled system).

The ConservativeRegridding.jl package is used to regrid fields between the atmosphere and ocean–sea ice components.

using Oceananigans, SpeedyWeather, NumericalEarth, ConservativeRegridding
using NCDatasets, CairoMakie
using Oceananigans.Units
using Printf, Statistics, Dates

[ Info: Oceananigans will use 8 threads

Ocean and sea-ice model configuration

The ocean and sea-ice are a simplified versions of the one-degree ocean-sea ice example.

The first step is to create the grid with realistic bathymetry.

Nx = 240
Ny = 120
Nz = 10
z = ExponentialDiscretization(Nz, -2000, 0)
grid = TripolarGrid(Oceananigans.CPU(); size=(Nx, Ny, Nz), z, halo=(6, 6, 5))

We regrid the bathymetry.

bottom_height = regrid_bathymetry(grid; major_basins=1, interpolation_passes=15)
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom_height); active_cells_map=true)

[ Info: Loading cached bathymetry from /storage5/github-action-runners/julia-depot/scratchspaces/904d977b-046a-4731-8b86-9235c0d1ef02/bathymetry_cache/bathymetry_240x120_0.0_360.0_-89.28571566522153_90.0_5a3b1e40.jld2

Now we can specify the numerical details and the closures for the ocean simulation.

momentum_advection   = VectorInvariant()
tracer_advection     = WENO(order=5)
free_surface         = SplitExplicitFreeSurface(grid; substeps=40)
catke_closure        = NumericalEarth.Oceans.default_ocean_closure()
eddy_closure         = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3)
viscous_closure      = Oceananigans.TurbulenceClosures.HorizontalScalarBiharmonicDiffusivity(ν=1e12)
closures             = (catke_closure, eddy_closure, viscous_closure)

The ocean simulation, complete with initial conditions for temperature and salinity from ECCO on Jan 1st, 1992.

ocean = ocean_simulation(grid;
                         momentum_advection,
                         tracer_advection,
                         free_surface,
                         closure = closures)

ecco_set = MetadataSet(:temperature, :salinity,
                       :sea_ice_thickness, :sea_ice_concentration;
                       dataset = ECCO4Monthly(),
                       date = DateTime(1992, 1, 1))
Oceananigans.set!(ocean.model, ecco_set)   # T, S

The sea-ice simulation, complete with initial conditions for sea-ice thickness and sea-ice concentration from ECCO.

sea_ice = sea_ice_simulation(grid, ocean; advection=WENO(order=5))

Oceananigans.set!(sea_ice.model, ecco_set)   # h, ℵ

Atmosphere model configuration

The atmosphere is provided by SpeedyWeather.jl. Here, we configure a T63L4 model with a 3-hour output interval. The atmosphere_simulation function takes care of building an atmosphere model with appropriate hooks so that NumericalEarth can compute inter-component fluxes.

nlayers = 4
spectral_grid = SpeedyWeather.SpectralGrid(; trunc=63, nlayers, Grid=FullClenshawGrid)
atmosphere = atmosphere_simulation(spectral_grid; output_interval=Hour(3))

Simulation{...} (16.14 MB)
├ variables::Variables{...} (13.00 MB)
└ model::SpeedyWeather.PrimitiveWetModel{...} (3.62 MB)

The atmosphere model already includes some initial conditions as described above:

atmosphere.model.initial_conditions

(vordiv = SpeedyWeather.ZonalWind{Float32} <: SpeedyWeather.AbstractInitialConditions
├ η₀::Float32 = 0.252
├ u₀::Float32 = 35.0
├ perturb_lat::Float32 = 40.0
├ perturb_lon::Float32 = 20.0
├ perturb_uₚ::Float32 = 1.0
└ perturb_radius::Float32 = 0.1, pres = SpeedyWeather.PressureOnOrography <: SpeedyWeather.AbstractInitialConditions
, temp = SpeedyWeather.JablonowskiTemperature{Float32} <: SpeedyWeather.AbstractInitialConditions
├ η₀::Float32 = 0.252
├ σ_tropopause::Float32 = 0.2
├ u₀::Float32 = 35.0
├ ΔT::Float32 = 0.0
└ Tmin::Float32 = 200.0, humid = SpeedyWeather.ConstantRelativeHumidity{Float32} <: SpeedyWeather.AbstractInitialConditions
└ relhumid_ref::Float32 = 0.7)

The coupled model

We are now ready to blend everything together. We need to specify the time step for the coupled model. We decide to step the global model every 2 atmosphere time steps (i.e., the ocean and the sea-ice are stepped every two atmosphere time steps).

Δt = 2 * convert(eltype(grid), atmosphere.model.time_stepping.Δt_sec)

We build the complete coupled earth_model and the coupled simulation. NumericalEarth still computes turbulent (sensible heat, water vapor) fluxes here using its own bulk-formula machinery and writes them back to Speedy. Top-level radiation, however, is not yet wired up against SpeedyWeather's upwelling-LW path, so we leave radiation = nothing for now.

earth_model = EarthSystemModel(; atmosphere, sea_ice, ocean)

EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: SpeedyWeather.Simulation
├── land: Nothing
├── sea_ice: SeaIceModel
├── ocean: HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfaces

Building and running the simulation

We are ready to build and run the coupled simulation. But before we do, we add callbacks to write outputs to disk every 3 hours.

earth = Oceananigans.Simulation(earth_model; Δt, stop_time=30days)
outputs = merge(ocean.model.velocities, ocean.model.tracers)
sea_ice_fields = merge(sea_ice.model.velocities, sea_ice.model.dynamics.auxiliaries.fields,
                       (; h=sea_ice.model.ice_thickness, ℵ=sea_ice.model.ice_concentration))

ocean.output_writers[:free_surf] = JLD2Writer(ocean.model, (; η=ocean.model.free_surface.displacement);
                                              overwrite_existing=true,
                                              schedule=TimeInterval(3hours),
                                              filename="ocean_free_surface.jld2")

ocean.output_writers[:surface] = JLD2Writer(ocean.model, outputs;
                                            overwrite_existing=true,
                                            schedule=TimeInterval(3hours),
                                            filename="ocean_surface_fields.jld2",
                                            indices=(:, :, grid.Nz))

sea_ice.output_writers[:fields] = JLD2Writer(sea_ice.model, sea_ice_fields;
                                             overwrite_existing=true,
                                             schedule=TimeInterval(3hours),
                                             filename="sea_ice_fields.jld2")

𝒬ᵀᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.sensible_heat
𝒬ᵛᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.latent_heat
τˣᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.x_momentum
τʸᵃᵒ = earth.model.interfaces.atmosphere_ocean_interface.fluxes.y_momentum
𝒬ᵀᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.sensible_heat
𝒬ᵛᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.latent_heat
τˣᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.x_momentum
τʸᵃⁱ = earth.model.interfaces.atmosphere_sea_ice_interface.fluxes.y_momentum
𝒬ⁱᵒ  = earth.model.interfaces.net_fluxes.sea_ice.bottom.heat
Jˢⁱᵒ  = earth.model.interfaces.sea_ice_ocean_interface.fluxes.salt
fluxes = (; 𝒬ᵀᵃᵒ, 𝒬ᵛᵃᵒ, τˣᵃᵒ, τʸᵃᵒ, 𝒬ᵀᵃⁱ, 𝒬ᵛᵃⁱ, τˣᵃⁱ, τʸᵃⁱ, 𝒬ⁱᵒ, Jˢⁱᵒ)

ocean.output_writers[:fluxes] = JLD2Writer(earth.model.ocean.model, fluxes;
                                           overwrite_existing=true,
                                           schedule=TimeInterval(3hours),
                                           filename="intercomponent_fluxes.jld2")

JLD2Writer scheduled on TimeInterval(3 hours):
├── filepath: intercomponent_fluxes.jld2
├── 10 outputs: (𝒬ᵀᵃᵒ, 𝒬ᵛᵃᵒ, τˣᵃᵒ, τʸᵃᵒ, 𝒬ᵀᵃⁱ, 𝒬ᵛᵃⁱ, τˣᵃⁱ, τʸᵃⁱ, 𝒬ⁱᵒ, Jˢⁱᵒ)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)

We also add a callback function that prints out a helpful progress message while the simulation runs.

wall_time = Ref(time_ns())

function progress(sim)
    atmos = sim.model.atmosphere
    ocean = sim.model.ocean

    ua, va     = atmos.variables.dynamics.u_mean_grid, atmos.variables.dynamics.v_mean_grid
    uo, vo, wo = ocean.model.velocities

    uamax = (maximum(abs, ua), maximum(abs, va))
    uomax = (maximum(abs, uo), maximum(abs, vo), maximum(abs, wo))

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

    msg1 = @sprintf("time: %s, iter: %d", prettytime(sim), iteration(sim))
    msg2 = @sprintf(", max|ua|: (%.1e, %.1e) m s⁻¹", uamax...)
    msg3 = @sprintf(", max|uo|: (%.1e, %.1e, %.1e) m s⁻¹", uomax...)
    msg4 = @sprintf(", wall time: %s \n", prettytime(step_time))

    @info msg1 * msg2 * msg3 * msg4

    wall_time[] = time_ns()

     return nothing
end

add_callback!(earth, progress, TimeInterval(2days))

Let's run the coupled model!

Oceananigans.run!(earth)

[ Info: Initializing simulation...
[ Info: time: 0 seconds, iter: 0, max|ua|: (0.0e+00, 0.0e+00) m s⁻¹, max|uo|: (0.0e+00, 0.0e+00, 0.0e+00) m s⁻¹, wall time: 2.782 minutes 
[ Info:     ... simulation initialization complete (6.448 seconds)
[ Info: Executing initial time step...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters /storage5/github-action-runners/julia-depot/packages/Oceananigans/pztwa/src/OutputWriters/jld2_writer.jl:234
[ Info:     ... initial time step complete (19.119 minutes).
[ Info: time: 2 days, iter: 72, max|ua|: (2.9e+01, 8.5e+00) m s⁻¹, max|uo|: (9.1e-01, 1.3e+00, 4.8e-03) m s⁻¹, wall time: 23.554 minutes 
[ Info: time: 4 days, iter: 144, max|ua|: (3.3e+01, 1.9e+01) m s⁻¹, max|uo|: (7.9e-01, 1.6e+00, 8.7e-03) m s⁻¹, wall time: 4.503 minutes 
[ Info: time: 6 days, iter: 216, max|ua|: (3.3e+01, 2.0e+01) m s⁻¹, max|uo|: (9.3e-01, 1.3e+00, 4.8e-03) m s⁻¹, wall time: 4.474 minutes 
[ Info: time: 8 days, iter: 288, max|ua|: (3.6e+01, 1.8e+01) m s⁻¹, max|uo|: (1.0e+00, 1.2e+00, 8.2e-03) m s⁻¹, wall time: 4.522 minutes 
[ Info: time: 10 days, iter: 360, max|ua|: (3.9e+01, 1.9e+01) m s⁻¹, max|uo|: (1.0e+00, 9.3e-01, 6.7e-03) m s⁻¹, wall time: 4.465 minutes 
[ Info: time: 12 days, iter: 432, max|ua|: (3.8e+01, 2.1e+01) m s⁻¹, max|uo|: (1.0e+00, 9.4e-01, 6.5e-03) m s⁻¹, wall time: 4.470 minutes 
[ Info: time: 14 days, iter: 504, max|ua|: (4.2e+01, 2.0e+01) m s⁻¹, max|uo|: (1.4e+00, 1.0e+00, 6.1e-03) m s⁻¹, wall time: 4.477 minutes 
[ Info: time: 16 days, iter: 576, max|ua|: (4.2e+01, 2.2e+01) m s⁻¹, max|uo|: (1.4e+00, 1.3e+00, 5.8e-03) m s⁻¹, wall time: 4.622 minutes 
[ Info: time: 18 days, iter: 648, max|ua|: (4.2e+01, 2.5e+01) m s⁻¹, max|uo|: (1.4e+00, 1.3e+00, 3.3e-03) m s⁻¹, wall time: 4.478 minutes 
[ Info: time: 20 days, iter: 720, max|ua|: (3.9e+01, 2.3e+01) m s⁻¹, max|uo|: (1.3e+00, 1.3e+00, 3.2e-03) m s⁻¹, wall time: 4.524 minutes 
[ Info: time: 22 days, iter: 792, max|ua|: (3.8e+01, 2.4e+01) m s⁻¹, max|uo|: (1.5e+00, 1.2e+00, 2.3e-03) m s⁻¹, wall time: 4.503 minutes 
[ Info: time: 24 days, iter: 864, max|ua|: (4.0e+01, 1.7e+01) m s⁻¹, max|uo|: (1.7e+00, 1.3e+00, 1.7e-03) m s⁻¹, wall time: 4.505 minutes 
[ Info: time: 26 days, iter: 936, max|ua|: (4.4e+01, 1.7e+01) m s⁻¹, max|uo|: (1.6e+00, 1.2e+00, 1.8e-03) m s⁻¹, wall time: 4.487 minutes 
[ Info: time: 28 days, iter: 1008, max|ua|: (4.2e+01, 2.8e+01) m s⁻¹, max|uo|: (1.7e+00, 1.2e+00, 2.2e-03) m s⁻¹, wall time: 4.513 minutes 
[ Info: Simulation is stopping after running for 1.444 hours.
[ Info: Simulation time 30 days equals or exceeds stop time 30 days.
[ Info: time: 30 days, iter: 1080, max|ua|: (3.9e+01, 2.6e+01) m s⁻¹, max|uo|: (1.6e+00, 1.2e+00, 2.4e-03) m s⁻¹, wall time: 4.511 minutes 

Visualizing the results

We can visualize some of the results. Here, we plot the surface speeds in the atmosphere, ocean, and sea-ice as well as the 2-meter temperature in the atmosphere, the sea surface temperature, and the sensible and latent heat fluxes at the atmosphere-ocean interface. SpeedyWeather outputs are stored in a NetCDF file located in the run_0001 folder, while ocean and sea-ice outputs are stored in JLD2 files that can be read by Oceananigans.jl using the FieldTimeSeries type.

SWO = Dataset("run_0001/output.nc")

Ta = reverse(SWO["temp"][:, :, nlayers, :], dims=2)
ua = reverse(SWO["u"][:, :, nlayers, :],    dims=2)
va = reverse(SWO["v"][:, :, nlayers, :],    dims=2)
sp = sqrt.(ua.^2 + va.^2)

SST = FieldTimeSeries("ocean_surface_fields.jld2", "T")
SSU = FieldTimeSeries("ocean_surface_fields.jld2", "u")
SSV = FieldTimeSeries("ocean_surface_fields.jld2", "v")

SIU = FieldTimeSeries("sea_ice_fields.jld2", "u")
SIV = FieldTimeSeries("sea_ice_fields.jld2", "v")
SIA = FieldTimeSeries("sea_ice_fields.jld2", "ℵ")

𝒬ᵀᵃᵒ = FieldTimeSeries("intercomponent_fluxes.jld2", "𝒬ᵀᵃᵒ")
𝒬ᵛᵃᵒ = FieldTimeSeries("intercomponent_fluxes.jld2", "𝒬ᵛᵃᵒ")

Nt = min(length(sp[1, 1, :]), length(𝒬ᵀᵃᵒ))
times = 𝒬ᵀᵃᵒ.times

uotmp = Oceananigans.Field{Face, Center, Nothing}(SST.grid)
votmp = Oceananigans.Field{Center, Face, Nothing}(SST.grid)

uitmp = Oceananigans.Field{Face,   Center, Nothing}(SST.grid)
vitmp = Oceananigans.Field{Center, Face,   Nothing}(SST.grid)
atmp  = Oceananigans.Field{Center, Center, Nothing}(SST.grid)

sotmp = Oceananigans.Field(sqrt(uotmp^2 + votmp^2))
sitmp = Oceananigans.Field(sqrt(uitmp^2 + vitmp^2) * atmp)

iter = Observable(1)

san = @lift sp[:, :, $iter]
son  = @lift begin
    Oceananigans.set!(uotmp, SSU[$iter])
    Oceananigans.set!(votmp, SSV[$iter])
    Oceananigans.compute!(sotmp)
    Oceananigans.interior(sotmp, :, :, 1)
end
ssn  = @lift begin
    Oceananigans.set!(uitmp, SIU[$iter])
    Oceananigans.set!(vitmp, SIV[$iter])
    Oceananigans.set!(atmp,  SIA[$iter])
    Oceananigans.compute!(sitmp)
    Oceananigans.interior(sitmp, :, :, 1)
end

fig = Figure(size = (1000, 1500))

ax1 = Axis(fig[1, 1], title = "Surface speed, atmosphere")
ax2 = Axis(fig[2, 1], title = "Surface speed, ocean")
ax3 = Axis(fig[3, 1], title = "Surface speed, sea-ice")

hm1 = heatmap!(ax1, san; colormap = :deep,  nan_color=:lightgray, colorrange = (0, 35))
hm2 = heatmap!(ax2, son; colormap = :magma, nan_color=:lightgray, colorrange = (0, 0.6))
hm3 = heatmap!(ax3, ssn; colormap = :ice,   nan_color=:lightgray, colorrange = (0, 0.6))

Colorbar(fig[1, 2], hm1, label="(m s⁻¹)")
Colorbar(fig[2, 2], hm2, label="(m s⁻¹)")
Colorbar(fig[3, 2], hm3, label="(m s⁻¹)")

for ax in (ax1, ax2, ax3)
    hidedecorations!(ax)
end

title = @lift string(prettytime(times[$iter] - times[1]))
Label(fig[0, :], title, fontsize=18)

record(fig, "surface_speeds.mp4", 1:Nt, framerate=8) do i
    iter[] = i
end

Tan = @lift Ta[:, :, $iter]
Ton = @lift interior(SST[$iter], :, :, 1)
𝒬ᵀn = @lift interior(𝒬ᵀᵃᵒ[$iter], :, :, 1)
𝒬ᵛn = @lift interior(𝒬ᵛᵃᵒ[$iter], :, :, 1)

fig = Figure(size = (1000, 2000))

ax1 = Axis(fig[1, 1], title = "2m Temperature, atmosphere")
ax2 = Axis(fig[2, 1], title = "Sea Surface Temperature")
ax3 = Axis(fig[3, 1], title = "Sensible heat flux")
ax4 = Axis(fig[4, 1], title = "Latent heat flux")

hm1 = heatmap!(ax1, Tan; colormap = :plasma, nan_color=:lightgray, colorrange = (-45, 30))
hm2 = heatmap!(ax2, Ton; colormap = :plasma, nan_color=:lightgray, colorrange = (-2, 32))
hm3 = heatmap!(ax3, 𝒬ᵀn; colormap = :balance, colorrange = (-200, 200),  nan_color=:lightgray)
hm4 = heatmap!(ax4, 𝒬ᵛn; colormap = :balance, colorrange = (-200, 200),  nan_color=:lightgray)

Colorbar(fig[1, 2], hm1, label="(ᵒC)")
Colorbar(fig[2, 2], hm2, label="(ᵒC)")
Colorbar(fig[3, 2], hm3, label="(W/m²)")
Colorbar(fig[4, 2], hm4, label="(W/m²)")

for ax in (ax1, ax2, ax3, ax4)
    hidedecorations!(ax)
end

title = @lift string(prettytime(times[$iter] - times[1]))
Label(fig[0, :], title, fontsize=18)

record(fig, "surface_temperature_and_heat_flux.mp4", 1:Nt, framerate=8) do i
    iter[] = i
end

This page was generated using Literate.jl.