Are Horizontal Viscosity Closures Numerically Compatible with WENO or VectorInvariant Schemes?

[kenflat2]

I found in my work that adding a horizontal viscosity closure such as horizontal_closure = HorizontalScalarDiffusivity(ν=1.2, κ=0.0) causes numerical instability in simulations that use WENO or VectorInvariant numerical schemes. In my case, including the closure makes the simulation either behave unphysically or blow up. This is counterintuitive to me because I see adding artificial viscosity as a general means to stabilize a fluid simulation.

Here is an example from one of my simulations (I can provide details about the setup if people are interested) where a horizontal closure as above is included. Here I plot snapshots of transects of zonal velocity u during my simulation

1 day in:

1 week in:

2 weeks in:

3 weeks in:

In the bottom right of the domain, a column of water becomes increasingly noisy. This noise eventually spreads and dominates the simulation. When I disable the horizontal closure, this instability does not form. Alternatively, if I make the horizontal viscosity 3 orders of magnitude larger than it is here, it damps out everything and the simulation is very well behaved (if a bit boring).

My best guess is that WENO is a diffusive scheme so adding extra diffusivity on top is not anticipated by the formulation and causes problems. Is this view correct or might something else be going on?

[glwagner]

Are you obeying the diffusive CFL constraint on the time-step? The diffusive cfl is something like kappa * dt / dx^2, where dx is the minimum horizontal spacing in the simulation.

[kenflat2]

The other quadratic bottom drag implementation I used was this:

function u_immersed_bottom_drag(i, j, k, grid, clock, fields)
    function speedᶠᶜᶜ(i, j, k, grid, fields) 
        return (fields.u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, fields.v))^0.5
    end

    μ = 1e-3

    return -μ * fields.u[i, j, k] * speedᶠᶜᶜ(i, j, k, grid, fields) 
end

function v_immersed_bottom_drag(i, j, k, grid, clock, fields)
    function speedᶜᶠᶜ(i, j, k, grid, fields) 
        return (fields.v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, fields.u))^0.5
    end

    μ = 1e-3
    return -μ * fields.v[i, j, k] * speedᶜᶠᶜ(i, j, k, grid, fields)
end

drag_u = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form=true)
drag_v = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form=true)

u_immersed_bc = ImmersedBoundaryCondition(bottom = drag_u)
v_immersed_bc = ImmersedBoundaryCondition(bottom = drag_v)

u_bcs = FieldBoundaryConditions(immersed=u_immersed_bc , top=FluxBoundaryCondition(surface_momentum_flux, discrete_form=true))
v_bcs = FieldBoundaryConditions(immersed=v_immersed_bc)

and it leads to this error

ERROR: InexactError: trunc(Int64, NaN)
Stacktrace:
  [1] trunc
    @ ./float.jl:905 [inlined]
  [2] ceil
    @ ./float.jl:384 [inlined]
  [3] calculate_substeps
    @ ~/.julia/packages/Oceananigans/u5TxZ/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl:55 [inlined]
  [4] step_free_surface!(free_surface::SplitExplicitFreeSurface{…}, model::HydrostaticFreeSurfaceModel{…}, baroclinic_timestepper::Oceananigans.TimeSteppers.QuasiAdamsBashforth2TimeStepper{…}, Δt::Float64)
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurfaces ~/.julia/packages/Oceananigans/u5TxZ/src/Models/HydrostaticFreeSurfaceModels/SplitExplicitFreeSurfaces/step_split_explicit_free_surface.jl:135
  [5] ab2_step!
    @ ~/.julia/packages/Oceananigans/u5TxZ/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl:24 [inlined]
  [6] time_step!(model::HydrostaticFreeSurfaceModel{…}, Δt::Float64; callbacks::Tuple{}, euler::Bool)
    @ Oceananigans.TimeSteppers ~/.julia/packages/Oceananigans/u5TxZ/src/TimeSteppers/quasi_adams_bashforth_2.jl:99
  [7] time_step!(sim::Simulation{…})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/u5TxZ/src/Simulations/run.jl:146
  [8] run!(sim::Simulation{…}; pickup::Bool)
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/u5TxZ/src/Simulations/run.jl:104
  [9] run!(sim::Simulation{…})
    @ Oceananigans.Simulations ~/.julia/packages/Oceananigans/u5TxZ/src/Simulations/run.jl:91
 [10] top-level scope
    @ REPL[253]:1
Some type information was truncated. Use `show(err)` to see complete types.

[kenflat2]

Here is the surface velocity field after 1 month into my latest simulation. The resolution is a quarter of a degree, I updated the CFL to 0.1, set the maximum time step to 2 minutes, added the above quadratic bottom drag, and set the momentum order to 5 using WENOVectorInvariant(order=5).

We get odd behavior in the South of the domain. Some fine grain structure/noise in meridional velocity is present in the center (which is right above the immersed boundary). Also, there are wide bands of oscillating zonal velocity at the Southern boundary. In previous simulations, instabilities such as these develop and spread to take over the simulation. This one may be stable enough to tamp out the noise, but I am not sure yet.

[glwagner]

Your surface temperature flux looks wrong:

    return Δz * ρ₀ * cₚ * (T - Tₐ(y)) / λ

a flux of any quantity has unites velocity * quantity. So this is a heat flux, not a temperature flux. You may just want Δz * (T - Tₐ(y)) / λ.

The problem is quite complicated. It might help to strip out some of the complexities. What I like to do is to set up the geometry of the problem and then start with a freely decyaying intiail condition, eg baroclinic instability. Once that seems to work well, then we can add in the forcing. The last thing to add is the drag. The resolution still looks low. Try doubling it.

[kenflat2]

Indeed, when we double the resolution the noise at the bottom becomes a rich eddy field. Here is vorticity at the surface 1 month into a 1/16 degree simulation.

I suspect these eddies are due to baroclinic instability. Below is a plot of the zonally averaged temperature in this simulation (we use a linear equation of state, so this is buoyancy). The Southern part of the domain has almost no stratification due to the initial condition and the strong temperature forcing, hence the Rossby Radius of Deformation is very small - I'm finding it to be as low as 3km.

I believe the lower resolution simulation is simply struggling to capture these detailed dynamics.

[glwagner]

For sure, there is often a critical threshold with resolution, below which you get noise and above you get turbulence.

[kenflat2]

Then yes the CFL is met, I'm using the default radius of the earth.