Boundary condition as a function of eddy viscosity

[j-ostiguy]

Hello,

I would like to implement a boundary condition which is a function of the eddy viscosity. From the documentation, I know that model.diffusivity_fields.νₑ is calculated but I don't understand how to pass it into the boundary condition function. I have the following function and model:

#DIC flux boundary condition from pCO2 flux
@inline function dic_flux_fnct(x,y,t,alk,dic,S,T,νₑ)
    ϵ = 2*νₑ*avgS
    kg = A*Sc^-n*(ϵ*ν)^0.25 #gas transfer coeff m/s 
    uc2 = 1e3*1e-6 # /L/atm to /m^3/uatm 
    k0 = 3.4e-2*uc2 # mol/m^3 /uatm
    pCO2w_val = calc_pCO2_follows(alk,dic,S,T)
    flux = kg*k0*(pCO2w_val-pCO2a) #mol/m^2s
    return flux
end

dic_flux_bc = FluxBoundaryCondition(dic_flux_fnct, field_dependencies=(:alk,:dic,:S,:T,:νₑ))
dic_bcs = FieldBoundaryConditions(top = dic_flux_bc)

model = NonhydrostaticModel(; grid, buoyancy, 
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers = (:T, :S, :alk, :dic),
                            coriolis = FPlane(f=1e-4),
                            closure = SmagorinskyLilly(C=0.16, Pr=10, Cb=1/10),
                            boundary_conditions = (u=u_bcs, T=T_bcs, S=S_bcs, alk=alk_bcs, dic=dic_bcs))

I have tried several methods but they all give me errors. Is νₑ already a field that I can just pass in or do I have to compute/define it?

[glwagner]

Unfortunately this is not super simple to achieve at the moment.

The simplest thing to do is to create a new field

auxiliary_νₑ = CenterField(grid)

and use this field in the boundary condition (either through a parameter or referencing it as a global variable). Then, add a callback to your simulation that copies the eddy diffusivity into your auxiliary field during update_state!:

function copy_eddy_diffusivity!(sim)
    model = sim.model
    diffusivity_fields = model.diffusivity_fields
    parent(auxiliary_νₑ) .= parent(diffusivity_fields.νₑ)
    return nothing
end

using Oceananigans: UpdateStateCallsite
simulation.callbacks[:copy_eddy_diffusivity] = Callback(copy_eddy_diffusivity!, callsite=UpdateStateCallsite())

It might be worthwhile to brainstorm ways to streamline this.

Finally, you will have to use the discrete_form boundary condition, because you're going to have to index into auxiliary_νₑ[i, j, grid.Nz] (if the boundary condition is at the top).

[j-ostiguy]

When I try this solution, I get the error "type NonhydrostaticModel has no field model" which I think is from the model=sim.model line. Is that correct?

[glwagner]

This issue inspired me to open an issue for a utility that might make implementing models like this a little easier: #3099

Note the solution on #3099 still requires discrete_form=true. Extending #3099 to auxiliary fields might be possible though, just a little more work (not sure it's worth it, since this is already a niche application?)

[glwagner]

What is the physical motivation for using the eddy viscosity in the boundary condition here? I have seen some work in engineering contexts that leverages the LES viscosity (as well as additional information) to compute things like the wall stress. However, I don't think such wall models are commonly used. Certainly in oceanography, simpler wall models that depend (for example) on the wind stress / friction velocity are more prevalent? Or perhaps not!

One clear issue is that any LES viscosity is a little dubious near the ocean surface --- SmagorinskyLilly more so than most. The engineering applications I mentioned above are modeling fluxes in and out of solid boundaries, not wavy free surfaces with wave breaking, etc.

I'm curious, because I want to assess how important it is to support behavior like this...

@seamanticscience might be interested in this discussion too, since it seems to involve biogeochemistry!

[j-ostiguy]

I am looking at air-sea gas exchange and parameterizations for the gas transfer velocity. The most common parameterizations are based on wind speeds, but I am attempting to implement one as a function of the turbulent kinetic energy dissipation rate.

[glwagner]

I see, that is interesting!

For SmagorinkyLilly, you can try to calculate the eddy diffusivity directly in the boundary condition:

using Oceananigans.TurbulenceClosures: ΣᵢⱼΣᵢⱼᶜᶜᶜ

@inline function top_dic_flux(i, j, grid, clock, fields, parameters)
    k = grid.Nz # top grid cell
    Δ = parameters.Δ # precalculate filter width
    C = parameters.C # free parameter
    Σ² = ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, fields.u, fields.v, fields.w)
    νₑ = C * Δ^2 * sqrt(2Σ²)
    # continue
end

top_dic_bc = FluxBoundaryCondition(top_dic_flux, discrete_form=true, parameters=dic_flux_parameters)

For your application this is a nice approach, because it will give you a bit more control over the parameterization. In particular, there are a few tricky things about computing the eddy diffusivity in the top grid cell (where we don't have velocity gradient information across the surface). So you may wish to modify the way it is done for SmagorinskyLilly specifically.

Here is the source code for SmagorinskyLilly which you can use to reference:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

[glwagner]

I was a little confused by the current implementation so I refactored it, if that helps... #3102

[j-ostiguy]

Thanks for your help!