This notebook is studying the vortices shed by the WCR. Inspired by output from runew-06-spall-const-5d-4-shc

Hypotheses are:

  1. Barotropic instability of near shelfbreak jet - solve the barotropic instability problem for a Gaussian jet.
  2. WCR exports vorticity in bursts
  3. Vortex rollup of sheet exported by WCR

In [1]:
%matplotlib inline

from numpy import *
import h5py
import netCDF4
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['savefig.dpi'] = 2 * matplotlib.rcParams['savefig.dpi']
import numexpr as ne
from import loadmat

import pyroms not found. Remapping function will not be available

First, bring up plots

In [2]:
dirname = 'runs/topoeddy/runew-06-spall-const-5d-4-shc/'
fname = dirname +  ''
iname = dirname + '/config/'
edname = dirname + 'eddytrack.mat'

# This time index looks good
tind = 20

# fftconvolve is pretty slow...
def avg1(variable, axis):
    import scipy
    from scipy import signal
    cnv = scipy.ones(variable.ndim)
    cnv[axis] = 2
    cnv = scipy.ones(cnv)
    out = signal.fftconvolve(variable, cnv/2, mode='valid')
    return out

# open file
hisfile = h5py.File(fname,'r')
inifile = netCDF4.Dataset(iname,'r')

# get roms grid using pyroms
rgrid = pyroms.grid.get_ROMS_grid('id',hist_file=fname,grid_file=fname)
dz = diff(rgrid.vgrid.z_w[:,:,:], axis=0)

# get variables I need
ubar = hisfile['ubar'][tind,:,:]
vbar = hisfile['vbar'][tind,:,:]
beta = inifile.variables['phys.beta'][:]

# depth integrate - doesn't work because rvor is incorrect in the first place
#rvorbar = sum(rvor * avg1(avg1(dz/rgrid.vgrid.h,1),2), axis=0)

# Calculate rvorbar terms
Uy = diff(ubar,n=1,axis=0)
Vx = diff(vbar,n=1,axis=1)
rvorbar = Vx - Uy;

Load cartesian grid from file

Vorticity Export

Plot area-averaged vorticity exported onto the shelf and look at the pattern. I need to calculate pvflux = $ \frac{1}{H} \frac{1}{L_x}\int \int_{-H}^{\zeta} v q \, dz \, dx $ where, q is Ertel PV.

In [3]:
isb = 20 ; # must add code to figure this out
hsb = hisfile['h'][26,1]
Lx = amax(rgrid.hgrid.x_rho[:])

# get first eddy shedding
fname2 = dirname + ''
pvname = dirname + ''

eddy = loadmat(edname)['eddy']


zr = rgrid.vgrid.z_r[:]
zr = zr[:,0:isb+5,:]

pvfile = h5py.File(pvname,'r');

dV = avg1(1500 * 1500 * dz,0);
v = hisfile['v'][trange,:,0:isb+5,1:-1]
pv = pvfile['pv'][trange,:,0:isb+5,:]
rv = avg1(pvfile['rv'][trange,:,0:isb+6,:],3)

tpv = pvfile['ocean_time'][trange]

rvflux = sum(sum( rv[:,:,isb,:] * avg1(v[:,:,isb,:]*dz[:,isb,1:-1],1),axis=1),axis=1)/hsb/Lx
pvflux = sum(sum( pv[:,:,isb,:] * avg1(v[:,:,isb,:],1)*avg1(dz[:,isb,1:-1],0),axis=1),axis=1)/hsb/Lx

shelfvol = sum(dV[:,0:isb,:][:])
pvshelf = sum(sum(sum( pv[:,:,0:isb,:] * dV[:,0:isb,1:-1], axis=1), axis=1), axis=1)/shelfvol
rvshelf = sum(sum(sum( rv[:,:,0:isb,:] * dV[:,0:isb,1:-1], axis=1), axis=1), axis=1)/shelfvol

In [73]:
plt.xlabel('Time (days)')
plt.ylabel('Area - averaged fluxes')
plt.legend( ('Ertel PV','Rel Vor.'), loc='lower left')

plt.plot(tpv/86400,pvshelf - pvshelf[0], tpv/86400,rvshelf - rvshelf[0])
plt.ylabel('Volume average vorticity on shelf')
plt.xlabel('Time (days)');

Barotropic Instability

In [6]:


<matplotlib.colorbar.Colorbar instance at 0x5289050>

Now compare $U_{yy}$ to $\beta$

In [7]:
Uyy = diff(ubar,n=2,axis=0)/1500**2
plt.contourf((beta - Uyy)[0:100,:])


<matplotlib.colorbar.Colorbar instance at 0x51583b0>

In [337]:
plt.plot( (beta - Uyy[:,100]).transpose()*1e7 )
plt.xticks(list(plt.xticks()[0]) + [26,59])
plt.legend( ('V_x','U_y','beta - U_yy*1e7','U*0.1') )
plt.xlabel('Y co-ordinate')

<matplotlib.text.Text at 0x497cafd0>

Line plots of $U_y, \beta - U_{yy}, V_x$ . As expected, $V_x \ll U_y$, and so, I can idealize this as a barotropic along-shore jet that is pretty much centered over the slope. Shelfbreak is at 26 and slopebreak at 59

PV gradient $\beta - U_{y}$ crosses zero - Rayleigh's criterion is satisfied. The vorticity maximum is also within the domain, and not at the boundary or inifnity - hence, Fjortoft's criterion is also satisfied.

Barotropic instability is possible! Note that $\beta$ cannot be ignored.

For the theoretical derivation, let's assume that I have a Gaussian velocity profile

In [325]:
y = linspace(-pi/2,pi/2,120)

# parameters
U0 = -1
L = 5

# functional forms
U = cos(2*pi*y)**2 #U0 * exp(-(y/L)**2);
Uy = diff(U,n=1)/diff(y)
Uyy = diff(U,n=2)/avg1(diff(y)**2,0)
y1 = avg1(y,0)
y2 = avg1(y1,0)

# plot
plt.plot(y,U, 'c', y1,Uy,'g', y2,-Uyy,'r')
plt.legend( ('U','Uy','-Uyy') )

<matplotlib.legend.Legend at 0x41ddb750>

I need an analytical expression for the growth rate now. Then code it. Pick out appropriate values for the Gaussian assumption and then see if the estimates match.

In [ ]: