Geostrophic adjustment using Oceananigans.HydrostaticFreeSurfaceModel

This example demonstrates how to simulate the one-dimensional geostrophic adjustment of a free surface using Oceananigans.HydrostaticFreeSurfaceModel. Here, we solve the hydrostatic Boussinesq equations beneath a free surface with a small-amplitude about rest $z = 0$, with boundary conditions expanded around $z = 0$, and free surface dynamics linearized under the assumption $η / H \ll 1$, where $η$ is the free surface displacement, and $H$ is the total depth of the fluid.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, JLD2, Plots"

A one-dimensional domain

We use a one-dimensional domain of geophysical proportions,

using Oceananigans
using Oceananigans.Utils: kilometers

grid = RegularRectilinearGrid(size = (128, 1, 1),
                            x = (0, 1000kilometers), y = (0, 1), z = (-400, 0),
                            topology = (Bounded, Periodic, Bounded))

RegularRectilinearGrid{Float64, Bounded, Periodic, Bounded}
                   domain: x ∈ [0.0, 1.0e6], y ∈ [0.0, 1.0], z ∈ [-400.0, 0.0]
                 topology: (Bounded, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (128, 1, 1)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (7812.5, 1.0, 400.0)

and Coriolis parameter appropriate for the mid-latitudes on Earth,

coriolis = FPlane(f=1e-4)

FPlane{Float64}: f = 1.00e-04

Building a `HydrostaticFreeSurfaceModel`

We use grid and coriolis to build a simple HydrostaticFreeSurfaceModel,

using Oceananigans.Models: HydrostaticFreeSurfaceModel

model = HydrostaticFreeSurfaceModel(grid=grid, coriolis=coriolis)

HydrostaticFreeSurfaceModel{CPU, Float64}(time = 0 seconds, iteration = 0) 
├── grid: RegularRectilinearGrid{Float64, Bounded, Periodic, Bounded}(Nx=128, Ny=1, Nz=1)
├── tracers: (:T, :S)
├── closure: IsotropicDiffusivity{Float64,NamedTuple{(:T, :S),Tuple{Float64,Float64}}}
├── buoyancy: SeawaterBuoyancy{Float64,LinearEquationOfState{Float64},Nothing,Nothing}
└── coriolis: FPlane{Float64}

A geostrophic adjustment initial value problem

We pose a geostrophic adjustment problem that consists of a partially-geostrophic Gaussian height field complemented by a geostrophic $y$-velocity,

Gaussian(x, L) = exp(-x^2 / 2L^2)

U = 0.1 # geostrophic velocity
L = grid.Lx / 40 # Gaussian width
x₀ = grid.Lx / 4 # Gaussian center

vᵍ(x, y, z) = - U * (x - x₀) / L * Gaussian(x - x₀, L)

g = model.free_surface.gravitational_acceleration

η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

ηᵍ(x, y, z) = η₀ * Gaussian(x - x₀, L)

ηᵍ (generic function with 1 method)

We use an initial height field that's twice the geostrophic solution, thus superimposing a geostrophic and ageostrophic component in the free surface displacement field:

ηⁱ(x, y, z) = 2 * ηᵍ(x, y, z)

ηⁱ (generic function with 1 method)

We set the initial condition to $vᵍ$ and $ηⁱ$,

set!(model, v=vᵍ, η=ηⁱ)

Running a `Simulation`

We pick a time-step that resolves the surface dynamics,

gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed

wave_propagation_time_scale = model.grid.Δx / gravity_wave_speed

simulation = Simulation(model, Δt = 0.1 * wave_propagation_time_scale, stop_iteration = 1000)

Simulation{Oceananigans.Models.HydrostaticFreeSurfaceModels.HydrostaticFreeSurfaceModel{CPU, Float64}}
├── Model clock: time = 0 seconds, iteration = 0 
├── Next time step (Float64): 12.474 seconds 
├── Iteration interval: 1
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0 seconds, wall time limit: Inf
├── Stop time: Inf years, stop iteration: 1000
├── Diagnostics: OrderedCollections.OrderedDict with 1 entry:
│   └── nan_checker => NaNChecker
└── Output writers: OrderedCollections.OrderedDict with no entries

Output

We output the velocity field and free surface displacement,

output_fields = merge(model.velocities, (η=model.free_surface.η,))

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = IterationInterval(10),
                                                      prefix = "geostrophic_adjustment",
                                                      force = true)

run!(simulation)

[ Info: Simulation is stopping. Model iteration 1000 has hit or exceeded simulation stop iteration 1000.

Visualizing the results

using JLD2, Plots, Printf, Oceananigans.Grids
using Oceananigans.Utils: hours

xη = xw = xv = xnodes(model.free_surface.η)
xu = xnodes(model.velocities.u)

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)

    u = file["timeseries/u/$iter"][:, 1, 1]
    v = file["timeseries/v/$iter"][:, 1, 1]
    η = file["timeseries/η/$iter"][:, 1, 1]
    t = file["timeseries/t/$iter"]

    titlestr = @sprintf("Geostrophic adjustment at t = %.1f hours", t / hours)

    v_plot = plot(xv / kilometers, v, linewidth = 2, title = titlestr,
                  label = "", xlabel = "x (km)", ylabel = "v (m s⁻¹)", ylims = (-U, U))

    u_plot = plot(xu / kilometers, u, linewidth = 2,
                  label = "", xlabel = "x (km)", ylabel = "u (m s⁻¹)", ylims = (-2e-3, 2e-3))

    η_plot = plot(xη / kilometers, η, linewidth = 2,
                  label = "", xlabel = "x (km)", ylabel = "η (m)", ylims = (-η₀/10, 2η₀))

    plot(v_plot, u_plot, η_plot, layout = (3, 1), size = (800, 600))
end

close(file)

This page was generated using Literate.jl.