This notebook looks at the vertical eddy viscosity/diffusivity during a deep water renewal event in late August 2003.
This is a 10 day simulation with daily averaged output.
In [1]:
%matplotlib inline
from matplotlib import pylab
import matplotlib.pyplot as plt
import netCDF4 as NC
import numpy as np
from salishsea_tools import (nc_tools,viz_tools)
In [2]:
# Load the data. Path name can be changed to look at different data.
path = '/data/nsoontie/MEOPAR/SalishSea/results/BBL/eddy_noBBL/SalishSea_1d_20030819_20030828_grid_T.nc'
f = NC.Dataset(path,'r');
sal=f.variables['vosaline']
depths = f.variables['deptht']
T_lat = f.variables['nav_lat']
T_lon = f.variables['nav_lon']
#Loading eddy viscosity/diffusivity data on the vertical grid
path = '/data/nsoontie/MEOPAR/SalishSea/results/BBL/eddy_noBBL/SalishSea_1d_20030819_20030828_grid_W.nc'
f = NC.Dataset(path,'r');
avm=f.variables['ve_eddy_visc']
avd = f.variables['ve_eddy_diff'] #
W=f.variables['vovecrtz']
W_lat = f.variables['nav_lat']
W_lon = f.variables['nav_lon']
depthw = f.variables['depthw']
#Loading data on the ugrid
path = '/data/nsoontie/MEOPAR/SalishSea/results/BBL/eddy_noBBL/SalishSea_1d_20030819_20030828_grid_U.nc'
f = NC.Dataset(path,'r');
U=f.variables['vozocrtx']
U_lat = f.variables['nav_lat']
U_lon = f.variables['nav_lon']
#Loading data on the ugrid
path = '/data/nsoontie/MEOPAR/SalishSea/results/BBL/eddy_noBBL/SalishSea_1d_20030819_20030828_grid_V.nc'
f = NC.Dataset(path,'r');
V=f.variables['vomecrty']
V_lat = f.variables['nav_lat']
V_lon = f.variables['nav_lon']
grid = NC.Dataset('/data/nsoontie/MEOPAR/NEMO-forcing/grid/bathy_meter_SalishSea2.nc')
bathy=grid.variables['Bathymetry']
Thalweg¶
Plotting salinity and eddy viscosity/diffusivity along thalweg over time. Daily average outputs over 10 days.
In [3]:
lines = np.loadtxt('/data/nsoontie/MEOPAR/tools/analysis_tools/thalweg.txt', delimiter=" ", unpack=False)
lines = lines.astype(int)
thalweg_lon = T_lon[lines[:,0],lines[:,1]]
thalweg_lat = T_lat[lines[:,0],lines[:,1]]
ds=np.arange(0,lines.shape[0],1);
vs=np.arange(34,27.5,-0.5);
XX_T,ZZ_T = np.meshgrid(ds,-depths[:])
XX_W,ZZ_W = np.meshgrid(ds,-depthw[:])
Below I will plot the salinity field, and eddy viscosity/diffusivity along the thalweg. Each row of plots corresponds to a daily average. A deep water renewal event occurs around t=6 (the 30.5 salinity contour is slightly elevated on the left hand slope of the Strait of Georgia basin).
In [4]:
smin=28; smax=34
emin=-4; emax=2
(fig,axs)=plt.subplots(10,3,figsize=(25,25),sharey=True)
ts=np.arange(10)
for t,axS,axV,axD in zip(ts,axs[:,0],axs[:,1],axs[:,2]):
#salinity
salP=sal[t,:,lines[:,0],lines[:,1]];
mu =salP == 0; salP= np.ma.array(salP,mask=mu)
mesh=axS.pcolormesh(XX_T,ZZ_T,salP,vmin=smin,vmax=smax,cmap='rainbow')
CS=axS.contour(XX_T,ZZ_T,salP,vs, colors='black')
axS.clabel(CS,fontsize=9, inline=1)
axS.set_title('t = ' +str(t))
axS.set_ylim([-450,0]); axS.set_xlim([0,700])
#viscosity
avmP=avm[t,:,lines[:,0],lines[:,1]];
mu =avmP == 0; avmP= np.ma.array(avmP,mask=mu)
mesh=axV.pcolormesh(XX_W,ZZ_W,np.log10(avmP),vmin=emin, vmax=emax,cmap='hot')
axV.set_title('t = ' +str(t))
axV.set_ylim([-450,0]); axV.set_xlim([0,700])
#diffusivity
avdP=avd[t,:,lines[:,0],lines[:,1]];
mu =avdP == 0; avdP= np.ma.array(avdP,mask=mu)
mesh=axD.pcolormesh(XX_W,ZZ_W,np.log10(avdP),vmin=emin,vmax=emax,cmap='hot')
axD.set_title('t = ' +str(t))
axD.set_ylim([-450,0]); axD.set_xlim([0,700])
fig.subplots_adjust(wspace=0.05)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.82, 0.25, 0.02, 0.5])
cbar=fig.colorbar(mesh, cax=cbar_ax)
cbar.set_label('Log of Vertical eddy coefficient (m^2/s)')
print 'Left: Salinity'
print 'Middle: Eddy Viscosity'
print 'Right: Eddy Diffusivity'
-c:20: RuntimeWarning: divide by zero encountered in log10 -c:27: RuntimeWarning: divide by zero encountered in log10
Left: Salinity Middle: Eddy Viscosity Right: Eddy Diffusivity
The path of the talweg is shown below.
In [5]:
fig,ax=plt.subplots(1,1,figsize=(4,8))
ax.pcolormesh(T_lon[:], T_lat[:], bathy, cmap='winter_r')
ax.plot(thalweg_lon,thalweg_lat,linestyle='-', marker='+', color='red',
label='Thalweg Points')
legend = ax.legend(loc='best', fancybox=True, framealpha=0.25)
ax.axis([-126,-122,47,51])
Out[5]:
[-126, -122, 47, 51]
Profiles¶
Next, look at vertical profiles of the salinity and eddy coefficients. First, pick out a few interesting spots on the map.
In [6]:
fig,axs=plt.subplots(1,2,figsize=(8,4))
cmap = plt.get_cmap('winter_r')
land_colour = 'burlywood'
cmap.set_bad(land_colour)
ax=axs[0]
mesh=ax.pcolormesh(bathy[:],cmap=cmap)
fig.colorbar(mesh,ax=ax)
ax.set_xlim([200,350]); ax.set_ylim([200,500])
ax.set_title('Bathymetry')
ii=[300,280,300,270,245,240]
jj=[400,400,350,340,340,300]
colors=['b','g','r','c','m','y']
for (i,j,c) in zip(ii,jj,colors):
ax.plot(i,j,'o',color=c)
ax=axs[1]
t=9; depth=0;
salP=sal[t,:,:,:]
mu =salP == 0; salPl= np.ma.array(salP,mask=mu)
cmap = plt.get_cmap('YlGnBu')
cmap.set_bad(land_colour)
mesh=ax.pcolormesh(salPl[depth,:,:],cmap=cmap)
cbar=fig.colorbar(mesh,ax=ax)
ax.set_xlim([200,350]); ax.set_ylim([200,500])
ax.set_title('Salinity')
for (i,j,c) in zip(ii,jj,colors):
ax.plot(i,j,'o',color=c)
Below: salinity and eddy viscosity/diffusivity profiles at the above markers. The colour of the salinity line corresponds to the color of the marker above.
These are plots from the last output time displayed in the thalweg plots.
In [7]:
fig,axs=plt.subplots(2,3,figsize=(10,8))
for (i,j,ax,c) in zip(ii,jj,axs.flat,colors):
ax1=ax
ax1.plot(sal[t,:,j,i],depths[:],'-', color=c)
ax1.set_xlabel('Salinity',color=c)
ax1.set_ylabel('depth (m)')
ax1.set_xlim([0,34]); ax1.set_ylim([0,-450])
for tl in ax1.get_xticklabels():
tl.set_color(c)
ax2=ax1.twiny()
ax2.plot(avm[t,:,j,i],depthw[:],'-k', label='Eddy viscosity')
ax2.plot(avd[t,:,j,i],depthw[:],'--k', label='Eddy diffusivity')
ax2.set_xlabel('Eddy coefficient (m^2/s)')
ax2.set_xlim([0,18]);
ax2.legend(loc=3)
fig.subplots_adjust(hspace=0.4)
fig.subplots_adjust(wspace=0.4)
River Plume¶
Look at the eddy viscosity/diffusivity in the river plume.
In [8]:
fig,axs=plt.subplots(1,2,figsize=(8,4))
cmap = plt.get_cmap('winter_r')
land_colour = 'burlywood'
cmap.set_bad(land_colour)
ax=axs[0]
mesh=ax.pcolormesh(bathy[:],cmap=cmap)
fig.colorbar(mesh,ax=ax)
ax.set_xlim([200,350]); ax.set_ylim([200,500])
ax.set_title('Bathymetry')
ii=[300,300,300,280,310,290]
jj=[400,450,425,420,420,400]
colors=['b','g','r','c','m','y']
for (i,j,c) in zip(ii,jj,colors):
ax.plot(i,j,'o',color=c)
ax=axs[1]
t=9; depth=0;
salP=sal[t,:,:,:]
mu =salP == 0; salPl= np.ma.array(salP,mask=mu)
cmap = plt.get_cmap('YlGnBu')
cmap.set_bad(land_colour)
mesh=ax.pcolormesh(salPl[depth,:,:],cmap=cmap)
cbar=fig.colorbar(mesh,ax=ax)
ax.set_xlim([200,350]); ax.set_ylim([200,500])
ax.set_title('Salinity')
for (i,j,c) in zip(ii,jj,colors):
ax.plot(i,j,'o',color=c)
In [9]:
def profile(i,j):
#plots salinity and eddy coefficent profiles at i,j grid points and t time
fig,axs=plt.subplots(2,3,figsize=(15,10))
depth=0
ax=axs[0,0]
cmap = plt.get_cmap('winter_r')
land_colour = 'burlywood'
cmap.set_bad(land_colour)
mesh=ax.pcolormesh(bathy[:],cmap=cmap)
fig.colorbar(mesh,ax=ax)
ax.set_xlim([200,350]); ax.set_ylim([200,500])
ax.set_title('Bathymetry')
ax.plot(i,j,'o',color='r')
for t in range(10):
ax=axs[1,0]
ax.plot(sal[t,:,j,i],depths[:],'-+',label='t=' +str(t) +'d')
ax.set_xlabel('Salinity')
ax.set_ylabel('depth (m)')
ax.set_xlim([0,34]); ax.set_ylim([10,0])
ax.legend(loc=0)
ax=axs[1,1]
ax.plot(avm[t,:,j,i],depthw[:],'-+', label=t)
ax.set_xlabel('Eddy viscosity (m^2/s)')
ax.set_xlim([-1,10]); ax.set_ylim([10,0])
ax=axs[1,2]
ax.plot(avd[t,:,j,i],depthw[:],'-+', label=t)
ax.set_xlabel('Eddy diffusivity (m^2/s)')
ax.set_xlim([-1,10]); ax.set_ylim([10,0])
return fig, axs
In [10]:
k=4; t=0; i=ii[k]; j=jj[k]
fig,ax=profile(i,j)
In [11]:
k=2; i=ii[k]; j=jj[k]
fig,ax=profile(i,j)
In [12]:
k=4; i=ii[k]; j=jj[k]+10
fig,ax=profile(i,j)
In [13]:
k=0;i=ii[k]; j=jj[k]+10
fig,ax=profile(i,j)
In [14]:
k=5; i=ii[k]; j=jj[k]+10
fig,ax=profile(i,j)
In [15]:
ii=[300,280,300,270,245,240]
jj=[400,400,350,340,340,300]
In [16]:
k=1; i=ii[k];
fig,ax=profile(i,j)
In [17]:
k=3; i=ii[k];
fig,ax=profile(i,j)
In [18]:
k=3; i=ii[k]; j=jj[k]
fig,ax=profile(i,j)
In [19]:
k=4; i=ii[k]; j=jj[k]
fig,ax=profile(i,j)
In [20]:
k=5; i=ii[k]; j=jj[k]
fig,ax=profile(i,j)
In the upper layer (less than 3 m), there is often an increase in eddy diffusicity and viscosity on t=4. This happens at most of the profiles selected. It also looks like this is that start of a deep water renewal event.
Look at tidal cycle.