Implementing drag boundary conditions based on similarity theory

[glwagner]

I think it's time to discuss the implementation of an abstraction for implementing drag forces on immersed boundaries. Curious to get feedback on this and also corrections if anything is wrong, plus additional considerations that I may be missing.

Background

First a bit of background on why we impose quadratic drag boundary conditions on solid surfaces. Similarity theory supposes that shear is a function of only distance $d$ from a wall (here $d$ will also be a coordinate increasing away from the wall), such that in a simple 1D situation,

$$ \partial_d u = \frac{u_\star}{\varkappa d} $$

where $u_\star$ is the friction velocity defined such that $-u_\star^2 = \tau$ where $\tau$ is the kinematic stress, $\varkappa$ is the Von Karman constant. This expression can be integrated from an "inner layer thickness" $\ell$ (often called the roughness length, although the word "roughness" may be something of a misnomer) to obtain a model for the near-wall velocity profile:

$$ u(d) = \frac{u_\star}{\varkappa} \log \frac{d}{\ell} $$

This formula may then be inverted to find the stress $\tau$ as a function of the near-wall velocity at some particular distance $d = d_0$:

$$ \tau = - u_\star^2 = - c(d_0) u(d_0)^2 $$

where $c = \varkappa^2 / log(d_0/\ell)^2$ is often called the "drag coefficient" --- again a bit misleading since it actually depends on $d_0$, ie how far we are from the wall, and so it isn't really a "coefficient" in the context of this theory (in other contexts, the drag coefficient is a non-dimensional number that characterizes the bulk drag on an object and in that case its more appropriate). But anyways. When the distance is fixed (ie when we evaluate this on a fixed mesh) --- and the roughness length is given --- then it's a constant.

The roughness length is not always constant. Above a wavy free surface, we often use the Charnock relation $\ell = c_g u_\star^2 / g$ where $g$ is gravitational acceleration and $c_g$ is the Charnock parameter, often taken to be $c_g = 0.011$. Above a smooth surface, laboratory measurements indicate $\ell = 0.11 \nu / u_\star$ (the reference for this is hard to pin down, but it dates to laboratory experiments in the 1930s, see Turbulent transfer near the interface from "Atmosphere-Ocean interaction" by Kraus and Businger).

With this background let's talk about how these parameterizations might be implemented in Oceananigans...

Constant roughness case

For an isotropic grid, constant roughness length case --- a lot of simplification! --- we have to write something like

using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ, ℑxzᶠᵃᶜ, ℑxzᶜᵃᶠ, ℑyzᵃᶠᶜ, ℑyzᵃᶜᶠ

@inline ϕ²(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k]^2
@inline speedᶠᶜᶜ(i, j, k, grid, u, v, w) = @inbounds sqrt(u[i, j, k]^2 + ℑxyᶠᶜᵃ(i, j, k, grid, ϕ², v) + ℑxzᶠᵃᶜ(i, j, k, grid, ϕ², w))
@inline speedᶜᶠᶜ(i, j, k, grid, u, v, w) = @inbounds sqrt(v[i, j, k]^2 + ℑxyᶜᶠᵃ(i, j, k, grid, ϕ², u) + ℑyzᵃᶠᶜ(i, j, k, grid, ϕ², w))
@inline speedᶜᶜᶠ(i, j, k, grid, u, v, w) = @inbounds sqrt(w[i, j, k]^2 + ℑxzᶜᵃᶠ(i, j, k, grid, ϕ², u) + ℑyzᵃᶜᶠ(i, j, k, grid, ϕ², v))
@inline u_drag(i, j, k, grid, clock, f, cᴰ) = @inbounds - cᴰ * f.u[i, j, k] * speedᶠᶜᶜ(i, j, k, grid, f.u, f.v, f.w)
@inline v_drag(i, j, k, grid, clock, f, cᴰ) = @inbounds - cᴰ * f.v[i, j, k] * speedᶜᶠᶜ(i, j, k, grid, f.u, f.v, f.w)
@inline w_drag(i, j, k, grid, clock, f, cᴰ) = @inbounds - cᴰ * f.w[i, j, k] * speedᶜᶜᶠ(i, j, k, grid, f.u, f.v, f.w)

ϰ = 0.4 # "Von Karman" constant
ℓ = 1e-4 # simply an estimated roughness length
d = minimum_xspacing(grid) / 2 # assuming an isotropic grid, using half cell-width as "distance to the wall"
cᴰ = (ϰ / log(d / ℓ))^2

u_drag_bc = FluxBoundaryCondition(u_drag, discrete_form=true, parameters=cᴰ)
v_drag_bc = FluxBoundaryCondition(v_drag, discrete_form=true, parameters=cᴰ)
w_drag_bc = FluxBoundaryCondition(w_drag, discrete_form=true, parameters=cᴰ)

Note, this only really makes sense on an isotropic grid and I think to strictly treat anisotropic grids we need to manually construct ImmersedBoundaryCondition.

Even with the simplifications we've made, it's still a lot of code that has to be repeated every time somebody wants to implemented a drag law on an immersed boundary. So we could easily motivate implementing an abstraction drag_boundary_conditions that returns the bcs for u, v, w, eg:

u_drag_bc, v_drag_bc, w_drag_bc = drag_boundary_conditions(grid; roughness_length=1e-4, von_karman_constant=0.4)

The utility could also property deal with grid stretching. And maybe omit Flat directions.

Roughness length computed from similarity theory?

In the case that we would like to use the smooth wall approximation $\ell = c_\nu \nu / u_\star$, we have to solve a transcendental equation to find the drag coefficient coefficient at every grid point (one could also formulate this as computing the flux). This could be implemented with callbacks, etc, but to formulate a nice interface for users for this we might actually have to add some kind of update_boundary_conditions! feature to update_state! that could precompute the drag coefficient and/or fluxes.

In terms of a path forward, I think we could simply start with the constant roughness length case, and perhaps make the utility a bit general so that users could also directly specify a drag coefficient if desired. For example in hydrostatic cases we often specify the drag coefficient directly, and we also usually omit the dependence on w.

There are even more considerations one might consider... for example, in a finite volume model (and in the code above) the distance to the wall is taken as the "center of the cell", ie half the cell thickness. But this is not really consistent with the finite volume framework, and better approaches have been proposed (eg see Nishizawa and Kitamura 2018).

Note that if we want to compute fluxes across the air-sea interface, it's probably better to use ClimaOcean (tools for building coupled and realistic air-sea, air-ice-sea models) for that. But ClimaOcean is not going to support the algorithmically simpler, yet geometrically more complicated case of computing momentum fluxes into complex solid objects. I think this is in scope for Oceananigans since it is not a concept in coupling between two fluids, like the air-sea case handled by ClimaOcean is.

[jagoosw]

This sounds really useful!

I kept reusing some code for wall (just with a flat-bottomed domain) and wind stress so I put some bits together in this package and I came to the same conclusion of having functions that make the boundary conditions so I could just write:

stress_boundary_conditions = WallStressBoundaryConditions()

boundary_conditions = (u = FieldBoundaryConditions(bottom = stress_boundary_conditions.u),
                       v = FieldBoundaryConditions(bottom = stress_boundary_conditions.v))

The way I did it was to make a type for WindStress which could have different models for the drag coefficient and parameters. The only other thing I learned from doing it was that it was very useful to precompute roughness lengths for different speeds rather than compute them on the fly since they're well behaved functions that I could interpolate easily.

(There is some inefficiency in the code that I did the wind and wall stress separately, even though they're basically the same thing, but a kind of merge of how I did each would probably be the most general way todo it because I put more work into the wind stress and just have the very simple flat bottomed law of the wall wall stress. Also, the code is not well-tested and probably won't work on non-rectilinear grids).

[glwagner]

Huh, I wasn't aware of that work! I'm amazed at how much is in there. Some of the air-sea stuff is overlapping with work we've done at ClimaOcean.

I had assumed that we would want to implement bulk formula in the context of ClimaOcean's OceanSeaIceModel. But that package is taking a different approach, so maybe we should consider that. How should we proceed? I do think it's best to work on a single tool rather than copy or replicate efforts in many places.

As for wall stress models for rough boundaries, given what I see in that package I feel it makes sense to design some generalizable infrastructure in Oceananigans for implementing / testing different wall models. What I've written above is merely one choice I think.

[glwagner]

(There is some inefficiency in the code that I did the wind and wall stress separately, even though they're basically the same thing, but a kind of merge of how I did each would probably be the most general way todo it because I put more work into the wind stress and just have the very simple flat bottomed law of the wall wall stress.

Do you mean computational inefficiency, or inefficiency in terms of code length? For flat surfaces, we've designed an interface in ClimaOcean that depends on a roughness length callable object. This allows one to implement new models for the roughness length (which limit to wavy surface or smooth walls, combinations of the two, or other models). More generally our intent is definitely to encapsulate the parts that are not worth repeating elsewhere, and let users / external packages proliferate to do the interesting parts like designing new roughness lengths. Was that your vision too?