Allow divergent Lagrangian-mean velocity by subtracting Stokes drift before pressure correction

[glwagner]

This experimental PR tweaks the pressure correction algorithm so that Lagrangian-mean velocities can be divergent. This might be convenient in practice, since otherwise users are tasked with correctly extracting the solenoidal part of the Stokes drift (eg Vanneste and Young 2022).

The simple change to the algorithm is basically to omit the dt_uS term on the RHS, and instead add the uS contribution manually. The pressure correction is performed on uE. For this to work users have to provide uS rather than dt_uS.

We could merge this, but if we do I want to change the user interface a bit. Basically I think users should either provide dt_uS or uS but not both. That allows the user to select the "mode" (divergnece-free uL, or divergence-free uE).

The new algorithm is illustrated by the following script which disproves a claim in McIntyre (1981) about the ultimate fate of oceanic momentum beneath forced surface waves (McIntyre claimed that the momentum propagates away in shallow water waves, but when the ocean is deep compared to the scale of the packet, the momentum actually stays put at the forcing site).

surface_wave_induced_flow.mp4

using Oceananigans
using GLMakie

ϵ = 0.1
λ = 60 # meters
g = 9.81

const k = 2π / λ

c = sqrt(g / k)
const δ = 500
const cᵍ = c / 2
const Uˢ = ϵ^2 * c

@inline A(ξ) = exp(- ξ^2 / 2δ^2)
@inline A′(ξ) = - ξ / δ^2 * A(ξ)
@inline A′′(ξ) = (ξ^2 / δ^2 - 1) * A(ξ) / δ^2

# Write the Stokes drift as
#
# uˢ(x, z, t) = A(x - cᵍ * t) * ûˢ(z)
#
# which describes a wave packet propogating with speed cᵍ. This implies

@inline    ûˢ(z)       = Uˢ * exp(2k * z)
@inline    uˢ(x, z, t) =         A(x - cᵍ * t) * ûˢ(z)
@inline ∂z_uˢ(x, z, t) =   2k *  A(x - cᵍ * t) * ûˢ(z)
@inline ∂t_uˢ(x, z, t) = - cᵍ * A′(x - cᵍ * t) * ûˢ(z)

# Note that if uˢ represents the solenoidal component of the Stokes drift,
# then
#
# ```math
# ∂z_wˢ = - ∂x_uˢ = - A′ * ûˢ .
# ```
#
# Thus, after integrating from bottom up to ``z`` and assuming that ``w`` at
# the bottom vanishes, we find that
#
# ```math
# wˢ = - A′ / 2k * ûˢ
# ```
#
# and

@inline ∂x_wˢ(x, z, t) = -  1 / 2k * A′′(x - cᵍ * t) * ûˢ(z)
@inline ∂t_wˢ(x, z, t) = + cᵍ / 2k * A′′(x - cᵍ * t) * ûˢ(z)
@inline    wˢ(x, z, t) = -  1 / 2k *  A′(x - cᵍ * t)  * ûˢ(z)

stokes_drift = StokesDrift(; uˢ, wˢ, ∂z_uˢ, ∂x_wˢ)
coriolis = FPlane(f=0.1)

grid = RectilinearGrid(size = (512, 1024),
                       x = (-2kilometers, 7kilometers),
                       z = (-4096, 0),
                       topology = (Periodic, Flat, Bounded))

#model = NonhydrostaticModel(; grid, stokes_drift, coriolis, tracers=:b, buoyancy=BuoyancyTracer())
model = NonhydrostaticModel(; grid, stokes_drift, tracers=:b) #, coriolis, buoyancy=BuoyancyTracer())

# Set Lagrangian-mean flow equal to uˢ,
uᵢ(x, z) = uˢ(x, z, 0)

# And put in a stable stratification,
N² = 0
bᵢ(x, z) = N² * z
set!(model, u=uᵢ, b=bᵢ)

Δx = minimum_xspacing(grid)
Δt = 0.2 * Δx / cᵍ
simulation = Simulation(model; Δt, stop_iteration = 1800)

progress(sim) = @info string("Iter: ", iteration(sim), ", time: ", prettytime(sim))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

filename = "surface_wave_induced_flow.jld2"
outputs = model.velocities
simulation.output_writers[:jld2] = JLD2OutputWriter(model, outputs; filename,
                                                    schedule = IterationInterval(10),
                                                    overwrite_existing = true)

run!(simulation)

ut = FieldTimeSeries(filename, "u")
wt = FieldTimeSeries(filename, "w")

times = ut.times
Nt = length(times)

fig = Figure(size=(800, 800))
axU = Axis(fig[1, 1], xlabel="x (m)", ylabel="U = H⁻¹ ∫ u dz (m s⁻¹)")
axu = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)")
axw = Axis(fig[3, 1], xlabel="x (m)", ylabel="z (m)")

slider = Slider(fig[4, 1], range=1:Nt, startvalue=1)
n = slider.value

rowsize!(fig.layout, 1, Relative(0.2))

u = XFaceField(grid)
U = Field(Average(u, dims=3))

Un = @lift begin
    parent(u) .= parent(ut[$n])
    compute!(U)
    U
end

un = @lift interior(ut[$n], :, 1, :)
wn = @lift interior(wt[$n], :, 1, :)

xu, yu, zu = nodes(ut)
xw, yw, zw = nodes(wt)

lines!(axU, xu, U)
heatmap!(axu, xu, zu, un, colormap=:balance, colorrange=(-5e-6, 5e-6))
heatmap!(axw, xw, zw, wn, colormap=:balance, colorrange=(-5e-6, 5e-6))

display(fig)

record(fig, "surface_wave_induced_flow.mp4", 1:Nt, framerate=12) do nn
    n[] = nn
end

[glwagner]

@BrodiePearson @qingli411 tagging you two because you may be interested in this.

I think if we want to move forward with this, possibly we should change the interface for using StokesDrift and UniformStokesDrift so that users provide the Stokes drift itself rather than the time-dependence. This also create the option of computing derivatives of the Stokes drift / components of pseudovorticity, rather than requiring users to supply them.

Thoughts?

[qingli411]

Yes, I think asking users to provide the Stokes drift itself sounds good as it allows more flexibility and we can also make sure the derivatives are consistent. One minor point I can think of is that, is there a way to do the derivatives analytically in the code then evaluate their values at appropriate locations, or we will need to evaluate Stokes drift first then do the derivatives numerically? Probably doesn't matter except very close to the surface?

[glwagner]

In principle we can use automatic differentiation to compute analytical derivatives. But I don't know if we want to bake that kind of thing into the code.

Now I am thinking, maybe its nice to allow the derivatives to be specified analytically, but we can add a feature whereby they are computed from the Stokes drift in the case they aren't provided. So then both cases are supported.

[glwagner]

Also I am not sure what is the "right" thing to do in this case. For example if we want to compute a finite volume approximation to $\partial_z u^S$ then we integrate this over a grid cell:

$$ \partial_z u^S \big |_{z_k} = \frac{1}{\Delta_z} \int_a^b \partial_z u^S \mathrm{d} z = \frac{u^S(b) - u^S(a)}{\Delta_z} $$

where $b = z_k + \Delta_z/2$ and $a = z_k - \Delta_z/2$. So perhaps to compute this more accurately, we actually do want to use this formula?

The same argument can be used for $\partial_t u^S$, since to form the discrete approximation, we integrate the momentum equation in time.

[glwagner]

Okay, my new idea is actually to stop supporting specifying derivatives because it looks like computing things directly from $u^S$ is actually more accurate! It's also clear this is much friendlier for complicated cases. If your Stokes drift varies in x, y, t then you need to compute 9 extra terms. So my back of the envelope calculation suggests that this accelerates user set up time by 900% 😆

[BrodiePearson]

I agree with this approach (specify Stokes drift and then compute relevant derivatives via finite-difference methods - rather than supplying all the derivatives).

Since the wave-averaged forces are body forces, it seems your example above often provides the exact derivative needed. There would be some error if the Stokes drift varies horizontally. Even in this case, the errors would presumably be less than the current method of specifying local derivatives to be applied within a volume-averaged equation set.

[qingli411]

So should we specify the cell-averaged Stokes drift or Stokes drift at the cell center?

[BrodiePearson]

So should we specify the cell-averaged Stokes drift or Stokes drift at the cell center?

@qingli411 I think the Stokes drift would be specified as an analytical function $u^s(x, y, z, t)$, and then the values of Stokes drift at specific locations (like the depths in Greg's example) would be computed by Oceananigans using the analytical function.

[BrodiePearson]

Is this subtracting/adding the Stokes drift at a specific depth from a depth-averaged current? Perhaps this is the kind of issue you were alluding to @qingli411

[qingli411]

Yes, I think here we should use cell-averaged Stokes drift and in the finite volume approximation of $\partial_z u^S$ we will need the values at specific locations. When Stokes drift is specified as an analytical function, we can make sure this is done consistently. But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

[glwagner]

Hmm, it's a good idea to discuss the discretization. I think on a C-grid we probably want to have the Stokes drift at the same location as the velocity components. This would not be cell centers, but Face, Center, Center for u^S, Center, Face, Center for v^S, etc.

Probably, there is a correct discretization of the term $\tilde \omega \times u^L$ that conserves energy (where $\tilde \omega = \nabla \times u^S$ is the curl of the Stokes drift), similar to how there is a correct discretization of vorticity / Coriolis that conserves energy.

For example see how we do it for Coriolis on the sphere:

Oceananigans.jl/src/Coriolis/hydrostatic_spherical_coriolis.jl

Lines 112 to 119 in d8273e7

	@inline f_ℑx_vᶠᶠᵃ(i, j, k, grid, coriolis, v) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, v)
	@inline f_ℑy_uᶠᶠᵃ(i, j, k, grid, coriolis, u) = fᶠᶠᵃ(i, j, k, grid, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, u)
	@inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
	@inbounds - ℑyᵃᶜᵃ(i, j, k, grid, f_ℑx_vᶠᶠᵃ, coriolis, U[2]) / Δxᶠᶜᶜ(i, j, k, grid)
	@inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnergyConserving, U) =
	@inbounds + ℑxᶜᵃᵃ(i, j, k, grid, f_ℑy_uᶠᶠᵃ, coriolis, U[1]) / Δyᶜᶠᶜ(i, j, k, grid)

If we want to be able to naturally compute $\omega$ at the "C-grid vorticity locations" then I believe we want u^S at the same location as the velocity field? The three components of vorticity are at cff, fcf, and ffc for the x, y, z components of vorticity respectively.

[glwagner]

I think we can build a simple example to demonstrate conservation of energy (which is $(u^L)^2 / 2$) which will show that we've done the discretization correctly.

[glwagner]

But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

That's a fair point. Though it would still be possible to compute the Stokes drift numerically by integrating its derivatives, but maybe that's not so convenient.

[BrodiePearson]

But what I'm not sure about is that, if the Stokes drift cannot be specified as an analytical function, what values should we pass in? So maybe also keep the option to pass in the derivatives of Stokes drift?

That's a fair point. Though it would still be possible to compute the Stokes drift numerically by integrating its derivatives, but maybe that's not so convenient.

The pressure solver depends only on the divergence of the velocity field. If the Stokes drift gradients were provided, would we now need to specify every derivative since the terms needed for the pressure solver Eulerian-ization ($\partial_xu^s$, $\partial_yv^s$ and $\partial_zw^s$) are distinct from the ones needed for the wave-averaged forces:

Oceananigans.jl/src/StokesDrifts.jl

Lines 153 to 167 in 0fc5708

	struct StokesDrift{P, US, VS, WS, VX, WX, UY, WY, UZ, VZ, UT, VT, WT}
	uˢ :: US
	vˢ :: VS
	wˢ :: WS
	∂x_vˢ :: VX
	∂x_wˢ :: WX
	∂y_uˢ :: UY
	∂y_wˢ :: WY
	∂z_uˢ :: UZ
	∂z_vˢ :: VZ
	∂t_uˢ :: UT
	∂t_vˢ :: VT
	∂t_wˢ :: WT
	parameters :: P
	end

[qingli411]

OK. Maybe asking the users to specify every derivative is too complicated. I agree that the easiest approach is to specify the Stokes drift as $u^S(x,y,z,t)$ and compute the required values and derivatives internally. The difference between the Stokes drift at a specific depth and a depth-averaged Stokes drift around that depth should be small if we have sufficient resolution.

[glwagner]

Maybe asking the users to specify every derivative is too complicated.

It's what we do currently 😂

We don't have to choose to do anything specific. We can keep the existing structure, and add a new one wherein users specify the Stokes drift only. Then the test of time shall arbitrate.

Mainly I think it could save us time and effort in the end to have just one implementation. Then we know that we are all using the same code and more likely to find bugs / benefit from improvements together.

The difference between the Stokes drift at a specific depth and a depth-averaged Stokes drift around that depth should be small if we have sufficient resolution.

I do think that either way, sufficient resolution renders either discretization similar. But it would be interesting to understand the implications of inadequate resolution too.

[qingli411]

I think as @BrodiePearson pointed out we will need to specify three more spatial derivatives ($\partial_x u^s$, $\partial_y v^s$, and $\partial_z w^s$)...

I agree. I was not saying we should ignore the differences. But I cannot think of a clean way to specify the Stokes drift. I think we need both the Stokes drift at a specific depth and the depth-averaged Stokes drift around that depth. Maybe asking users to specify both is another option? I was also thinking of an option to specify the wave spectrum, or a simplified one with less frequency bins as Brandon did in MOM6. But maybe these different options can be wrapped into another layer of interface in the StokesDrift module, or defined by users. Internally in Oceananigans, I guess we need $u^s$ and $v^s$ defined as depth-averaged values at the locations of $u$ and $v$, $w^s$ defined at the location of $w$, spatial derivatives should also be depth-averaged over a grid cell?

[glwagner]

We can compute the finite volume derivatives of the Stokes drift field in the same way that we compute the discrete finite volume derivatives of the prognostic velocity field. The simplest stencil is second-order accurate (it looks just like a finite difference stencil). You could be right though that the derivatives may be known with higher accuracy if users are allowed to supply them.

It wouldn't be too much work to put an interface where the needed derivatives may be optionally provided, and otherwise are computed with a finite volume method.

[BrodiePearson]

This is a good discussion, but I think I am falling into a rabbit hole. I am swinging towards @qingli411 's suggestion of the optimal solution being to specify both the Stokes drift at specific depths and the depth-averaged Stokes drift for the following reasons:

1. The Stokes drift at specific depths allows more accurate calculation of the depth-averaged Stokes drift gradients needed in the wave-vortex force term. As @glwagner demonstrated.
2. The depth-averaged Stokes drift allows accurate conversion between depth-averaged Eulerian and depth-averaged Lagrangian currents in the new pressure solver formulation. This would also come into play for the Coriolis term if Oceananigans didn't use the Lagrangian formulation (hooray!)
3. Specifying both would be the most direct way to reduce discretization errors, but at the cost of requiring the user to specify more fields.

This got me thinking more about discretization/finite-volume errors in traditional (Eulerian) LES and Oceananigans.

• I think traditional LES typically take in the Stokes drift at specific depths. This allows accurate estimation of the wave-vortex force (see above), but I suspect most also implement this local Stokes drift directly in the Stokes-Coriolis term, creating a discretization error source (should use the depth-averaged Stokes drift).
• I think the current Oceananigans implementation is precise for the Coriolis term and pressure solver, since the Lagrangian formulation never requires separation of the Stokes and Eulerian flow components. However, the specification of Stokes drift gradients at specific depths creates a discretization error source in the wave-vortex force.
• The new pressure formulation adds another source of discretization error if the supplied Stokes drift is at specific depths. If the supplied Stokes drift is the depth-averaged Stokes drift, I think it is precise. How would this method be applied if only the Stokes drift gradients are supplied?
• The above suggests that some level of discretization error is considered acceptable within the wavy-LES community, but supplying both gives an opportunity to remove another source of this error.

[glwagner]

A few preliminary thoughts: if you have a function that computes the Stokes drift, you can (if you want) also compute or pre-compute the depth-average of that field. So for Stokes drift from functions I don't know if we would ever need anything more than Stokes drift; we can compute volume averaged derivatives and volume-averaged Stokes drift from the function, if needed.

But also note that throughout the code --- for initial conditions, forcing, and boundary conditions --- we always approximate volume-averages of functions by evaluating the function in the middle of the cell. There is an issue somewhere about using more accurate quadrature but nobody has found a need for it, it would seem, so it isn't implemented...

I think we should start with something simple and build up more features / more accurate approximations as needed (eg if someone would like to test the impact of an approximation and write a paper about it, that's a good reason to make the UI more powerful and allow more aspects of the Stokes drift to be set). Note, improving the accuracy with which Stokes drift + derivatives are computed is good, but if the Stokes drift is only marginally resolved, then one might also want to worry about the accuracy of the dynamics themselves. Increasing the accuracy of the Stokes drift computation may only go so far if the equations of motion themselves are underresolved. Could be an interesting topic for research though.