2. Is there a minus sign missing in the downgradient computation? Turbulent fluxes are typically defined as minus the diffusivity multiplied by the gradient, but the code does not have this minus sign (in the plots below we flipped the sign of the diagnosed subfilter fluxes):

There's no minus signing missing in the sense that the quantity output[:νₑ_∂z_u] is equal to νₑ * ∂z(u).

But, I do agree that the name of the function that yields a dictionary of these quantities is misleading.

I think it does make sense to either 1. change the name of the function (but to what?) or 2. find a name for νₑ_∂z_u and flip the sign. For example we could have

 averages = Dict( 
                 :τ¹³ => AveragedField(- ∂z(u) * νₑ, dims=(1, 2)), 
                 :τ²³ => AveragedField(- ∂z(v) * νₑ, dims=(1, 2)), 
                 :τ³³ => AveragedField(- ∂z(w) * νₑ, dims=(1, 2)) 
                ) 

The downside of this change is that "τ¹³" is maybe less explicit than νₑ_∂z_u. But the upside is that τ¹³ is the flux rather than it's negative.

Does Oceananigans.jl consider the resolved, subfilter, and boundary flux separate?

• We obtained the resolved flux (wt) from LESbrary.TurbulenceStatistics.first_through_third_order(model). This goes to zero at the surface (see left panel below) and we thought this may be due to the subfilter flux contribution.
• We obtained the subfilter temperature flux (wt_sgs) using LESbrary.TurbulenceStatistics.subfilter_tracer_fluxes(model), but this does not capture the missing flux at the surface. It also increases the flux below the mixed layer, which is resolution-dependent (presumably because the AMD closure has a resolution dependence - shown in the right panel below).
• We noticed Oceananigans.jl code that seems to explicitly add a separate term (wt_bc) due to boundary condition fluxes (https://github.com/CliMA/Oceananigans.jl/blob/2f02bb442960439aa5fbba77e59a67829cae8833/src/BoundaryConditions/apply_flux_bcs.jl#L150-L154) so we added that term explicitly. This fixes the missing surface flux.

Oceananigans.jl calculates all interior and boundary contributions to tendency terms (for tracer and momentum) simultaneously.

But I think the important question here is about diagnostics rather than Oceananigans source code calculations. When we perform diagnostics with Oceananigans, we are always computing terms after the fact. Oceananigans does not store components of the tendency that can be evaluated later.

The subfilter tracer fluxes are defined as

On boundaries, the value of this diagnostic depends on boundary conditions. If a Flux boundary condition is used, then the flux across the boundary is determined explicitly via functions like apply_y_north_bc! as you've noted. In this case the boundary-normal tracer gradient is set to zero; therefore ∂z(c) is zero at the top when the top has a FluxBoundaryCondition.

In this case, there is no contribution to boundary fluxes from the interior turbulence closure. The way to think about this model mathematically is that the user has specified a different model for fluxes across exterior boundaries than in the interior. The only way to recover those exterior fluxes is to save them explicitly as the model is running, or to write down their values if the fluxes are constant.

If a user specifies a Value or Gradient boundary condition on the other hand, then the interior flux model is used across the boundary. To set a flux in when using Gradient boundary conditions for example, the user would also have to ensure a value for the diffusivity on the boundary (either a constant diffusivity or eddy diffusivity). If this method is used then the diagnostic ∂z(c) * κₑ should indeed return the boundary flux (provided there is only one turbulence closure present...)

Thanks for your answers, @glwagner!

We are trying to validate Oceananigans against several other LES schemes before using it for a larger project, so we wanted to make sure we were comparing apples to apples. It sounds like we were, and the distinction for Flux vs. Value/Gradient boundary conditions could also be useful for us going forward.

I see your point on the dictionary/variable names vs. the function name. If you do want to retain the dictionary terms, perhaps the functions could have alternative names like subfilter_viscous_terms and subfilter_diffusive_terms