How to calculate the momentum flux?

[Tinydog8]

Hi all!

I'm trying to calculate the momentum flux term (both resolved and sgs), but the result seems strange near the bounded grids.
Here is the initial data of u'w' (Center)and νₑ∂z(u) (z Face). Should I replace the bounded data with my forcing? (Such as Qᵘ)

u_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(Qᵘ), bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

Thanks!

[glwagner]

Hi! It depends on what the boundary conditions are, so can you post more of the code? For example, what is drag_bc_u and drag_bc_v?

That said we are lacking a general utility for this, and I think we need one as it would be useful to a lot of people.

EDIT: I may have been confused. Are you trying to calculate u'w' or are you trying to calculate the momentum flux at the boundary?

[glwagner]

Did you fix the problem with the staggered location?

[Tinydog8]

Yes, they are located at (Face, Center, Face) before the average operation now.

[Tinydog8]

u = model.velocities.u
v = model.velocities.v
w = model.velocities.w

νₑ = model.diffusivity_fields.νₑ
b = model.tracers.b

U = Field(Average(u, dims=(1, 2)))
V = Field(Average(v, dims=(1, 2)))
W = Field(Average(w, dims=(1, 2)))


momentum_flux_uw_op = @at (Face,Center,Face) (u-U)*(w-W)
momentum_flux_uw = Average(momentum_flux_uw_op,dims=(1, 2))
sgs_flux_uw_op = @at (Face,Center,Face) -νₑ*∂z(u)
sgs_flux_uw = Average(sgs_flux_uw_op,dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2OutputWriter(model, (; U, V,momentum_flux_uw,sgs_flux_uw),
                     schedule = AveragedTimeInterval(output_interval, window=2minutes),
                     filename = "langmuir_turbulence_averages.jld2",
                     overwrite_existing = true)


run!(simulation)

I calculate the momentum flux of the Langmuir Turbulence case in Oceananigans, the result is shown below. Do I have any incorrect operations here? (I have added the boundary condition to the top-most grid point here)

Thanks a lot!

[glwagner]

Your code looks right --- though note that W is probably zero, so you don't need to write w - W. If you do need to subtract a horizontal average from w, then your domain must be heterogeneous and W = Field(Average(w, dims=(1, 2))) wouldn't be correct anyways (you'd need to average over only part of the domain to define a decomposition of w).

However, your resolution looks coarse, perhaps 2m at the top? (You haven't included the code that makes the grid.) You probably want to be 0.5m or smaller. You might try using a grid stretched towards the top if the surface matters to you.

One subtle point is that this momentum flux calculation is an approximation if you are using a high-order advection scheme. As the resolution converges I don't think this would be an issue, but maybe it matters here.

We could try averaging the "true" momentum flux. It's a bit of work but possible.

[glwagner]

You could try this:

using Oceananigans.Advection: advective_momentum_flux_Wu

wu_op = KernelFunctionOperation{Face, Center, Face}(advective_momentum_flux_Wu, grid, model.advection, w, u)
wu = Field(Average(wu_op, dims=(1, 2))

I believe this returns $\overline{w'u'}$ --- since $\overline{w} = 0$, we have $\overline{w u} = \overline{w' (U + u')} = \overline{w' u'}$.

[glwagner]

PS I reformatted your original post to use block quotes for the code which fixed the line spacing

[ankitsb]

Hi @Tinydog8 and @glwagner,

Very interesting discussion. Any update on this?