Surface wave effects via the Craik-Leibovich or wave-averaged Navier-Stokes formulation

[glwagner]

Eventually, we will need to represent the effects of surface waves in Oceananigans.jl.

The best way to do that is to introduce a term that models the 'wave-averaged' effects of surface waves on near-surface motions. I propose that we use the Lagrangian-mean interpretation for our velocity field in adding this term, rather than an Eulerian-mean interpretation. The Lagrangian-mean interpretation has the advantage of

1. Avoiding "accidental" initialization of large near-inertial oscillations due to initial conditions that are out of equilibrium with the surface wave field (eg the initial conditions in McWilliams et al 1997);

2. Requiring the least modification of physics, since we have to only introduce two terms (or just one for steady surface wave fields) to the momentum equations (as opposed to two terms in the momentum equations and one term in the tracer equations in the Eulerian-mean formulation);

3. Arguably using the Lagrangian-mean interpretation means we can avoid modifying our subgrid turbulence closure, since surface wave terms do not affect the turbulent kinetic energy balance;

4. The pressure field retains its original interpretation.

The surface wave field interacts with interior dynamics through its Stokes drift field, which is generally prescribed to boundary layer LES. To implement the surface wave term in the Lagrangian-mean formulation, we require gradients of the Stokes drift velocity associated with the surface wave field.

To start the best way to implement this, at least to start, is through functions. This might look something like

u_stokes(x, y, z, t) = U * exp(2*k*z)

model = Model(stokes_drift=StokesDrift(u=u_stokes), ... )

In this pattern, we have to take gradients of the function u_stokes either numerically or analytically (?)

Alternatively, we can ask users to supply the gradients of the Stokes drift of the surface wave field, though this is a bit more complicated and requires more inputs if the Stokes drift varies along coordinates other than z.

We can either call the new field model.stokes_drift or model.surface_waves.

We can validate our implementation by reproducing results in McWilliams et al (1997), Polton and Belcher (2007), and Skyllingstad and Denbo (1995).

cc @johncmarshall54

[glwagner]

Another verification experiment we can consider is

https://journals.ametsoc.org/doi/full/10.1175/1520-0485%282004%29034%3C0720%3ALESOTO%3E2.0.CO%3B2

cc @ali-ramadhan

[glwagner]

If we want to require specification of gradients of the surface wave field, we can use

uz_stokes(x, y, z, t) = 2k * U * exp(2k*z)

model = Model(surface_waves=SurfaceWavesStokesDrift(uz=uz_stokes), ...)

We should also probably use SurfaceWavesStokesDrift regardless to be explicit that the Stokes drift field is associated with surface waves.

Perhaps model.surface_waves is also preferred for this reason.

[ali-ramadhan]

@glwagner Has this been resolved with PRs #519 and #571?

[glwagner]

Yes.