Eddy-driven jet variability using barotropic model

[KristianJS]

Milan basically implemented this StochasticStir branch at my request, in order for me to explore in more detail the ideas in Barnes and Hartmann (2011). In particular, I wanted to understand better for what parameter choices the presence of a subtropical jet can induce bimodality of the eddy-driven jet latitude.

He's asked me to create an issue describing the experiments I've done and anything I learned. The basic model I've been running is as follows:

stirlat=50
substrength=20
stirstrength=7
spectral_grid = SpectralGrid(trunc=42,nlev=1)
initial_conditions = StartFromRest()
stochastic_stirring = StochasticStirring(spectral_grid,latitude=stirlat,strength=stirstrength*1e-11)
sub_drag = JetDrag(spectral_grid,u₀=substrength)
output = OutputWriter(spectral_grid,Barotropic,output_dt=Hour(24),output_vars=[:u,:v,:vor])

model = BarotropicModel(;spectral_grid,initial_conditions,drag=sub_drag,output,forcing=stochastic_stirring)
simulation = initialize!(model)
run!(simulation_with_drag,period=Day(N),output=true)

Here the stirstrength parameter (determining the strength of the stochastic stirring) has been changed from its initial default value of 1e-11 to 7e-11. The Barnes et al. papers are not always great at clarifying parameter values, but eventually I found it clearly stated that it should be 7 in this 2010 paper. Setting it to 7 makes the eddy-driven jet have a realistic looking strength of ~10m/s, so I recommend making this the new default value.

I've then explored how the variability changes as a function of the two other parameters: substrength, which sets the strength of the subtropical jet (20m/s in Barnes and Hartmann 2011), and stirlat, which determines the latitude at which the eddy-driven jet will form. I varied each of these parameters across a wide range, and for each choice simulated ~30 years. The results are summarised in this contour plot:

The filled contour is the speed of the eddy-driven jet (stirstrength=7 always), and the black dots indicate where clear bimodality of the jet latitude is detected (essentially when two peaks of a minimum width and height can clearly be detected).

Note that here the 'eddy-driven winds' are defined by subtracting away the fixed, time-invariant subtropical jet profile from the zonal winds (which otherwise exhibit bimodality in a trivial manner, since the subtropical jet is forced to sit at the same level as the eddy-driven jet in the barotropic model). The jet speed and latitudes are computed by taking zonal averages (across the northern hemisphere) and looking for the latitude at which the maximum occurs (the speed being the maximum itself).

For context, the following figures show explicit examples of the PDFs one obtains for two selected vertical and horizontal slices:

These results show that the model is behaving in a manner consistent with Barnes and Hartmann's results, and also clarify further the role of the strength of the subtropical jet. Excellent! I'm very happy with this. I aim to next explore varying the width of the stochastic stirring, as suggested by both Tim Palmer and Tim Woollings.

[KristianJS]

In terms of supporting the repo itself, I can assist in adding comments to the code to aid readability. This is probably a good idea to help make this branch pedagogically useful.

[KristianJS]

Milan asked for the exact code. Here it is:

using SpeedyWeather, Dates

###################################################################################
#
# Definition of the 'stochastic stirring' forcing term
#
###################################################################################

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

    # DIMENSIONS from SpectralGrid
    "Spectral resolution as max degree of spherical harmonics"
    trunc::Int

    "Number of latitude rings, used for latitudinal mask"
    nlat::Int

    # OPTIONS
    "Decorrelation time scale τ [days]"
    decorrelation_time::Float64 = 2

    "Stirring strength A [1/s²]"
    strength::Float64 = 1e-11

    "Stirring latitude [˚N]"
    latitude::Float64 = 45

    "Stirring width [˚]"
    width::Float64 = 24

    "Wavenumbers which will be stirred"
    lmin::Int = 8
    lmax::Int = 12
    mmin::Int = 4

# TO BE INITIALISED
    "Stochastic stirring term S"
    S::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)

    "a = A*sqrt(1 - exp(-2dt/τ)), the noise factor times the stirring strength [1/s²]"
    a::Base.RefValue{NF} = Ref(zero(NF))

    "b = exp(-dt/τ), the auto-regressive factor [1]"
    b::Base.RefValue{NF} = Ref(zero(NF))

    "Latitudinal mask, confined to mid-latitude storm track by default [1]"
    lat_mask::Vector{NF} = zeros(NF,nlat)
end

function StochasticStirring(SG::SpectralGrid;kwargs...)
    (;trunc,Grid,nlat_half) = SG
    nlat = RingGrids.get_nlat(Grid,nlat_half)
    return StochasticStirring{SG.NF}(;trunc,nlat,kwargs...)
end

function SpeedyWeather.initialize!( forcing::StochasticStirring,
                                    model::SpeedyWeather.ModelSetup)

    # precompute forcing strength, scale with radius^2 as is the vorticity equation
    (;radius) = model.spectral_grid
    A = radius^2 * forcing.strength

    # precompute noise and auto-regressive factor, packed in RefValue for mutability
    dt = model.time_stepping.Δt_sec
    τ = forcing.decorrelation_time*24*3600   # convert to seconds
    forcing.a[] = A*sqrt(1 - exp(-2dt/τ))
    forcing.b[] = exp(-dt/τ)

    # precompute the latitudinal mask
    (;Grid,nlat_half) = model.spectral_grid
    latd = RingGrids.get_latd(Grid,nlat_half)

    for j in eachindex(forcing.lat_mask)
        # Gaussian centred at forcing.latitude of width forcing.width
        forcing.lat_mask[j] = exp(-(forcing.latitude-latd[j])^2/forcing.width^2*2)
    end

    return nothing
end

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

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

    # noise and auto-regressive factors
    a = forcing.a[]    # = sqrt(1 - exp(-2dt/τ))
    b = forcing.b[]    # = exp(-dt/τ)

    (;S) = forcing
    lmax,mmax = size(S)    # 1-based index

    for m in 1:mmax
        for l in m:lmax
            if forcing.lmin <= l <= forcing.lmax && m >= forcing.mmin
                # Barnes and Hartmann, 2011 Eq. 2
                Qi = 2rand(Complex{NF}) - (1 + im)   # ~ [-1,1] in complex
                S[l,m] = a*Qi + b*S[l,m]
            else
                S[l,m] = 0
            end
        end
    end

    # to grid-point space
    S_grid = diagn.dynamics_variables.a_grid
    SpeedyTransforms.gridded!(S_grid,S,spectral_transform)

    # mask everything but mid-latitudes
    RingGrids._scale_lat!(S_grid,forcing.lat_mask)

    # back to spectral space
    (;vor_tend) = diagn.tendencies

    # write into vor_tend by transforming S_grid to spectral
    SpeedyTransforms.spectral!(vor_tend,S_grid,spectral_transform)
    SpeedyTransforms.spectral_truncation!(vor_tend)    # set lmax+1 to zero

    return nothing
end

###################################################################################
#
# Definition of the 'subtropical jet' drag term
#
###################################################################################

Base.@kwdef struct JetDrag{NF} <: SpeedyWeather.AbstractDrag{NF}

    # DIMENSIONS from SpectralGrid
    "Spectral resolution as max degree of spherical harmonics"
    trunc::Int

    # OPTIONS
    "Relaxation time scale τ"
    time_scale::Second = Day(6)

    "Jet strength [m/s]"
    u₀::Float64 = 20

    "latitude of Gaussian jet [˚N]"
    latitude::Float64 = 30

    "Width of Gaussian jet [˚]"
    width::Float64 = 6

    # TO BE INITIALISED
    "Relaxation back to reference vorticity"
    ζ₀::LowerTriangularMatrix{Complex{NF}} = zeros(LowerTriangularMatrix{Complex{NF}},trunc+2,trunc+1)
end

function JetDrag(SG::SpectralGrid;kwargs...)
    return JetDrag{SG.NF}(;SG.trunc,kwargs...)
end

function SpeedyWeather.initialize!( drag::JetDrag,
                                    model::SpeedyWeather.ModelSetup)

    (;spectral_grid, geometry) = model
    (;Grid,NF,nlat_half) = spectral_grid
    u = zeros(Grid{NF},nlat_half)

    lat = model.geometry.latds

    for ij in RingGrids.eachgridpoint(u)
        u[ij] = drag.u₀ * exp(-(lat[ij]-drag.latitude)^2/(2*drag.width^2))
    end

    û = SpeedyTransforms.spectral(u,one_more_degree=true)
    v̂ = zero(û)
    SpeedyTransforms.curl!(drag.ζ₀,û,v̂,model.spectral_transform)
    return nothing
end

function SpeedyWeather.drag!(   diagn::SpeedyWeather.DiagnosticVariablesLayer,
                                drag::JetDrag,
                                time::DateTime,
                                model::SpeedyWeather.ModelSetup)

    (;vor_grid) = diagn.grid_variables
    (;vor_tend) = diagn.tendencies
    (;ζ₀) = drag

    vor = SpeedyTransforms.spectral(vor_grid,model.spectral_transform)

    (;radius) = model.spectral_grid
    r = radius/drag.time_scale.value
    #r = 1/drag.time_scale.value
    for lm in eachindex(vor,vor_tend,drag.ζ₀)
        vor_tend[lm] -= r*(vor[lm] - drag.ζ₀[lm])
    end

    SpeedyTransforms.spectral_truncation!(vor_tend)    # set lmax+1 to zero

    return nothing
end