That's Geophysical Fluid Dynamics.

Baroclinic Instability of 3-Layer Flow

In this case, we will set up a 1000 km by 1000 km domain with three layers and study the spectral energy transfers that arise in a baroclinically unstable flow. The setup is listed below:

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 = 300*Ti      # simulation time [s]

Layer one is the thinnest, while the bottom two layers are equally thick. The layers increase in density with depth (as expected), and the zonal velocities decrease with depth. Time stepping parameters are set based on the estimate of the most unstable e-folding time scale, which is set by the Rossby radius of deformation and greatest velocity.

The model is constructed as such:

kwargs = {'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': 5000, 'tavestart': Ti*10}
m = pyqg.LayeredModel(**kwargs)

Random initial conditions are assigned to the potential vorticity, and finally the model is run:

sig = 1.e-7
q1 = np.random.randn(m.nx,m.ny)[None,]
q2 = np.random.randn(m.nx,m.ny)[None,]
q3 = np.random.randn(m.nx,m.ny)[None,]
qi = sig * np.vstack([q1, q2, q3])
m.set_q(qi)
m.run()

Here is the picture of the PV distribution in the three layers before running:

And here is after:

When setting up the parameters, we gave the deformation scale as 15 km. pyqg has a method to compute the barotropic and baroclinic deformation radii. The baroclinic radii can be computed from the eigenvalues of the “stretching” matrix:

$$\mathbf{S} \vec{p}_n = -\frac{1}{R_n^2}\vec{p}_n$$

where \(R_n\) is the nth baroclinic deformation radius. In the two-layer case:

$$
\mathbf{S} = \frac{k_d^2}{1 + \delta}\begin{pmatrix}
-1 & 1 \\ \delta & -\delta
\end{pmatrix}
$$

where \(\delta = \frac{H_1}{H_2} \). The eigenvalues are 0 and \(-(1+\delta)\), and the latter value yields the first baroclinic deformation radius \(k_d\). The barotropic radius is simply \(\frac{\sqrt{gH}}{f_0}\), which is about 1600 km.

In the three-layer case (and above), the matrix is tridiagonal:

$$
\mathbf{S} = \begin{pmatrix}
-\frac{f_0^2}{g’_1 H_1} & \frac{f_0^2}{g’_1 H_1} & 0 \\
\frac{f_0^2}{g’_1 H_2} & – (\frac{f_0^2}{g’_1 H_2} + \frac{f_0^2}{g’_2 H_2}) & \frac{f_0^2}{g’_2 H_2} \\
0 & \frac{f_0^2}{g’_2 H_3} & -\frac{f_0^2}{g’_2 H_3}
\end{pmatrix}
$$

where \(g’_k = g\frac{\rho_{k+1}-\rho_{k}}{\rho_k}\).

The eigenvectors \(\vec{p}_n\) are the three vertical modes. We can project the streamfunction m.p onto the modal bases to get the amplitudes of each mode:

pn = m.modal_projection(m.p)

Dividing these by \(U_1 L_d\) yields the normalized barotropic streamfunction and two baroclinic streamfunctions:

From here it is helpful to look at the energy spectra and the transfer of energy between different length scales. From the 2D kinetic energy spectrum provided in the diagnostics, the authors compute the radial (isotropic) kinetic and potential energy spectrum so that one can plot the quantities against a single wavenumber axis:

import pyqg.diagnostic_tools as tools
ke_bt, ke_bc1, ke_bc2 = m.get_diagnostic('KEspec_modal')
kr, modal_kespec_1 = tools.calc_ispec(m, ke_bt)
_,  modal_kespec_2 = tools.calc_ispec(m, ke_bc1)
_,  modal_kespec_3 = tools.calc_ispec(m, ke_bc2)

pe_bc1, pe_bc2 = m.get_diagnostic('PEspec_modal')
_,  modal_pespec_2 = tools.calc_ispec(m, pe_bc1)
_,  modal_pespec_3 = tools.calc_ispec(m, pe_bc2)

And lastly we can compute the fluxes of energy at different scales:

l = ['APEgenspec', 'APEflux', 'KEflux', 'KEspec']
vals = [m.get_diagnostic(x) for x in l]
APEgenspec, APEflux, KEflux = [tools.calc_ispec(m, v)[1] for v in vals[:-1]]
_, KEspec = tools.calc_ispec(m, vals[-1][1]*m.M**2)

ebud = [APEgenspec, APEflux, KEflux, -m.rek*(m.Hi[-1]/m.H)*KEspec]
resid = -np.vstack(ebud).sum(axis=0)
ebud.append(resid)
ebud_labels = ['APE gen', 'APE flux div.', 'KE flux div.', 'Diss.', 'Resid.']
[plt.semilogx(kr, term) for term in ebud]
plt.legend(ebud_labels, loc='upper right')
plt.xlim([m.kk.min(), m.kk.max()])
plt.xlabel(r'k (m$^{-1}$)'); 
plt.grid()
plt.title('Spectral Energy Transfers');

See Larichev & Held (1995) for a paper on baroclinic instability in a 2D model with similar energetic properties. APE is fluxed towards deformation length scales, where it is converted into KE. Because there is no beta effect, the KE cascades all the way up to the scale of the domain. The mechanical bottom drag removes the large scale KE.

Next we will model a fully barotropic flow, trying to reproduce the results of this James McWilliams paper.