Zanna & Bolton eddy parameterization

[milankl]

I've just implemented the Zanna & Bolton eddy parameterization

\frac{\partial \mathbf{u}}{\partial t} = ... + \kappa \nabla \cdot \mathbf{S}, \qquad \mathbf{S} =  \left( \begin{matrix} \zeta^2 - \zeta D & \zeta \bar{D} \\ \zeta \bar{D} & \zeta^2 + \zeta D \end{matrix} \right)

with $D = \partial_x v + \partial_y u$ the shear, and $\bar{D} = \partial_x u - \partial_y v$ the strain in vorticity-divergence formulation as (here now divergence $\sigma$ to not confuse it with Zanna and Bolton's notation)

\begin{aligned}
\frac{\partial \zeta}{\partial t} &= ... + \kappa \nabla \times \left( \nabla \cdot \mathbf{S} \right) \\
\frac{\partial \sigma}{\partial t} &= ... + \kappa \nabla \cdot \left( \nabla \cdot \mathbf{S} \right)
\end{aligned}

using

\begin{aligned}
\bar{D} &= 2\partial_xu - \sigma \\
D &= 2\partial_xv - \zeta
\end{aligned}

because that replaces the meridional derivative with the already available vorticity and divergence. Because $D,\bar{D}$ have to be available in grid-point space for the products with $\zeta$, we "cheat" and evaluate $\partial_xu, \partial_yu$ in grid-point space with centred 2nd order finite differences (which is easy to generalise to any ringgrid!) in order to avoid 4 transforms . In spherical coordinates with radius scaling we then have

\begin{aligned}
R\bar{D} &= \frac{2}{\cos\theta}\partial_\lambda u - R\sigma  \\
RD &= \frac{2}{\cos\theta}\partial_\lambda v - R\zeta
\end{aligned}

1. Get the tensor entries

On the grid, do (radius squared scaling omitted)

A = \frac{\zeta^2}{\cos^2\theta},\quad B = \frac{\zeta D}{\cos^2\theta},\quad C = \frac{\zeta \bar{D}}{\cos^2\theta}

Divide by cosine squared as we'll be taking two divergences/divergence+curl. Then to spectral $\hat{A},\hat{B},\hat{C}$.

2. Divergence of tensor

\begin{aligned}
\hat{G}_u &= \nabla \cdot (\hat{A} - \hat{B},\hat{C}) \\
\hat{G}_v &= \nabla \cdot (\hat{C},\hat{A} + \hat{B})
\end{aligned}

3. Divergence and curl of that for divergence and vorticity forcing

\begin{aligned}
\hat{F}_\zeta &= \nabla \times (\hat{G}_u,\hat{G}_v) \\
\hat{F}_\sigma &= \nabla \cdot (\hat{G}_u,\hat{G}_v)
\end{aligned}

Implemented as

Construction

Base.@kwdef struct ZannaBolton{NF} <: SpeedyWeather.AbstractForcing{NF}

    "Spectral resolution"
    trunc::Int

    "Strength of parameterization"
    κ::Float64 = -4.87e8

    "Work arrays: Elements on tensor"
    A::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)
    B::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)
    C::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)
    D::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)
end

ZannaBolton(SG::SpectralGrid;kwargs...) = ZannaBolton{SG.NF}(trunc=SG.trunc;kwargs...)

Initialization

function SpeedyWeather.initialize!(ZB::ZannaBolton,model::SpeedyWeather.ModelSetup)
    return nothing
end

Simulation

# function barrier
function SpeedyWeather.forcing!(diagn::SpeedyWeather.DiagnosticVariablesLayer,
                                forcing::ZannaBolton,
                                time::SpeedyWeather.DateTime,
                                model::SpeedyWeather.ModelSetup)
    forcing!(diagn,forcing,model.spectral_transform,model.geometry)
end

function forcing!(  diagn::SpeedyWeather.DiagnosticVariablesLayer,
                    forcing::ZannaBolton{NF},
                    spectral_transform::SpectralTransform,
                    geometry::SpeedyWeather.Geometry) where NF

    (;u_grid, v_grid, vor_grid, div_grid) = diagn.grid_variables
    κ = convert(NF,forcing.κ/geometry.radius^2)

    # A = ζ², half the trace of tensor S
    (;A) = forcing
    ζ² = diagn.dynamics_variables.a_grid
    for (j,ring) in enumerate(RingGrids.eachring(ζ²,vor_grid))
        κ_coslat⁻² = κ*geometry.coslat⁻²[j]
        for ij in ring
            ζ²[ij] = vor_grid[ij]^2 * κ_coslat⁻²
        end
    end
    SpeedyTransforms.spectral!(A,ζ²,spectral_transform)

    # B = ζD, the sign-changing other term on the diagonal and
    # C = ζD̄, the off diagonal term
    (;B,C) = forcing
    ζD = diagn.dynamics_variables.a_grid
    ζD̄ = diagn.dynamics_variables.b_grid
    for (j,ring) in enumerate(RingGrids.eachring(ζD,ζD̄,vor_grid,div_grid))

        n = length(ring)
        coslat⁻¹ = 2*geometry.coslat⁻¹[j]
        κ_coslat⁻² = κ*geometry.coslat⁻²[j]
        Δλ = convert(NF,2π)/n

        for (i,ij) in enumerate(ring)

            # centered finite difference gradient in zonal direction
            ij_west = ring[mod(i-1,1:n)]
            ij_east = ring[mod(i+1,1:n)]
            dudλ = (u_grid[ij_east] - u_grid[ij_west])/(2Δλ)
            dvdλ = (v_grid[ij_east] - v_grid[ij_west])/(2Δλ)

            ζD̄[ij] = κ_coslat⁻²*vor_grid[ij]*(2coslat⁻¹*dudλ - div_grid[ij])
            ζD[ij] = κ_coslat⁻²*vor_grid[ij]*(2coslat⁻¹*dvdλ - vor_grid[ij])
        end
    end
    SpeedyTransforms.spectral!(B,ζD,spectral_transform)
    SpeedyTransforms.spectral!(C,ζD̄,spectral_transform)

    # Divergence of the tensor
    AminusB = forcing.D
    AminusB .= A .- B
    Gu = diagn.dynamics_variables.a
    SpeedyTransforms.divergence!(Gu,AminusB,C,spectral_transform)

    AplusB = forcing.D
    AplusB .= A .+ B
    Gv = diagn.dynamics_variables.b
    SpeedyTransforms.divergence!(Gu,C,AplusB,spectral_transform)

    # Curl/Divergence of the ZannaBolton term
    (;vor_tend,div_tend) = diagn.tendencies
    SpeedyTransforms.curl!(vor_tend,Gu,Gv,spectral_transform)
    SpeedyTransforms.divergence!(div_tend,Gu,Gv,spectral_transform)
    return nothing
end

And then run with

forcing = ZannaBolton(spectral_grid,κ=-8e8)
model = ShallowWaterModel(;spectral_grid,forcing)

so about 2x stronger than default at T127 the blowup-test shows a self-amplifying eddy, which seems to inject energy into the flow. No more sophisticated test that the term is actually injecting energy into the system, but that sound somehow legit

@swilliamson7

[milankl]

@TomBolton Is this ☝🏼 behaviour to be expected if your parameterization is too strong? Note this is an atmospheric shallow water model setup and also I guess in general your $\kappa = -4.87e8$ strength should probably depend on the resolution too? I remember self-amplifying eddies also from backscatter experiments and it also makes dynamically sense to me if you end up putting more energy in than you take out!

[TomBolton]

No idea sorry! I haven't touched this stuff in years 😅

[milankl]

Haha, fair enough 😉

[milankl]

@swilliamson7 I'm closing this issue as I'm currently not working on this, but if you ever want to feel free to reopen!!