Internal wave example

In this example, we initialize an internal wave packet in two-dimensions and watch it propagate.

using Oceananigans, Plots, Printf

┌ Debug: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1242
┌ Debug: Precompiling GeometryTypes [4d00f742-c7ba-57c2-abde-4428a4b178cb]
└ @ Base loading.jl:1242
┌ Debug: Precompiling ColorTypes [3da002f7-5984-5a60-b8a6-cbb66c0b333f]
└ @ Base loading.jl:1242
┌ Debug: Precompiling FixedPointNumbers [53c48c17-4a7d-5ca2-90c5-79b7896eea93]
└ @ Base loading.jl:1242
┌ Debug: Precompiling FFMPEG [c87230d0-a227-11e9-1b43-d7ebe4e7570a]
└ @ Base loading.jl:1242
┌ Debug: Precompiling RecipesBase [3cdcf5f2-1ef4-517c-9805-6587b60abb01]
└ @ Base loading.jl:1242
┌ Debug: Precompiling PlotUtils [995b91a9-d308-5afd-9ec6-746e21dbc043]
└ @ Base loading.jl:1242
┌ Debug: Precompiling Colors [5ae59095-9a9b-59fe-a467-6f913c188581]
└ @ Base loading.jl:1242
┌ Debug: Precompiling PlotThemes [ccf2f8ad-2431-5c83-bf29-c5338b663b6a]
└ @ Base loading.jl:1242
┌ Debug: Precompiling Showoff [992d4aef-0814-514b-bc4d-f2e9a6c4116f]
└ @ Base loading.jl:1242
┌ Debug: Precompiling StatsBase [2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91]
└ @ Base loading.jl:1242
┌ Debug: Precompiling SortingAlgorithms [a2af1166-a08f-5f64-846c-94a0d3cef48c]
└ @ Base loading.jl:1242
┌ Debug: Precompiling NaNMath [77ba4419-2d1f-58cd-9bb1-8ffee604a2e3]
└ @ Base loading.jl:1242
┌ Debug: Precompiling Measures [442fdcdd-2543-5da2-b0f3-8c86c306513e]
└ @ Base loading.jl:1242
┌ Debug: Precompiling GR [28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71]
└ @ Base loading.jl:1242
.

Numerical, domain, and internal wave parameters

First, we pick some numerical and physical parameters for our model and its rotation rate.

Nx = 128 # resolution
Lx = 2π  # domain extent

We set up an internal wave with the pressure field

$ p(x, y, z, t) = a(x, z) cos(kx + mz - ω t) $.

where m is the vertical wavenumber, k is the horizontal wavenumber, ω is the wave frequncy, and a(x, z) is a Gaussian envelope.

# Non-dimensional internal wave parameters
m = 16      # vertical wavenumber
k = 1       # horizontal wavenumber
N = 1       # buoyancy frequency
f = 0.2     # inertial frequency

A Gaussian wavepacket

Next, we set up an initial condition corresponding to a propagating wave packet with a Gaussian envelope. The internal wave dispersion relation yields

ω² = (N^2 * k^2 + f^2 * m^2) / (k^2 + m^2)

# and thus
ω = sqrt(ω²)

The internal wave polarization relations follow from the linearized Boussinesq equations,

U = k * ω   / (ω^2 - f^2)
V = k * f   / (ω^2 - f^2)
W = m * ω   / (ω^2 - N^2)
B = m * N^2 / (ω^2 - N^2)

Finally, we set-up a small-amplitude, Gaussian envelope for the wave packet

# Some Gaussian parameters
A, x₀, z₀, δ = 1e-9, Lx/2, -Lx/2, Lx/15

# A Gaussian envelope
a(x, z) = A * exp( -( (x - x₀)^2 + (z - z₀)^2 ) / 2δ^2 )

Create initial condition functions

u₀(x, y, z) = a(x, z) * U * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * V * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * W * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z

We are now ready to instantiate our model on a uniform grid. We give the model a constant rotation rate with background vorticity f, use temperature as a buoyancy tracer, and use a small constant viscosity and diffusivity to stabilize the model.

model = Model(
        grid = RegularCartesianGrid(size=(Nx, 1, Nx), length=(Lx, Lx, Lx)),
     closure = ConstantIsotropicDiffusivity(ν=1e-6, κ=1e-6),
    coriolis = FPlane(f=f),
     tracers = :b,
    buoyancy = BuoyancyTracer()
)

We initialize the velocity and buoyancy fields with our internal wave initial condition.

set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

A wave packet on the loose

Finally, we release the packet and watch it go!

anim = @animate for i=1:100
    time_step!(model, Nt = 20, Δt = 0.001 * 2π/ω)

    x, z = model.grid.xC, model.grid.zF
    w = model.velocities.w
    heatmap(x, z, w.data[1:Nx, 1, 1:Nx+1]', title=@sprintf("t = %.2f", model.clock.time),
            xlabel="x", ylabel="z", c=:balance, clims=(-1e-8, 1e-8))
end

This page was generated using Literate.jl.