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"

The HydrostaticFreeSurfaceModel is still "experimental". This means some of its features are not exported, such as the ImplicitFreeSurface, and must be brought into scope manually:

using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface

## A one-dimensional domain

We use a one-dimensional domain of geophysical proportions,

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

We always have to include the z-direction for HydrostaticFreeSurfaceModel, even if the model is "barotropic" with just one grid point in z.

We deploy a 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,

model = HydrostaticFreeSurfaceModel(grid = grid,
coriolis = coriolis,
free_surface = ImplicitFreeSurface())
HydrostaticFreeSurfaceModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: RegularRectilinearGrid{Float64, Bounded, Flat, Bounded}(Nx=128, Ny=1, Nz=1)
├── tracers: (:T, :S)
├── closure: Nothing
├── buoyancy: Buoyancy{SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}, Oceananigans.Grids.ZDirection}
└── 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) = η₀ * 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) = 2 * ηᵍ(x)
ηⁱ (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.1wave_propagation_time_scale, stop_iteration = 1000)
Simulation{typename(HydrostaticFreeSurfaceModel){typename(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: typename(OrderedCollections.OrderedDict) with 1 entry:
│   └── nan_checker => typename(NaNChecker)
└── Output writers: typename(OrderedCollections.OrderedDict) with no entries

## Output

We output the velocity field and free surface displacement,

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

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
schedule = IterationInterval(10),
field_slicer = nothing,
force = true)

run!(simulation)
[ Info: Updating model auxiliary state before the first time step...
[ Info:     ... updated in 1.797 ms.
[ Info: Executing first time step...
[ Info: Simulation is stopping. Model iteration 1000 has hit or exceeded simulation stop iteration 1000.

## Visualizing the results

using Oceananigans.OutputReaders: FieldTimeSeries
using Plots, Printf

xη = xw = xv = xnodes(v_timeseries)
xu = xnodes(u_timeseries)

t = u_timeseries.times

anim = @animate for i = 1:length(t)

u = interior(u_timeseries[i])[:, 1, 1]
v = interior(v_timeseries[i])[:, 1, 1]
η = interior(η_timeseries[i])[:, 1, 1]

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

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

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

η_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