Implementing Li et al's "theory wave" surface wave model

[glwagner]

We need to implement a surface wave model in our idealized LES cases. This model will output an "equilibrium" Stokes drift profile given a surface stress.

One possibility is to use a model proposed / described by Breivik et al., 2016 and Li et al. 2017. This model assumes that the surface wave spectrum is self-similar according to the Phillips spectrum to produce a Stokes drift profile. We're about to do some epic modeling so hold onto your hats --- no, really.

In the version of "theory wave" that does not account for directional spreading, the Stokes drift profile takes the form

Note that we require the vertical average of the Stokes shear for our finite volume model. This can be computed as the difference between the Stokes velocity at the top and bottom of a cell (let's take all the simple wins we can here).

This model has two parameters: the wave number at the spectral peak, k_p, and the surface Stokes drift u^S_Phil(0). Li et al propose

and

where V^S is the Stokes transport. One possibility is to use a scaling proposed by [Restrepo and McWilliams 1999] (https://journals.ametsoc.org/view/journals/phoc/29/10/1520-0485_1999_029_2523_twdoc_2.0.co_2.xml) (here T^S is Li et al's V^S):

Note that other aspects of that paper are extremely misleading so read with caution...

All of that closes one possible model for equilibrium waves.

A second possibility is to use the variant / extension of the above model proposed by Li et al. 2017 that accounts for directional spreading (ie not all waves propagate in the direction of the wind stress).

We might want to implement this in Oceananigans rather than LESbrary, but let's discuss.

[qingli411]

@glwagner It would be great to have this option implemented! It's basically a function of the wind speed. This estimate gives the depth averaged Stokes drift so you can get the layer averaged Stokes drift profile by applying this estimate to each layer depth then compute the difference. Here is the implementation in GOTM.

Another useful estimate of the Stokes drift is to use the wind-wave spectrum in Donelan et al. (1985). A modified version of it was used in the LESs of Harcourt and D'Asaro, (2008). See their Section 2b. Also see their Appendix B for computing the layer averaged Stokes drift. I used a simple version of it (without the directional spreading) in the LESs in Li and Fox-Kemper, (2017). It's also implemented in LANL's version of PALM.

I think it will be wonderful to have a Stokes drift module in either LESbrary or Oceananigans that contains different options to estimate the Stokes drift profile. Is this something you are thinking of?

[glwagner]

I think it will be wonderful to have a Stokes drift module in either LESbrary or Oceananigans that contains different options to estimate the Stokes drift profile. Is this something you are thinking of?

Yes for sure! It's probably simple enough to include directly in Oceananigans. We'd put it here:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/StokesDrift.jl

@qingli411 note that our input is the Stokes shear, so I think the most convenient thing to have is the Stokes drift at every level. We can difference that between the top and bottom of a cell to get the "finite volume averaged Stokes shear". What do you think about that?

Perhaps it makes sense to start with the Harcourt and D'Asaro 2008 model then? I'll look into that and document what I find here...

[glwagner]

The LANL-PALM reference is very useful, thanks!

[qingli411]

note that our input is the Stokes shear, so I think the most convenient thing to have is the Stokes drift at every level. We can difference that between the top and bottom of a cell to get the "finite volume averaged Stokes shear". What do you think about that?

I think it would be good to have the Stokes drift profile available somewhere. So even though only the Stokes shear is used in Oceananigans, I think it is still useful to get the Stokes drift at every level and then differentiate that to get the Stokes shear.

Perhaps it makes sense to start with the Harcourt and D'Asaro 2008 model then? I'll look into that and document what I find here...

In terms of implementation the "theory wave" is easier -- it's essentially a function of wind speed (the integration over the empirical wave spectrum is done analytically). We just need to apply it multiple time to different depths and then differentiate that to get the layer averaged Stokes drift. For the model used in Harcourt and D'Asaro (2008) we need to do the integration over the empirical wave spectrum numerically. And in addition to wind speed, it also depends on wave age. So I would probably start from the "theory wave". I guess it should be straightforward to convert the Fortran function in GOTM to Julia?

In terms of the resulting Stokes drift profile, "theory wave" usually gives a weaker Stokes shear under the same wind. The Donelan et al. (1985) spectrum used in Harcourt and D'Asaro (2008) is for pure wind-wave and gives a strong Stokes shear especially very close to the surface.

I have never used "theory wave" to drive an LES. But if we want to approximate the Stokes drift under a wide range of wave conditions with both wind waves and swell, I guess it might not be a bad option at least in an averaged sense?

[glwagner]

I think it would be good to have the Stokes drift profile available somewhere. So even though only the Stokes shear is used in Oceananigans, I think it is still useful to get the Stokes drift at every level and then differentiate that to get the Stokes shear.

I agree! Just to makes sure we're on the same page, it seems we need two separate functions to calculate both the Stokes velocity and Stokes shear averaged over a finite volume cell:

1. The finite-volume-averaged Stokes shear is the difference between a continuous expression for the Stokes velocity evaluated at the top and bottom of a cell
2. The finite-volume-averaged Stokes velocity is the difference between a continuous expression for the total-depth-averaged Stokes velocity at the top and bottom of a cell

In terms of implementation the "theory wave" is easier -- it's essentially a function of wind speed (the integration over the empirical wave spectrum is done analytically). We just need to apply it multiple time to different depths and then differentiate that to get the layer averaged Stokes drift. For the model used in Harcourt and D'Asaro (2008) we need to do the integration over the empirical wave spectrum numerically. And in addition to wind speed, it also depends on wave age. So I would probably start from the "theory wave".

Very useful information, thank you!

I guess it should be straightforward to convert the Fortran function in GOTM to Julia?

That depends strongly on my fortran-reading skills but yes, I think so! We still need to write down a continuous expression for the Stokes velocity though, I think. There's some difficulties because the Stokes velocity diverges near the surface? We can perhaps make some approximation there.

[glwagner]

@qingli411 I moved this to Oceananigans.jl because I think it's best to implement this here. It's relevant to both LES and GCM...

[glwagner]

Here's what I got:

This seems sort of reasonable at moderate wind speeds. At low wind speeds it presumably doesn't make sense to have a peak period of 2.5s ... ? At the higher wind speeds peak periods of 25.5 s are maybe not impossible.

It'd be nice to have a way of validating this. But maybe that's not easy short of actually running WW3.

Notes...

• The code to produce this is below
• Since Oceananigans doesn't yet have bulk formulae, we have to "invert" a toy bulk formula to obtain an "expected air speed" for a given kinematic momentum flux τ. I'm not sure whether it makes more sense to parameterize Stokes drift in terms of U10 or momentum flux (momentum fluxes saturate when waves stop growing at high wind speeds?)
• If we want to re-parameterize this "equilibrium Stokes drift" in terms of momentum flux, we'll presumably have to run a sweep with WW3

using SpecialFunctions
using GLMakie
using Printf

# Parameters
g = 9.81 # Gravitational acceleration for diagnostic wave property calculations
ρʷ = 1035 # Water density for toy bulk formula
ρᵃ = 1.225 # Air density for toy bulk formula
Cᵈ = 1e-3 # Drag coefficient for toy bulk formula
Ckᵖ = 0.167 # Peak wavenumber scaling parameter
Cuˢ = 0.016 # Scaling between wind stress and surface stokes drift
A_McWilliamsRestrepo = 5.1e-4 # parameter with units "inverse acceleration"
                              # that relates wind stress to Stokes transport

# Peak wave number and surface Stokes drift calculation
U₁₀(τ) = sqrt(τ * ρʷ / (Cᵈ * ρᵃ))
surface_Stokes_drift(τ) = Cuˢ * U₁₀(τ)
Stokes_transport(τ) = A_McWilliamsRestrepo * U₁₀(τ)^3

@inline function peak_wavenumber(τ)
    uˢ₀ = surface_Stokes_drift(τ)
    Vˢ = Stokes_transport(τ)
    return Ckᵖ * uˢ₀ / Vˢ
end

# Stokes drift profile calculation
T₁(k, z) = exp(2k * z)
T₂(k, z) = sqrt(2π * k * abs(z)) * erfc(sqrt(2k * abs(z)))

function Stokes_drift(τ, z)
    kᵖ = peak_wavenumber(τ)
    uˢ₀ = surface_Stokes_drift(τ)
    return uˢ₀ * (T₁(kᵖ, z) + T₂(kᵖ, z))
end

# For example...
z = -256:0.1:0

fig = Figure()
ax = Axis(fig[1, 1], xlabel="uˢ (m s⁻¹)", ylabel = "z (m)")
for τ in [1e-5, 1e-4, 2e-4, 1e-3]
    kᵖ = peak_wavenumber(τ)
    σᵖ = sqrt(g * kᵖ)
    label = @sprintf("U₁₀ = %.1f (m s⁻¹), τ = %.0e m² s⁻², λᵖ = %d m, Tᵖ = %.1f s",
                     U₁₀(τ), τ, 2π/kᵖ, 2π/σᵖ)
    uˢ = Stokes_drift.(τ, z)
    lines!(ax, uˢ, z, ; label, linewidth=3)
end

Legend(fig[0, 1], ax, tellwidth=false, tellheight=true)

display(fig)

[glwagner]

Another option is to use a model proposed by Pizzo et al. 2019 and described in the appendix of Lenain and Pizzo 2020

This is conveniently expressed in terms of the friction velocity / momentum flux rather than wind speed at 10 meters.

[qingli411]

function Stokes_drift(τ, z)
    kᵖ = peak_wavenumber(τ)
    uˢ₀ = surface_Stokes_drift(τ)
    return uˢ₀ * (T₁(kᵖ, z) + T₂(kᵖ, z))
end

I think it should be T₁(kᵖ, z) - T₂(kᵖ, z) in the above code? It might explain why the maximum Stokes drift in the figure is below the surface rather than at the surface.

It'd be nice to have a way of validating this.

One way is to use the observation at Ocean Station Papa. The UW APL group has very nice Waverider mooring observations there. Fig. 7 and Fig.8a,b of this paper is a comparison of the Stokes drift between the "theory wave" and OS Papa in one year. Here is the half-hourly Stokes drift data I computed from the wave spectrum at OS Papa used in that paper (without directional spreading). It was computed on a 1 m vertical grid in the upper 36 m and a much coarser grid below. But I can easily compute it on a finer grid if we need.

Since Oceananigans doesn't yet have bulk formulae, we have to "invert" a toy bulk formula to obtain an "expected air speed" for a given kinematic momentum flux τ. I'm not sure whether it makes more sense to parameterize Stokes drift in terms of U10 or momentum flux (momentum fluxes saturate when waves stop growing at high wind speeds?)

The reason we used U10 in the "theory wave" is that the empirical relations we used to estimate the peak wave number are based on wind speed. We need a bulk formula to convert these relations to momentum flux. The relation between U10 and surface Stokes drift does vary with U10 (see, e.g., Fig. 6 of Rascle & Ardhuin 2013). But Langmuir number also seems to vary with U10 even for fully developed waves (Below is an example showing the relation between the turbulent Langmuir number and U10 assuming Pierson-Moskowitz spectrum). So I'm not sure which way is better...

Another option is to use a model proposed by Pizzo et al. 2019 and described in the appendix of Lenain and Pizzo 2020

I haven't tried to compare with their estimate of the surface Stokes drift but it would be interesting to see!

[glwagner]

I think it should be T₁(kᵖ, z) - T₂(kᵖ, z) in the above code?

Thanks @qingli411 ! Fixed:

I talked to Nick Pizzo and spent some time reading Lenain and Melville 2017 and Lenain and Pizzo 2020 and I'm left with a few questions.

First and foremost is: how is it possible to formulate a stationary relationship between Stokes drift and transport as in McWilliams and Restrepo? In the experiments we run, we impose a constant momentum flux / constant wind on top of an initially quiescent boundary layer. Thus it seems the most appropriate model for the peak wavenumber would be one that accounts for the fetch dependence; eg, the peak wavenumber decreases in time as the wind continues to blow, consistent with the concept of a constant wind that starts blowing impulsively over an ocean at rest.

Two other details in those papers are:

• The relationships are formulated in terms of friction velocity rather than atmospheric wind speed, so we avoid introducing bulk formula which I think is an advantage.
• Lenain and Pizzo 2020 claim that their correction to the surface wave spectrum changes the surface Stokes drift (and I suppose Stokes shear) by 15-20%. I'm not sure how much this would imprint on our results but it's something to think about.

[glwagner]

Here's a few notes on a model with fetch dependence.

Time-dependence of the spectral peak

Under persistent winds, the sea state continuously evolves. This is captured by the concept of fetch, which is more or less the distance over which the wind blows. The fetch dependence of the sea state is strikingly captured by this figure from JONSWAP:

which shows the fetch dependence of the surface wave spectrum: as the fetch gets longer (as the winds blow for a longer period of time), the peak wavenumber decreases (the waves get longer) and the spectrum spreads out.

Lenain and Melville 2017 derive fetch relationships between peak frequency and non-dimensional fetch in their equations 11 and 12. Putting ourselves in a reference frame moving at the phase speed of the peak wavenumber, and rewriting these relationships in terms of peak wavenumber yields

The above is a model for a time-dependent peak wavenumber for a given constant friction velocity. This is an alternative to a constant peak wavenumber scaling.

[glwagner]

EDIT: there were a number of mistakes in the original post here. I've updated it.

As for the profiles, there are some interesting differences between the "Brevik" and "Lenain" formulations. In particular, the Lenain profile depends on the air friction velocity u★, while the Brevik does not. This is because Lenain divide the spectrum into saturation and equilibrium ranges, and the transition wavenumber between the two depends on the friction velocity (it scales with 1 / u★^2. Lenain also introduce an "isotropic wavenumber", which is an upper cutoff above which short waves do not contribute to the net Stokes drift. This isotropic wavenumber also scales with 1 / u★^2.

In all cases the Lenain model has strong shear at the surface. The effect is less pronounced on a (relatively coarse) grid with dz = 1 m, but still significant. The shear is also not monotonic in u★, possibly because the lower-frequency components start to contribute more at higher wind speeds, reducing the shear relative to the surface value of the Stokes drift.

Check out the results:

Code:

using SpecialFunctions
using GLMakie
using Printf

# Parameters
Cᵝ = 0.105  # Toba's constant
Cʳ = 9.7e-3 # Transition wavenumber parameter, Lenain and Pizzo 2020 eq 4
Cⁱ = 0.072  # Cutoff / isotropic wavenumber parameter
            # exp(π/2 - θ₀) / γ) from Lenain and Pizzo 2020 Appendix A
Cᴮ = 7e-3   # Saturation constant
 g = 9.81   # m s⁻², gravitational acceleration
ρʷ = 1024   # kg m⁻³, water density
ρᵃ = 1.225  # kg m⁻³, air density

#####
##### Stokes drift profile from Brevik et al 2016
#####

T₁(k, z) = exp(2k * z)
T₂(k, z) = sqrt(2π * k * abs(z)) * erfc(sqrt(2k * abs(z)))
brevik(k, z) = T₁(k, z) - T₂(k, z)

#####
##### Stokes drift profile from Lenain and Pizzo 2020
#####

kⁿ(u★) = Cʳ * g / u★^2 # Transition wavenumber
kⁱ(u★) = Cⁱ * g / u★^2 # Isotropic wavenumber / upper wavenumber cutoff

# Contribution from equilibrium spectrum (kᵖ < k < kⁿ)
ℓᵉ(k, z) = expinti(2k * z)
Lᵉ(kᵖ, z) = ℓᵉ(kⁿ(u★), z) - ℓᵉ(kᵖ, z)

# Contribution from saturation spectrum (kⁿ < k < kⁱ)
ℓˢ(k, z) = 2 / √(k) * (exp(2k * z) + √(2π * k * abs(z)) * erf(√(2k*abs(z))))
Lˢ(k, u★, z) = ℓˢ(kⁿ(u★), z) - ℓˢ(kⁱ(u★), z)
lenain(k, u★, z) = u★ * Cᵝ * Lᵉ(k, z) + 2Cᴮ * sqrt(g) * Lˢ(k, u★, z)

# For example...
fig = Figure()

λ = 200    # m, peak wavelength
k = 2π / λ # m⁻¹, peak wavenumber

Δz = 0.1
zᵇ = -32 + Δz/2 # center of bottom cell
zᵗ = -Δz/2 # center of top cell
zᶜ = zᵇ:Δz:zᵗ
zᶠ = @. (zᶜ[1:end-1] + zᶜ[2:end]) / 2

B = brevik.(k, zᶜ)
E = exp.(2k * zᶜ)

B ./= B[end]
E ./= E[end]

dz_B = diff(B) / Δz
dz_E = diff(E) / Δz

ax_p = Axis(fig[1, 1], xlabel="uˢ normalized to uˢ(z=0) = 1", ylabel = "z (m)")
ax_d = Axis(fig[1, 2], xlabel="∂z uˢ", ylabel = "z (m)")
ax_c = Axis(fig[1, 3], xlabel="∂z uˢ (Brevik) - ∂z uˢ (Lenain)", ylabel = "z (m)")


for τʷ = (1e-5, 1e-4, 1e-3)
    u★ = sqrt(τʷ * ρʷ / ρᵃ) # m s⁻¹

    L = lenain.(k, u★, zᶜ)
    L ./= L[end]
    dz_L = diff(L) / Δz

    suffix = @sprintf("u★ = %.2f m s⁻¹", u★)
    label = "Lenain with " * suffix
    lines!(ax_p, L, zᶜ; label, linewidth=3)
    lines!(ax_d, dz_L, zᶠ; label, linewidth=3)
    lines!(ax_c, dz_B .- dz_L, zᶠ; label=suffix, linewidth=3)
end

lines!(ax_p, B, zᶜ; label="Brevik", linewidth=3)
lines!(ax_p, E, zᶜ; label="Exponential", linewidth=3)
lines!(ax_d, dz_B, zᶠ; label="Brevik", linewidth=3)
lines!(ax_d, dz_E, zᶠ; label="Exponential", linewidth=3)

axislegend(ax_d, position=:rb)
axislegend(ax_c, position=:lb)

display(fig)

The "Lenain" profile is derived from

[glwagner]

Here's a comparison of the surface Stokes drift from the two models as a function of the kinematic (water) momentum flux:

[glwagner]

The reason we used U10 in the "theory wave" is that the empirical relations we used to estimate the peak wave number are based on wind speed. We need a bulk formula to convert these relations to momentum flux. The relation between U10 and surface Stokes drift does vary with U10 (see, e.g., Fig. 6 of Rascle & Ardhuin 2013). But Langmuir number also seems to vary with U10 even for fully developed waves (Below is an example showing the relation between the turbulent Langmuir number and U10 assuming Pierson-Moskowitz spectrum). So I'm not sure which way is better...

@qingli411 curious your thoughts on Lenain and Pizzo (2020), who provide a model in terms of friction velocity. I think their opinion is that bulk formula tend to introduce further errors into an already approximate parameterization. I suppose this model wasn't available for prior work in 2017 so perhaps the situation is different now? One unique aspect of Lenain and Pizzo (2020) is that they can co-observe the air-sea flux directly (from R/P FLIP) along with high quality measurements of wave properties.