Boundary condition not respected

[Swarnali-D]

I am setting up a simple buoyancy driven flow as:

using Oceananigans
using Oceananigans.Units: minute, minutes, hours
using Oceananigans.Solvers

grid = RectilinearGrid(GPU(),size=(50,50,150), extent=(150,150,90))
u_boundary_conditions = FieldBoundaryConditions(top = GradientBoundaryCondition(0),bottom = ValueBoundaryCondition(0))
v_boundary_conditions = FieldBoundaryConditions(top = GradientBoundaryCondition(0),bottom = ValueBoundaryCondition(0))

Qᵇ = 4.79e-9 # m² s⁻³, surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient

b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),bottom = GradientBoundaryCondition(N²))

model = NonhydrostaticModel(; grid,advection = WENO(),timestepper = :RungeKutta3,tracers = :b, buoyancy = BuoyancyTracer(), closure = ScalarDiffusivity(ν=1e-5,κ=1e-5), boundary_conditions = (u=u_boundary_conditions,v=v_boundary_conditions, b=b_boundary_conditions))

Ξ(z) = randn() * exp(z / 4)
initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z : N² * (-initial_mixed_layer_depth)

bᵢ(x, y, z) = stratification(z) + 1e-5 * Ξ(z) * N² * model.grid.Lz
set!(model, b=bᵢ)

b = model.tracers.b
∂z_b = ∂z(b)

and saving b and ∂z_b files.

During post processing, I'm getting ∂z_b values at (x,y,t) (say [25,25,:,t]) as:
151-element OffsetArray(::Vector{Float64}, 1:151) with eltype Float64 with indices 1:151: 1.9360000000000167e-5 1.9331962255287063e-5 1.9404637478112465e-5 1.938270326901103e-5 1.938021456330323e-5 1.9378401685660445e-5 1.9376282715224047e-5 1.9373811762737012e-5 1.937098753611043e-5 1.936781668513787e-5 1.9364316400258945e-5 1.936051751846351e-5 1.9356467931115113e-5 ⋮ -1.77240753562381e-7 -3.3917049937827665e-7 -4.447103124575229e-7 -4.921364411953411e-7 -3.6356250571294884e-7 -8.160312936912172e-7 -1.611274977807302e-6 -2.065157744712824e-6 -2.3156400622184894e-6 -3.453931430638144e-6 -5.738190339956502e-6 0.0

whereas I should get Qᵇ /κ = 4.79e-4 at the top instead of 0.0

[jagoosw]

Hi @Swarnali-D,

The algorithm for flux boundary conditions sets the gradient across the boundary to zero and applies the flux as a tendency term in the boundary adjacent cell, so this is the expected behaviour:
https://clima.github.io/OceananigansDocumentation/stable/numerical_implementation/boundary_conditions/#Flux-boundary-conditions

[Swarnali-D]

Then what boundary condition do I impose so that I get ∂z_b= 4.79e-4 at the top?

[glwagner]

Have you tried using GradientBoundaryCondition(4.79e-4)?

PS please use triple backticks ("```") for code blocks

[Swarnali-D]

This works! Thanks!