Internal tide over a seamount
In this example, we show how an internal tide is generated from a barotropic tidal flow sloshing back and forth over a seamount.
Install dependencies
First let's make sure we have all required packages installed.
using Pkg
pkg"add Oceananigans, CairoMakie"using Oceananigans
using Oceananigans.UnitsGrid
We create an ImmersedBoundaryGrid wrapped around an underlying two-dimensional RectilinearGrid that is periodic in
Nx, Nz = 256, 128
H, L = 2kilometers, 1000kilometers
underlying_grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4),
x = (-L, L), z = (-H, 0),
topology = (Periodic, Flat, Bounded))256×1×128 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 4×0×4 halo
├── Periodic x ∈ [-1.0e6, 1.0e6) regularly spaced with Δx=7812.5
├── Flat y
└── Bounded z ∈ [-2000.0, 0.0] regularly spaced with Δz=15.625Now we can create the non-trivial bathymetry. We use PartialCellBottom that gets as input either (i) a two-dimensional function whose arguments are the grid's native horizontal coordinates and it returns the
In this example we'd like to have a Gaussian hill at the center of the domain.
h₀ = 250meters
σ = 20kilometers
bottom(x) = -H + h₀ * exp(-x^2 / 2σ^2)
grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(bottom))256×1×128 ImmersedBoundaryGrid{Float32, Periodic, Flat, Bounded} on CPU with 4×0×4 halo:
├── immersed_boundary: PartialCellBottom(mean(zb)=-1993.75, min(zb)=-2000.0, max(zb)=-1754.72, ϵ=0.2)
├── underlying_grid: 256×1×128 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 4×0×4 halo
├── Periodic x ∈ [-1.0e6, 1.0e6) regularly spaced with Δx=7812.5
├── Flat y
└── Bounded z ∈ [-2000.0, 0.0] regularly spaced with Δz=15.625Let's see what the domain with the bathymetry looks like.
x = xnodes(grid, Center())
bottom_height = interior(bottom_height_field(grid), :, 1, 1)
using CairoMakie
fig = Figure(size = (700, 200))
ax = Axis(fig[1, 1],
xlabel = "x [km]",
ylabel = "z [m]",
limits = ((-L, L) ./ kilometers, (-H, 0)))
band!(ax, x / kilometers, bottom_height, zero(x), color = :mediumblue)
fig
Now we want to add a barotropic tide forcing. For example, to add the lunar semi-diurnal
where
The excursion parameter is a nondimensional number that expresses the ratio of the tidal excursion to the width of the hill,
We prescribe the excursion parameter which, in turn, implies a tidal velocity
Now we have a way to find the tidal forcing amplitude that corresponds to a given excursion parameter. The Coriolis frequency is needed, so we start by constructing an
coriolis = FPlane(latitude = -45)
f = coriolis.f-0.000103126076f0Now we have everything we require to construct the tidal forcing given a value of the excursion parameter.
T₂ = 12.421hours
ω₂ = 2π / T₂
ϵ = 0.1
U₂ = ϵ * ω₂ * σ
F₀ = U₂ * (ω₂^2 - f^2) / ω₂
@inline tidal_forcing(x, z, t, p) = p.F₀ * sin(p.ω₂ * t)
u_forcing = Forcing(tidal_forcing, parameters=(; F₀, ω₂))ContinuousForcing{@NamedTuple{F₀::Float64, ω₂::Float64}}
├── func: tidal_forcing (generic function with 1 method)
├── parameters: (F₀ = 1.8218612216084723e-5, ω₂ = 0.00014051439111137024)
└── field dependencies: ()Model
We build a HydrostaticFreeSurfaceModel:
model = HydrostaticFreeSurfaceModel(grid; coriolis,
buoyancy = BuoyancyTracer(),
tracers = :b,
momentum_advection = WENO(),
tracer_advection = WENO(),
forcing = (; u = u_forcing))HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── grid: 256×1×128 ImmersedBoundaryGrid{Float32, Periodic, Flat, Bounded} on CPU with 4×0×4 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: b
├── closure: Nothing
├── buoyancy: BuoyancyTracer with ĝ = NegativeZDirection()
├── free surface: SplitExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
│ └── substepping: FixedTimeStepSize(39.049 seconds)
├── advection scheme:
│ ├── momentum: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
│ └── b: WENO{3, Float32, Oceananigans.Utils.BackendOptimizedDivision}(order=5)
├── vertical_coordinate: ZCoordinate
└── coriolis: FPlane{Oceananigans.Advection.EnstrophyConserving{Float32}, Float32}We initialize the model with the tidal flow and a linear stratification.
Nᵢ² = 1e-4 # s⁻²
bᵢ(x, z) = Nᵢ² * z
set!(model, u=U₂, b=bᵢ)Now let's build a Simulation.
simulation = Simulation(model, Δt=5minutes, stop_time=4days)Simulation of HydrostaticFreeSurfaceModel{CPU, ImmersedBoundaryGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 5 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: 4 days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│ └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entriesWe add a callback to print a message about how the simulation is going,
using Printf
progress(sim) = @info @sprintf("Iter: %d, time: %s, wall time: %s, max|w|: %6.3e m s⁻¹",
iteration(sim), prettytime(sim), prettytime(sim.run_wall_time),
maximum(abs, sim.model.velocities.w))
add_callback!(simulation, progress, IterationInterval(200))Diagnostics/Output
Add some diagnostics. Instead of
b = model.tracers.b
u, v, w = model.velocities
U = Field(Average(u))
u′ = u - U
N² = ∂z(b)
filename = "internal_tide.jld2"
simulation.output_writers[:fields] = JLD2Writer(model, (; u, u′, w, b, N²); filename,
schedule = TimeInterval(30minutes),
overwrite_files = true)JLD2Writer scheduled on TimeInterval(30 minutes):
├── filepath: internal_tide.jld2
├── 5 outputs: (u, u′, w, b, N²)
├── array_type: Array{Float32}
├── including: [:coriolis, :buoyancy, :closure]
├── file_splitting: NoFileSplitting
└── file size: 0 bytes (file not yet created)We are ready – let's run!
run!(simulation)[ Info: Initializing simulation...
[ Info: Iter: 0, time: 0 seconds, wall time: 0 seconds, max|w|: 2.066e-03 m s⁻¹
[ Info: ... simulation initialization complete (4.943 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (3.154 seconds).
[ Info: Iter: 200, time: 16.667 hours, wall time: 10.657 seconds, max|w|: 5.711e-03 m s⁻¹
[ Info: Iter: 400, time: 1.389 days, wall time: 13.460 seconds, max|w|: 5.070e-03 m s⁻¹
[ Info: Iter: 600, time: 2.083 days, wall time: 16.181 seconds, max|w|: 2.500e-03 m s⁻¹
[ Info: Iter: 800, time: 2.778 days, wall time: 19.157 seconds, max|w|: 5.701e-03 m s⁻¹
[ Info: Iter: 1000, time: 3.472 days, wall time: 21.792 seconds, max|w|: 3.912e-03 m s⁻¹
[ Info: Simulation is stopping after running for 23.968 seconds.
[ Info: Simulation time 4 days equals or exceeds stop time 4 days.Load output
First, we load the saved velocities and stratification output as FieldTimeSerieses.
u′_timeseries = FieldTimeSeries(filename, "u′")
w_timeseries = FieldTimeSeries(filename, "w")
N²_timeseries = FieldTimeSeries(filename, "N²")
u′max = maximum(abs, u′_timeseries[end])
wmax = maximum(abs, w_timeseries[end])
times = u′_timeseries.timesVisualize
Now we can visualize our results! We use CairoMakie here. On a system with OpenGL using GLMakie is more convenient as figures will be displayed on the screen.
We use Makie's Observable to animate the data. To dive into how Observables work we refer to Makie.jl's Documentation.
n = Observable(1)
title = @lift @sprintf("t = %1.2f days = %1.2f T₂", times[$n] / day, times[$n] / T₂)
u′ₙ = @lift u′_timeseries[$n]
wₙ = @lift w_timeseries[$n]
N²ₙ = @lift N²_timeseries[$n]
axis_kwargs = (xlabel = "x [m]",
ylabel = "z [m]",
limits = ((-L, L), (-H, 0)),
titlesize = 20)
fig = Figure(size = (700, 900))
fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false)
ax_u = Axis(fig[2, 1]; title = "u'-velocity", axis_kwargs...)
hm_u = heatmap!(ax_u, u′ₙ; nan_color=:gray, colorrange=(-u′max, u′max), colormap=:balance)
Colorbar(fig[2, 2], hm_u, label = "m s⁻¹")
ax_w = Axis(fig[3, 1]; title = "w-velocity", axis_kwargs...)
hm_w = heatmap!(ax_w, wₙ; nan_color=:gray, colorrange=(-wmax, wmax), colormap=:balance)
Colorbar(fig[3, 2], hm_w, label = "m s⁻¹")
ax_N² = Axis(fig[4, 1]; title = "stratification N²", axis_kwargs...)
hm_N² = heatmap!(ax_N², N²ₙ; nan_color=:gray, colorrange=(0.9Nᵢ², 1.1Nᵢ²), colormap=:magma)
Colorbar(fig[4, 2], hm_N², label = "s⁻²")
fig
Finally, we can record a movie.
@info "Making an animation from saved data..."
frames = 1:length(times)
record(fig, "internal_tide.mp4", frames, framerate=16) do i
n[] = i
end[ Info: Making an animation from saved data...Julia version and environment information
This example was executed with the following version of Julia:
using InteractiveUtils: versioninfo
versioninfo()Julia Version 1.13.0
Commit d1c37793dd2 (2026-09-09 19:00 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 128 × AMD EPYC 9374F 32-Core Processor
WORD_SIZE: 64
LLVM: libLLVM-20.1.8 (ORCJIT, znver4)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 128 virtual cores)
Environment:
JULIA_LOAD_PATH = @:@v#.#:@stdlib
JULIA_DEPOT_PATH = /var/lib/buildkite-agent/.julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/local/share/julia:/var/lib/buildkite-agent/.julia/juliaup/julia-1.13.0+0.x64.linux.gnu/share/julia
JULIA_VERSION_ENZYME = 1.13.0
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
LD_LIBRARY_PATH =
JULIA_MAX_NUM_PRECOMPILE_FILES = 24
JULIA_VERSION = 1.13.0
JULIA_CUDA_USE_COMPAT = false
JULIA_PROJECT = /var/lib/buildkite-agent/Oceananigans.jl-34931/docs/
JULIA_DEBUG = LiterateThese were the top-level packages installed in the environment:
import Pkg
Pkg.status()Status `~/Oceananigans.jl-34931/docs/Project.toml`
[79e6a3ab] Adapt v4.7.4
⌃ [052768ef] CUDA v6.2.2
[13f3f980] CairoMakie v0.15.16
[8e50ac2c] ConservativeRegridding v0.2.16
⌅ [e30172f5] Documenter v1.17.0
[daee34ce] DocumenterCitations v1.5.0
[4710194d] DocumenterVitepress v0.3.8
[7da242da] Enzyme v0.13.216
[033835bb] JLD2 v0.6.7
⌅ [63c18a36] KernelAbstractions v0.9.45
[98b081ad] Literate v2.21.0
[da04e1cc] MPI v0.20.29
[85f8d34a] NCDatasets v0.14.15
[9e8cae18] Oceananigans v0.114.0 `..`
[429524aa] Optim v2.3.2
[f27b6e38] Polynomials v4.1.3
⌅ [3c362404] Reactant v0.2.292
[6038ab10] Rotations v1.7.1
[d496a93d] SeawaterPolynomials v0.3.11
[09ab397b] StructArrays v0.7.4
[bdfc003b] TimesDates v0.3.3
[0a941bbe] Zarr v0.11.2
[b77e0a4c] InteractiveUtils v1.11.0
[37e2e46d] LinearAlgebra v1.13.0
[44cfe95a] Pkg v1.13.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`This page was generated using Literate.jl.