Hi all,
I'm running a modification of the two-dimensional turbulence example on a BetaPlane, for which I need to use Bounded rather than periodic topology in the y direction. My boundary conditions are free slip in u and impenetrable in v. I find that when I set β to be quite large (~10, non-dimensional), a strong jet develops at the northern boundary, which isn't physical. I'll attach my code and a figure showing this below. Some things I've noticed:

• The jet has the opposite sign to β, but stays at the northern boundary when β changes sign
• When I change the y range from (0, 2π) to (-2π, 0), the jet moves to the southern boundary, so it's associated with non-zero Coriolis parameter.
• On an FPlane, the same jets don't appear since the dynamics are different, but when f is large there is a discontinuity in the zonal velocity at the grid cell next to the boundary. This effect seems to be linked to large Coriolis parameter in general, though it's most obvious on a BetaPlane.

Could this be related to the non-zero Neumann boundary condition for pressure when the Coriolis parameter is non-zero?

using Oceananigans
using NCDatasets
using Statistics
using CUDA

grid = RectilinearGrid(GPU(), size=(128,256), x = (0, 2*pi), y = (0, 2*pi), topology=(Periodic, Bounded, Flat))

filename = "two_dim_beta_turbulence"

coriolis = BetaPlane(f₀=0, β=10)

model = NonhydrostaticModel(; grid,
                            advection = UpwindBiased(order=5),
                            coriolis = coriolis,
                            )

# Initialise                            
u, v, w = model.velocities

uᵢ = rand(size(u)...)
vᵢ = rand(size(v)...)

uᵢ .-= mean(uᵢ)
vᵢ .-= mean(vᵢ)

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

simulation = Simulation(model, Δt=0.05, stop_time=200)

wizard = TimeStepWizard(cfl=0.7, max_change=1.1, max_Δt=0.5)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

simulation.output_writers[:fields] = NetCDFWriter(model, (; u, v),
                                                schedule = TimeInterval(50),
                                                filename = filename * ".nc",
                                                overwrite_existing = true)


run!(simulation)   

The figure shows the zonal velocity at time = 200 when (left) the y range is (-2 π, 0) and (right) (0, 2 π).

I'm using Oceananigans v0.100.7, and running on a NVIDIA RTX A5000 GPU (problem remains on CPU).

julia> versioninfo()
Julia Version 1.10.9
Commit 5595d20a287 (2025-03-10 12:51 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Core(TM) i9-14900K
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, goldmont)
Threads: 1 default, 0 interactive, 1 GC (on 32 virtual cores)