Please help me understand behaviour with 0m depth wall and periodic boundaries.

[Erik-Johnson-Dal]

Hello,

I am modifying the internal tide example to attempt to make a 2d simulation of a harbour basin. I implemented a sponge layer that forces a sin function over a 1km region on the east (right hand boundary). With smoothed realistic bathymetry, the model was behaving very strangely, so I've stepped back to try to understand how bathymetry going to 0m at one edge affects the model (which currently has periodic conditions), and I'm noticing interesting behaviour that I'm not sure I understand.

Here is the code I'm running to simulate the sponge and a flat bottom:

using Oceananigans
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: PartialCellBottom

Nx, Nz = 256, 128

L = 21kilometers

H = 70
@inline bottom(x) = -H # + H*(x<1kilometer)

underlying_grid = RectilinearGrid(CPU(),size = (Nx, Nz),
                                  x = (0.0,L),
                                  z = (-H, 0),
                                  halo = (4, 4),
                                  topology = (Oceananigans.Periodic, Oceananigans.Flat, Oceananigans.Bounded))

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom))

coriolis = FPlane(latitude = 45)

T₂ = 12.421hours
ω₂ = 2π / T₂ # radians/sec
u₀ = 0.1 # tidal impulse m/s
τ_max=3hours
τ_min=1hour
sponge_width=1kilometer

@inline tau_func(x,p) = min((p.τ_max-p.τ_min)*(p.L-x)/p.sponge_width+p.τ_min,p.τ_max) # tau linear function at sponge
@inline sponge_env_func(x,p) = (tanh(10/p.sponge_width*(x-p.L+p.sponge_width))+1)/2 # sponge tanh envelope
@inline u_forcing_func(x, z, t, u, p) = sponge_env_func(x,p)*((p.u₀*sin(p.ω₂*t))-u)/tau_func(x,p)

u_forcing = Forcing(u_forcing_func,field_dependencies = :u,parameters=(;u₀ = u₀,ω₂ = ω₂, τ_max=τ_max, τ_min=τ_min, sponge_width = sponge_width, L = L))

# ## Model

# We built a `HydrostaticFreeSurfaceModel`:

model = HydrostaticFreeSurfaceModel(; grid, coriolis,
                                      buoyancy = BuoyancyTracer(),
                                      tracers = :b,
                                      momentum_advection = WENO(),
                                      tracer_advection = WENO(),
                                      forcing = (; u = u_forcing))

# We initialize the model with the zero flow and a linear stratification.

uᵢ(x, z) = 0.0

Nᵢ² = 1e-4  # [s⁻²] initial buoyancy frequency / stratification
bᵢ(x, z) = Nᵢ² * z

set!(model, u=uᵢ, b=bᵢ)

# Now let's build a `Simulation`.

Δt = 1minutes
stop_time = 4days

simulation = Simulation(model; Δt, stop_time)

# We add a callback to print a message about how the simulation is going,

using Printf

wall_clock = Ref(time_ns())

function progress(sim)
    elapsed = 1e-9 * (time_ns() - wall_clock[])

    msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹\n",
                   iteration(sim), prettytime(sim), prettytime(elapsed),
                   maximum(abs, sim.model.velocities.w))

    wall_clock[] = time_ns()

    @info msg

    return nothing
end

add_callback!(simulation, progress, name=:progress, IterationInterval(200))
nothing #hide

# ## Diagnostics/Output

# Add some diagnostics. Instead of ``u`` we save the deviation of ``u`` from its instantaneous
# domain average, ``u' = u - (L_x H)^{-1} \int u \, \mathrm{d}x \mathrm{d}z``. We also save
# the stratification ``N^2 = \partial_z b``.

b = model.tracers.b
u, v, w = model.velocities
η = model.free_surface.η
u_forcing = model.forcing.u

N² = ∂z(b)

filename = "flatBath"
save_fields_interval = 10minutes

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, w, b, η, N²);
                                                      filename,
                                                      schedule = TimeInterval(save_fields_interval),
                                                      overwrite_existing = true)

# We are ready -- let's run!

run!(simulation)

This runs through and looks quite reasonable for a flat bottom. However, if I modify the function bottom(x) to be = -H + H*(x<1kilometer) I have noticed something interesting. When the factor is equal to H (the bathymetry goes to exactly 0m in the first kilometer) there seems to be "spillover", and you can see velocity moving over the wall. However, if you change the factor to be much larger than H (say 100 meters, causing the "wall" to be 30m high), there is essentially no movement in the entire domain. Is this a case of a raised free surface allowing water to "spill" over the wall when it's 0m?

Following this, how would one go about closing one end of the simulation to model a tidal flow in a basin? Should I set the x-boundaries to bounded and give a long tail to the right of the sponge to give room for the model to relax the velocities back down to zero?

The 0m tall wall case after 2 days (note that the free surface height is plotted and scaled by 10^4):

The 30m tall wall case after 2 days (note that the free surface height is plotted and scaled by 10^4):

I can also share the actual basin bathymetry example if you think it'd help!

Cheers,
Erik

[glwagner]

However, if you change the factor to be much larger than H (say 100 meters, causing the "wall" to be 30m high), there is essentially no movement in the entire domain. Is this a case of a raised free surface allowing water to "spill" over the wall when it's 0m?

I'm surprised by this behavior, it seems wrong to me. If the bottom height is 0, then the whole column should be masked. Maybe there is some kind of round-off error, not sure, so that the free surface is not masked in the first place but is masked in the second case. It may help to check how high your "wall" needs to be. Is it identically 0 or machine eps less than 0? Does the wall height need to be machine eps greater than 0 for the free surface to be masked?

Can you check what happens if you use a Bounded topology in x?

Also, based on your description it sounds like you probably want to be using some kind of open boundary condition on the right side of the domain. Open boundary conditions are quite experimental right now, but you may be able to use @jagoosw's implementation on #3977 . Open boundary conditions can only be implemented with a Bounded topology (or put another way, Periodic topology requires Periodic boundary conditions, open boundary conditions cannot be used with a periodic topology).

[Erik-Johnson-Dal]

If you make the wall 0+epsilon, it prevents the "spill over". Similarly, if you use a GridFittedBottom instead of PartialCell, there is no flow over the closed end, as when you use a bounded x.

So, it ultimately looks like a strange behaviour limited to a very specific set of circumstances (exactly 0m at a wall, with periodic conditions in x, and a partialcellbottom).

[glwagner]

Ah, I missed that the issue is with PartialCellBottom. There may be something to fix there.