Baroclinic instability of a 2D front based on Stone (1971)

Introduction

Baroclinic instability (BCI) arises when a rotating, stratified fluid has tilted density surfaces, enabling eddies to tap available potential energy and convert it to kinetic energy. Stone (1971) (Stone, 1971) investigated non-hydrostatic effects on BCI using Eady’s framework. He found that as the Richardson number decreases, the wavelength of the most unstable mode increases while the growth rate diminishes relative to predictions from the quasigeostrophic (QG) approximation.

The basic state is given by

\[\begin{align} B(y, z) &= \text{Ri}\, z - y, \\ U(y, z) &= z - {1}/{2}, \end{align}\]

where $\text{Ri}$ is the Richardson number. We aim to analyze the stability of the above basic state against small perturbations. The perturbation variables are defined as

\[\begin{align} \mathbf{u}(x, y, z, t) &= (u, v, \epsilon w)(x, y, z, t), \\ p(x, y, z, t) &= p(x, y, z, t), \\ b(x, y, z, t) &= b(x, y, z, t), \end{align}\]

where $\epsilon$ is the aspect ratio, $\mathbf{u}$ is the velocity perturbation, $p$ is the pressure perturbation, and $b$ is the buoyancy perturbation.

Governing equations

The resulting nondimensional, linearized Boussinesq equations of motion under the $f$-plane approximation are given by

\[\begin{align} \frac{D \mathbf{u}}{Dt} + \Big(v \frac{\partial U}{\partial y} + w \frac{\partial U}{\partial z} \Big) \hat{x} + \hat{z} \times \mathbf{u} &= -\nabla p + \frac{1}{\epsilon} b \hat{z} + \text{E} \nabla^2 \mathbf{u}, \\ \frac{Db}{Dt} + v \frac{\partial B}{\partial y} + w \frac{\partial B}{\partial z} &= \frac{\text{E}}{\text{Pr}} \nabla^2 b, \\ \nabla \cdot \mathbf{u} &= 0, \end{align}\]

where $\text{E}$ is the Ekman number, $\text{Pr}$ is the Prandtl number, and

\[\begin{align} D/Dt \equiv \partial/\partial t + U (\partial/\partial x) \end{align}\]

is the material derivative. The operators:

\[\nabla \equiv (\partial/\partial x, \partial/\partial y, (1/\epsilon) \partial/\partial z),\]

\[\nabla^2 \equiv \partial^2/\partial x^2 + \partial^2/\partial y^2 + (1/\epsilon^2) \partial^2/ \partial z^2,\]

\[ \nabla_h^2 \equiv \partial^2 /\partial x^2 + \partial^2/\partial y^2.\]

To eliminate pressure, following (Teed et al., 2010), we apply the operator $\hat{z} \cdot \nabla \times \nabla \times$ and $\hat{z} \cdot \nabla \times$ to the above momentum equation. This procedure yields governing equations of three perturbation variables, the vertical velocity $w$, the vertical vorticity $\zeta \, (=\hat{z} \cdot \nabla \times \mathbf{u})$, and the buoyancy $b$

\[\begin{align} \frac{D}{Dt}\nabla^2 {w} + \frac{1}{\epsilon^2} \frac{\partial \zeta}{\partial z} &= \frac{1}{\epsilon^2} \nabla_h^2 b + \text{E} \nabla^4 w, \\ \frac{D \zeta}{Dt} - \frac{\partial U}{\partial z}\frac{\partial w}{\partial y} - \frac{\partial w}{\partial z} &= \text{E} \nabla^2 \zeta, \\ \frac{Db}{Dt} + v \frac{\partial B}{\partial y} + w \frac{\partial B}{\partial z} &= \frac{\text{E}}{\text{Pr}} \nabla^2 b, \end{align}\]

where

\[ \nabla_h^2 \equiv \partial^2 /\partial x^2 + \partial^2/\partial y^2.\]

The benefit of using the above set of equations is that it enables us to examine the instability at an along-front wavenumber $k \to 0$.

Normal mode solutions

Next we consider normal-mode perturbation solutions in the form of

\[\begin{align} [w, \zeta, b](x,y,z,t) = \mathfrak{R}\big([\tilde{w}, \tilde{\zeta}, \tilde{b}](y, z) e^{i kx + \sigma t}\big), \end{align}\]

where the symbol $\mathfrak{R}$ denotes the real part and a variable with tilde' denotes an eigenfunction. We consider that the variable\sigma` is complex with $\sigma=\sigma_r + i \sigma_i$. The real part represents the growth rate, and the imaginary part shows the frequency of the perturbation.

Finally, the following system of differential equations is obtained,

\[\begin{align} (i k U - \text{E} \mathcal{D}^2) \mathcal{D}^2 \tilde{w} + \epsilon^{-2} \partial_z \tilde{\zeta} - \epsilon^{-2} \mathcal{D}_h^2 \tilde{b} &= -\sigma \mathcal{D}^2 \tilde{w}, \\ - \partial_z U \partial_y \tilde{w} - \partial_z \tilde{w} + \left(ik U - \text{E} \mathcal{D}^2 \right) \tilde{\zeta} &= -\sigma \tilde{\zeta}, \\ \partial_z B \tilde{w} + \partial_y B \tilde{v} + \left[ik U - \text{E} \mathcal{D}^2 \right] \tilde{b} &= -\sigma \tilde{b}, \end{align}\]

where

\[\begin{align} \mathcal{D}^4 = (\mathcal{D}^2 )^2 = \big(\partial_y^2 + (1/\epsilon^2)\partial_z^2 - k^2\big)^2, \,\,\,\, \text{and} \,\, \mathcal{D}_h^2 = (\partial_y^2 - k^2). \end{align}\]

The eigenfunctions $\tilde{u}$, $\tilde{v}$ are related to $\tilde{w}$, $\tilde{\zeta}$ by the relations

\[\begin{align} -\mathcal{D}_h^2 \tilde{u} &= i k \partial_{z} \tilde{w} + \partial_y \tilde{\zeta}, \\ -\mathcal{D}_h^2 \tilde{v} &= \partial_{yz} \tilde{w} - i k \tilde{\zeta}. \end{align}\]

Boundary conditions

We choose periodic boundary conditions in the $y$-direction and free-slip, rigid lid, with zero buoyancy flux in the $z$ direction, i.e.,

\[\begin{align} \tilde{w} = \partial_{zz} \tilde{w} = \partial_z \tilde{\zeta} = \partial_z \tilde{b} = 0, \,\,\,\,\,\,\, \text{at} \,\,\, {z}=0, 1. \end{align}\]

Generalized eigenvalue problem (GEVP)

The above sets of equations with the boundary conditions can be expressed as a standard generalized eigenvalue problem (GEVP),

\[\begin{align} AX= λBX, \end{align}\]

where $\lambda$ is the eigenvalue, and $X$ is the eigenvector. The matrices $A$ and $B$ are given by

\[\begin{align} A &= \begin{bmatrix} \epsilon^2(i k U \mathcal{D}^2 -\text{E} \mathcal{D}^4) & \partial_z & -\mathcal{D}_h^2 \\ -\partial_z U \partial_y - \partial_z & i k U - \text{E} \mathcal{D}^2 & 0 \\ \partial_z B - \partial_y B (\mathcal{D}_h^1)^{-1} \partial_{yz} & ik \partial_y B (\mathcal{D}_h^2)^{-1} & ikU - \text{E} \mathcal{D}^2 \end{bmatrix}, \,\,\,\,\,\,\, B &= \begin{bmatrix} \epsilon^2 \mathcal{D}^2 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \end{align}\]

Numerical Implementation

To implement the above GEVP in a numerical code, we need to actually write following sets of equations:

\[\begin{align} A &= \begin{bmatrix} \epsilon^2(i k \operatorname{diagm}(U) \mathcal{D}^{2D} - \text{E} \mathcal{D}^{4D}) & -{D}_z^D & \mathcal{D}^{2y} \otimes I- k^2 Iₙ \\ -\operatorname{diagm}(\partial_z U) \mathcal{D}^y & i k \operatorname{diagm}(U) - \text{E} \mathcal{D}^{2N} & 0_n \\ \operatorname{diagm}(\partial_z B) - \operatorname{diagm}(\partial_y B) H \mathcal{D}^{yzD} & ik \operatorname{diagm}(\partial_y B) H & ik \operatorname{diagm}(U) - \text{E} \mathcal{D}^{2N} \end{bmatrix}, \end{align}\]

where $H$ is the inverse of the horizontal Laplacian operator $(\mathcal{D}_h^2)^{-1}$, and $\operatorname{diagm}(\phi)$ is a diagonal matrix with the elements of any vector $\phi$ on its diagonal. The right-hand side operator is discretized as

\[\begin{align} B &= \begin{bmatrix} \epsilon^2 \mathcal{D}^{2D} & 0_n & 0_n \\ 0_n & I_n & 0_n \\ 0_n & 0_n & I_n \end{bmatrix}. \end{align}\]

where $I_n$ is the identity matrix of size $(n \times n)$, where $n=N_y N_z$, $N_y$ and $N_z$ are the number of grid points in the $y$ and $z$ directions respectively. $0_n$ is the zero matrix of size $(n \times n)$. The differential operator matrices are given by

\[\begin{align} {D}^{2D} &= \mathcal{D}_y^2 \otimes {I}_z + {I}_y \otimes \mathcal{D}_z^{2D} - k^2 {I}_n, \\ {D}^{2N} &= \mathcal{D}_y^2 \otimes {I}_z + {I}_y \otimes \mathcal{D}_z^{2N} - k^2 {I}_n, \\ {D}^{4D} &= \mathcal{D}_y^4 \otimes {I}_z + {I}_y \otimes \mathcal{D}_z^{4D} + k^4 {I}_n - 2 k^2 {D}_y^2 \otimes {I}_z - 2 k^2 {I}_y \otimes {D}_z^{2D} + 2 {D}_y^2 \otimes {D}_z^{2D}, \\ {H} &= (\mathcal{D}_y^2 \otimes {I}_z - k^2 {I}_n)^{-1}, \end{align}\]

where $\otimes$ is the Kronecker product. ${I}_y$ and ${I}_z$ are identity matrices of size $(N_y \times N_y)$ and $(N_z \times N_z)$ respectively, and ${I}={I}_y \otimes {I}_z$. The superscripts $D$ and $N$ in the operator matrices denote the type of boundary conditions applied ($D$ for Dirichlet or $N$ for Neumann). $\mathcal{D}_y$, $\mathcal{D}_y^2$ and $\mathcal{D}_y^4$ are the first, second and fourth order Fourier differentiation matrix of size of $(N_y \times N_y)$, respectively. $\mathcal{D}_z$, $\mathcal{D}_z^2$ and $\mathcal{D}_z^4$ are the first, second and fourth order Chebyshev differentiation matrix of size of $(N_z \times N_z)$, respectively.

Load required packages

using LazyGrids
using LinearAlgebra
using Printf
using SparseArrays
using JLD2
using Parameters: @with_kw

using BiGSTARS
using BiGSTARS: AbstractParams
using BiGSTARS: Problem, OperatorI, TwoDGrid

Parameters

@with_kw mutable struct Params{T} <: AbstractParams
    L::T                = 1.0           # horizontal domain size
    H::T                = 1.0           # vertical domain size
    Ri::T               = 1.0           # the Richardson number
    ε::T                = 0.1           # aspect ratio ε ≡ H/L
    k::T                = 0.1           # along-front wavenumber
    E::T                = 1.0e-8        # the Ekman number
    Ny::Int64           = 24            # no. of y-grid points
    Nz::Int64           = 20            # no. of z-grid points
    w_bc::String        = "rigid_lid"   # boundary condition for vertical velocity
    ζ_bc::String        = "free_slip"   # boundary condition for vertical vorticity
    b_bc::String        = "zero_flux"   # boundary condition for buoyancy
    eig_solver::String  = "arpack"      # eigenvalue solver
end

Basic state

function basic_state(grid, params)

    Y, Z = ndgrid(grid.y, grid.z)
    Y    = transpose(Y)
    Z    = transpose(Z)

    # Define the basic state
    B   = @. 1.0 * params.Ri * Z - Y    # buoyancy
    U   = @. 1.0 * Z - 0.5 * params.H   # along-front velocity

    # Calculate all the 1st, 2nd and yz derivatives in 2D grids
    bs = compute_derivatives(U, B, grid.y; grid.Dᶻ, grid.D²ᶻ, gridtype = :All)
    precompute!(bs; which = :All)   # eager cache, returns bs itself
    @assert bs.U === U              # originals live in the same object
    @assert bs.B === B

    return bs
end

Constructing Generalized EVP

function generalized_EigValProb(prob, grid, params)

    # basic state
    bs = basic_state(grid, params)

    N  = params.Ny * params.Nz
    I⁰ = sparse(Matrix(1.0I, N, N))  # Identity matrix
    s₁ = size(I⁰, 1);
    s₂ = size(I⁰, 2);

    # the horizontal Laplacian operator:  ∇ₕ² = ∂ʸʸ - k²
    ∇ₕ² = (1.0 * prob.D²ʸ - 1.0 * params.k^2 * I⁰)

    # inverse of the horizontal Laplacian operator
    H = inverse_Lap_hor(∇ₕ²)

    # Construct the 4th order derivative
    D⁴ᴰ = (1.0 * prob.D⁴ʸ
        + 1.0/params.ε^4 * prob.D⁴ᶻᴰ
        + 1.0 * params.k^4 * I⁰
        - 2.0 * params.k^2 * prob.D²ʸ
        - 2.0/params.ε^2 * params.k^2 * prob.D²ᶻᴰ
        + 2.0/params.ε^2 * prob.D²ʸ²ᶻᴰ)

    # Construct the 2nd order derivative
    D²ᴰ = (1.0/params.ε^2 * prob.D²ᶻᴰ + 1.0 * ∇ₕ²) # with Dirichlet BC
    D²ᴺ = (1.0/params.ε^2 * prob.D²ᶻᴺ + 1.0 * ∇ₕ²) # with Neumann BC

    # See `Numerical Implementation' section for the theory
    # ──────────────────────────────────────────────────────────────────────────────
    # 1) Now define your 3×3 block-rows in a NamedTuple of 3-tuples
    # ──────────────────────────────────────────────────────────────────────────────
    # Construct the matrix `A`
    Ablocks = (
        w = (  # w-equation: [ε²(ikDiagM(U) - ED⁴ᴰ)], [Dᶻᴺ], [–∇ₕ²]
                sparse(complex.(-params.E * D⁴ᴰ + 1.0im * params.k * DiagM(bs.U) * D²ᴰ) * params.ε^2),
                sparse(complex.(prob.Dᶻᴺ)),
                sparse(complex.(-∇ₕ²))
        ),
        ζ = (  # ζ-equation: [DiagM(∂ᶻU)Dʸ - Dᶻᴰ], [kDiagM(U) – ED²ᴺ], [zero]
                sparse(complex.(-DiagM(bs.∂ᶻU) * prob.Dʸ - prob.Dᶻᴰ)),
                sparse(complex.(1.0im *params.k * DiagM(bs.U) * I⁰ - params.E * D²ᴺ)),
                spzeros(ComplexF64, s₁, s₂)
        ),
        b = (  # b-equation: [DiagM(∂ᶻB) – DiagM(∂ʸB) H Dʸᶻᴰ], [ikDiagM(∂ʸB)], [–ED²ᴺ + ikDiagM(U)]
                sparse(complex.(DiagM(bs.∂ᶻB) * I⁰ - DiagM(bs.∂ʸB) * H * prob.Dʸᶻᴰ)),
                sparse(1.0im * params.k * DiagM(bs.∂ʸB) * H),
                sparse(-params.E * D²ᴺ + 1.0im * params.k * DiagM(bs.U))
        )
    )

    # Construct the matrix `B`
    Bblocks = (
        w = (  # w-equation mass: [–ε²∂²], [zero], [zero]
                sparse(-params.ε^2 * D²ᴰ),
                spzeros(Float64, s₁, s₂),
                spzeros(Float64, s₁, s₂)
        ),
        ζ = (  # ζ-equation mass: [zero], [–I], [zero]
                spzeros(Float64, s₁, s₂),
                sparse(-I⁰),
                spzeros(Float64, s₁, s₂)
        ),
        b = (  # b-equation mass: [zero], [zero], [–I]
                spzeros(Float64, s₁, s₂),
                spzeros(Float64, s₁, s₂),
                sparse(-I⁰)
        )
    )

    # ──────────────────────────────────────────────────────────────────────────────
    # 2) Assemble the block-row matrices into a GEVPMatrices object
    # ──────────────────────────────────────────────────────────────────────────────
    gevp = GEVPMatrices(Ablocks, Bblocks)

    # ──────────────────────────────────────────────────────────────────────────────
    # 3) And now you have exactly:
    #    gevp.A, gevp.B                    → full sparse matrices
    #    gevp.As.w, gevp.As.ζ, gevp.As.b   → each block-row view of matrix A
    #    gevp.Bs.w, gevp.Bs.ζ, gevp.Bs.b   → each block-row view of matrix B
    # ──────────────────────────────────────────────────────────────────────────────

    return gevp.A, gevp.B
end

Eigenvalue solver

function EigSolver(prob, grid, params, σ₀)

    A, B = generalized_EigValProb(prob, grid, params)

    # Construct the eigenvalue solver
    # Methods available: :Krylov, :Arnoldi (by default), :Arpack
    # Here we are looking for largest growth rate (real part of eigenvalue)
    solver = EigenSolver(A, B; σ₀=σ₀, method=:Krylov, nev=1, which=:LR, sortby=:R)
    solve!(solver)
    λ, Χ = get_results(solver)
    print_summary(solver)

    # Print the largest growth rate
    @printf "largest growth rate : %1.4e%+1.4eim\n" real(λ[1]) imag(λ[1])

    return λ[1], Χ[:,1]
end

Solving the problem

function solve_Stone1971(k::Float64)

    # Calling problem parameters
    params = Params{Float64}()

    # Construct grid and derivative operators
    grid  = TwoDGrid(params)

    # Construct the necessary operator
    ops  = OperatorI(params)
    prob = Problem(grid, ops)

    # update the wavenumber
    params.k = k

    # initial guess for the growth rate
    σ₀   = 0.02

    λ, X = EigSolver(prob, grid, params, σ₀)

    # saving the result to file "stone_ms_eigenval.jld2" for the most unstable mode
    jldsave("stone_ms_eigenval.jld2";
            y=grid.y, z=grid.z, k=params.k,
            λ=λ, X=X);

    # Analytical solution of Stone (1971) for the growth rate
    cnst = 1.0 + 1.0 * params.Ri + 5.0 * params.ε^2 * params.k^2 / 42.0
    λₜ = 1.0 / (2.0 * √3.0) * (params.k - 2.0 / 15.0 * params.k^3 * cnst)

    @printf "Analytical solution of Stone (1971) for the growth rate: %f \n" λₜ

    return abs(λ.re - λₜ) < 1e-3

end

Result

solve_Stone1971(0.1) # growth rate is at k=0.1

(attempt  1/16) trying σ = 0.024000 with Krylov
  ✓ converged: λ₁ = 0.028805 + -0.000000i
(attempt  2/16) trying σ = 0.024800 with Krylov
  ✓ converged: λ₁ = 0.028805 + -0.000000i
  ✓ successive eigenvalues converged: |Δλ| = 5.57e-11 < 1.00e-05
EigenSolver Results Summary
                                                                                                                        
   Method: Krylov
   Converged: ✅ Yes
   Final shift: 0.0248
   Total time: 8.129s
   Attempts: 2
   ├─ Eigenvalues (9) ───────────────
   │ i │ λ (Re  Im·i)         │
   │───┼──────────────────────│
   │ 1 │  0.028805 -0.000000i │
   │ 2 │  0.028596 -0.000000i │
   │ 3 │  0.028569 +0.000000i │
   │ 4 │  0.027985 -0.000000i │
   │ 5 │  0.027897 +0.000000i │
   │ 6 │  0.026823 +0.000000i │
   │ 7 │  0.026754 +0.000000i │
   │ 8 │  0.025380 -0.000000i │
   │ 9 │  0.024868 -0.000000i │
   └──────────────────────────┘
                                                                                                                        
largest growth rate : 2.8805e-02-1.6004e-11im
Analytical solution of Stone (1971) for the growth rate: 0.028791 

This page was generated using Literate.jl.