Fully developed baroclinic instability of a 3-layer flow

import numpy as np
from numpy import pi
from matplotlib import pyplot as plt

import pyqg
from pyqg import diagnostic_tools as tools

Set up

L =  1000.e3     # length scale of box    [m]
Ld = 15.e3       # deformation scale      [m]
kd = 1./Ld       # deformation wavenumber [m^-1]
Nx = 64          # number of grid points

H1 = 500.        # layer 1 thickness  [m]
H2 = 1750.       # layer 2
H3 = 1750.       # layer 3

U1 = 0.05          # layer 1 zonal velocity [m/s]
U2 = 0.025         # layer 2
U3 = 0.00          # layer 3

rho1 = 1025.
rho2 = 1025.275
rho3 = 1025.640

rek = 1.e-7       # linear bottom drag coeff.  [s^-1]
f0  = 0.0001236812857687059 # coriolis param [s^-1]
beta = 1.2130692965249345e-11 # planetary vorticity gradient [m^-1 s^-1]

Ti = Ld/(abs(U1))  # estimate of most unstable e-folding time scale [s]
dt = Ti/200.   # time-step [s]
tmax = 500*Ti      # simulation time [s]

m = pyqg.LayeredModel(nx=Nx, nz=3, U = [U1,U2,U3],V = [0.,0.,0.],L=L,f=f0,beta=beta,
                         H = [H1,H2,H3], rho=[rho1,rho2,rho3],rek=rek,
                        dt=dt,tmax=tmax, twrite=10000, tavestart=Ti*200)

INFO:  Logger initialized

Initial condition

sig = 1.e-7
qi = sig*np.vstack([np.random.randn(m.nx,m.ny)[np.newaxis,],
                    np.random.randn(m.nx,m.ny)[np.newaxis,],
                    np.random.randn(m.nx,m.ny)[np.newaxis,]])
m.set_q(qi)

Run the model

m.run()

INFO: Step: 10000, Time: 1.50e+07, KE: 2.38e-04, CFL: 0.009
INFO: Step: 20000, Time: 3.00e+07, KE: 3.59e-02, CFL: 0.126
INFO: Step: 30000, Time: 4.50e+07, KE: 1.70e-01, CFL: 0.191
INFO: Step: 40000, Time: 6.00e+07, KE: 3.59e-01, CFL: 0.213
INFO: Step: 50000, Time: 7.50e+07, KE: 1.91e-01, CFL: 0.192
INFO: Step: 60000, Time: 9.00e+07, KE: 2.24e-01, CFL: 0.207
INFO: Step: 70000, Time: 1.05e+08, KE: 5.46e-01, CFL: 0.307
INFO: Step: 80000, Time: 1.20e+08, KE: 4.94e-01, CFL: 0.252
INFO: Step: 90000, Time: 1.35e+08, KE: 5.62e-01, CFL: 0.210
INFO: Step: 100000, Time: 1.50e+08, KE: 2.70e-01, CFL: 0.218

Xarray Dataset

Notice that the conversion to an xarray dataset requires xarray to be installed on your machine. See here for installation instructions: http://xarray.pydata.org/en/stable/getting-started-guide/installing.html#instructions

ds = m.to_dataset()
ds

<xarray.Dataset>
Dimensions:            (time: 1, lev: 3, y: 64, x: 64, l: 64, k: 33, lev_mid: 2)
Coordinates:
  * time               (time) float64 1.5e+08
  * lev                (lev) int64 1 2 3
  * lev_mid            (lev_mid) float64 1.5 2.5
  * x                  (x) float64 7.812e+03 2.344e+04 ... 9.766e+05 9.922e+05
  * y                  (y) float64 7.812e+03 2.344e+04 ... 9.766e+05 9.922e+05
  * l                  (l) float64 0.0 6.283e-06 ... -1.257e-05 -6.283e-06
  * k                  (k) float64 0.0 6.283e-06 ... 0.0001948 0.0002011
Data variables: (12/34)
    q                  (time, lev, y, x) float64 0.0001383 ... -8.462e-06
    u                  (time, lev, y, x) float64 0.0413 -0.05882 ... -0.1616
    v                  (time, lev, y, x) float64 -0.357 -0.3027 ... -0.1369
    ufull              (time, lev, y, x) float64 0.0913 -0.008824 ... -0.1616
    vfull              (time, lev, y, x) float64 -0.357 -0.3027 ... -0.1369
    qh                 (time, lev, l, k) complex128 (9.956089311189431e-06+0j...
    ...                 ...
    APEgenspec         (time, l, k) float64 0.0 1.724e-08 ... 5.211e-50
    KEspec_modal       (time, lev, l, k) float64 0.0 0.1616 ... 2.567e-41
    PEspec_modal       (time, lev_mid, l, k) float64 0.0 0.02742 ... 9.974e-42
    APEspec            (time, l, k) float64 0.0 0.04085 ... 2.995e-32 1.117e-41
    KEflux_div         (time, l, k) float64 0.0 6.588e-09 ... 1.817e-24 3.27e-29
    APEflux_div        (time, l, k) float64 0.0 -1.748e-08 ... 1.176e-29
Attributes: (12/24)
    pyqg:beta:       1.2130692965249345e-11
    pyqg:delta:      None
    pyqg:dt:         1500.0
    pyqg:filterfac:  23.6
    pyqg:L:          1000000.0
    pyqg:M:          4096
    ...              ...
    pyqg:tc:         100000
    pyqg:tmax:       150000000.0
    pyqg:twrite:     10000
    pyqg:W:          1000000.0
    title:           pyqg: Python Quasigeostrophic Model
    reference:       https://pyqg.readthedocs.io/en/latest/index.html

Snapshots

PV = ds.q + ds.Qy * ds.y
PV['x'] = ds.x/ds.attrs['pyqg:rd']; PV.x.attrs = {'long_name': r'$x/L_d$'}
PV['y'] = ds.y/ds.attrs['pyqg:rd']; PV.y.attrs = {'long_name': r'$y/L_d$'}

plt.figure(figsize=(18,4))

plt.subplot(131)
PV.sel(lev=1).plot(cmap='Spectral_r')

plt.subplot(132)
PV.sel(lev=2).plot(cmap='Spectral_r')

plt.subplot(133)
PV.sel(lev=3).plot(cmap='Spectral_r');

pyqg has a built-in method that computes the vertical modes. It is stored as an attribute in the Dataset.

print(f"The first baroclinic deformation radius is {ds.attrs['pyqg:radii'][1]/1.e3} km")
print(f"The second baroclinic deformation radius is {ds.attrs['pyqg:radii'][2]/1.e3} km")

The first baroclinic deformation radius is 15.375382785987185 km
The second baroclinic deformation radius is 7.975516271996243 km

We can project the solution onto the modes

pn = m.modal_projection(m.p)

plt.figure(figsize=(18,4))

plt.subplot(131)
plt.pcolormesh(m.x/m.rd, m.y/m.rd, pn[0]/(U1*Ld), cmap='Spectral_r', shading='auto')
plt.xlabel(r'$x/L_d$')
plt.ylabel(r'$y/L_d$')
plt.colorbar()
plt.title('Barotropic streamfunction')

plt.subplot(132)
plt.pcolormesh(m.x/m.rd, m.y/m.rd, pn[1]/(U1*Ld), cmap='Spectral_r', shading='auto')
plt.xlabel(r'$x/L_d$')
plt.ylabel(r'$y/L_d$')
plt.colorbar()
plt.title('1st baroclinic streamfunction')

plt.subplot(133)
plt.pcolormesh(m.x/m.rd, m.y/m.rd, pn[2]/(U1*Ld), cmap='Spectral_r', shading='auto')
plt.xlabel(r'$x/L_d$')
plt.ylabel(r'$y/L_d$')
plt.colorbar()
plt.title('2nd baroclinic streamfunction');

Diagnostics

kr, kespec_1 = tools.calc_ispec(m, ds.KEspec.sel(lev=1).data.squeeze())
_ , kespec_2 = tools.calc_ispec(m, ds.KEspec.sel(lev=2).data.squeeze())
_ , kespec_3 = tools.calc_ispec(m, ds.KEspec.sel(lev=3).data.squeeze())

plt.loglog(kr, kespec_1, '.-' )
plt.loglog(kr, kespec_2, '.-' )
plt.loglog(kr, kespec_3, '.-' )

plt.legend(['layer 1','layer 2', 'layer 3'], loc='lower left')
plt.ylim([1e-12,4e-6]);
plt.xlabel(r'k (m$^{-1}$)'); plt.grid()
plt.title('Kinetic Energy Spectrum');

By default the modal KE and PE spectra are also calculated

kr, modal_kespec_1 = tools.calc_ispec(m, ds.KEspec_modal.sel(lev=1).data.squeeze())
_,  modal_kespec_2 = tools.calc_ispec(m, ds.KEspec_modal.sel(lev=2).data.squeeze())
_,  modal_kespec_3 = tools.calc_ispec(m, ds.KEspec_modal.sel(lev=3).data.squeeze())

_,  modal_pespec_2 = tools.calc_ispec(m, ds.PEspec_modal.sel(lev_mid=1.5).data.squeeze())
_,  modal_pespec_3 = tools.calc_ispec(m, ds.PEspec_modal.sel(lev_mid=2.5).data.squeeze())

plt.figure(figsize=(15,5))

plt.subplot(121)
plt.loglog(kr, modal_kespec_1, '.-')
plt.loglog(kr, modal_kespec_2, '.-')
plt.loglog(kr, modal_kespec_3, '.-')

plt.legend(['barotropic ','1st baroclinic', '2nd baroclinic'], loc='lower left')
plt.ylim([1e-12,1e-5]);
plt.xlabel(r'k (m$^{-1}$)'); plt.grid()
plt.title('Kinetic Energy Spectra');


plt.subplot(122)
plt.loglog(kr, modal_pespec_2, '.-')
plt.loglog(kr, modal_pespec_3, '.-')

plt.legend(['1st baroclinic', '2nd baroclinic'], loc='lower left')
plt.ylim([1e-12,1e-5]);
plt.xlabel(r'k (m$^{-1}$)'); plt.grid()
plt.title('Potential Energy Spectra');

_, APEgenspec = tools.calc_ispec(m, ds.APEgenspec.data.squeeze())
_, APEflux    = tools.calc_ispec(m, ds.APEflux_div.data.squeeze())
_, KEflux     = tools.calc_ispec(m, ds.KEflux_div.data.squeeze())
_, Dissspec   = tools.calc_ispec(m, ds.Dissspec.squeeze().data)
_, KEfrictionspec = tools.calc_ispec(m, ds.KEfrictionspec.squeeze().data)

ebud = [ APEgenspec,
         APEflux,
         KEflux,
         Dissspec,
         KEfrictionspec]
ebud.append(-np.vstack(ebud).sum(axis=0))
ebud_labels = ['APE gen','APE flux div.','KE flux div.','Dissipation','Friction','Resid.']
[plt.semilogx(kr, term) for term in ebud]
plt.legend(ebud_labels, loc='lower left', ncol=2)
plt.xlabel(r'k (m$^{-1}$)'); plt.grid()
plt.title('Spectral Energy Transfers');

_, ENSflux    = tools.calc_ispec(m, ds.ENSflux.data.squeeze())
_, ENSgenspec = tools.calc_ispec(m, ds.ENSgenspec.data.squeeze())
_, ENSfrictionspec = tools.calc_ispec(m, ds.ENSfrictionspec.data.squeeze())
_, ENSDissspec = tools.calc_ispec(m, ds.ENSDissspec.data.squeeze())

ebud = [ ENSgenspec,
         ENSflux,
         ENSDissspec,
         ENSfrictionspec]
ebud.append(-np.vstack(ebud).sum(axis=0))
ebud_labels = ['ENS gen','ENS flux div.','Dissipation','Friction','Resid.']
[plt.semilogx(kr, term) for term in ebud]
plt.legend(ebud_labels, loc='lower left')
plt.xlabel(r'k (m$^{-1}$)'); plt.grid()
plt.title('Spectral Enstrophy Transfer');

The dynamics here is similar to the reference experiment of Larichev & Held (1995). The APE generated through baroclinic instability is fluxed towards deformation length scales, where it is converted into KE. The KE the experiments and inverse tranfer, cascading up to the scale of the domain. The mechanical bottom drag essentially removes the large scale KE.