Divergent initial velocity fields cause tracer and buoyancy non-conservation with zstar

[NoraLoose]

While working on #4395 and #4403, we tested initializing simulations with random velocity fields using the zstar vertical coordinate. This resulted in tracer and buoyancy non-conservation due to the barotropic velocity field being divergent at the start. The divergence causes non-zero vertical velocities at the surface, which immediately breaks conservation.

Proposed Solutions

• Short-term: Emit a warning if the user attempts to initialize a zstar simulation with a divergent horizontal velocity field.
• Long-term: Develop a strategy to correct arbitrary initial velocity fields to ensure they are exactly divergence-free if zstar is used. This would be helpful in all cases where users don't spin up from rest, for example if they:
  • initialize with random velocities
  • initialize with velocities regridded from other simulations (which are unlikely to be perfectly divergence-free)

[glwagner]

Is it also possible that there is some inconsistency between the free surface and the velocity field divergence?

There is a similar issue with the NonhydrostaticModel which is why we subtract out the divergent part of the initial condition:

Oceananigans.jl/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl

Lines 52 to 57 in dfe0b6e

	if enforce_incompressibility
	FT = eltype(model.grid)
	calculate_pressure_correction!(model, one(FT))
	pressure_correct_velocities!(model, one(FT))
	update_state!(model; compute_tendencies = false)
	end

Possibly we need to do something similar involving a free surface spin up / fake time-step to ensure consistency there before starting a simulation, not sure.

[NoraLoose]

The problem seems to be specific to the SplitExplicitFreeSurface, see this comment.

[NoraLoose]

Great, we can probably just do what is done for the NonhyrdostaticModel.

[glwagner]

Great, we can probably just do what is done for the NonhyrdostaticModel.

Ah perhaps, but note that the divergence in the nonhydrostatic model is 3D whereas (I think) you may be referring to 2D divergence here (which is valid, but there is a vertical velocity implied). I wonder if the right solution is not to do exactly the same thing but some kind of analog, eg compute the free surface displacement consistent with the horizontal divergence / vertical velocity, or something...

[simone-silvestri]

Yeah, I think we can compute ∂t_σ such that

           δx(Δy U) + δy(Δx V)       ∇ ⋅ U
∂t_σ = - --------------------- = - --------
                  Az ⋅ H               H 

(where U and V are the barotropic velocities), and then the corresponding η associated with ∂t_σ.
Note that the continuity equation reads

  wᵏ⁺¹ - wᵏ      δx(Ax u) + δy(Ay v)     Δr ∂t_σ
 ---------- = - --------------------- - ----------
     Δz                 vol                 Δz