<jupyter_start><jupyter_text>Water volume and tracer mass balance through the faces of the shelf boxThe plot at the end of this notebook shows a nice balanced shelf. Although the transport is slightly offshore throughout the simulation it can be explained by the drop in sea surface height.<jupyter_code>import cmocean as cmo
import matplotlib.pyplot as plt
%matplotlib inline
from netCDF4 import Dataset
import numpy as np
import seaborn as sns
import canyon_tools.readout_tools as rout 
import canyon_tools.metrics_tools as mpt
sns.set_context('notebook')
sns.set_style('whitegrid')
def unstagger_xarray(qty, index):
    """Interpolate u, v, or w component values to values at grid cell centres.
    
    Named indexing requires that input arrays are XArray DataArrays.

    :arg qty: u, v, or w component values
    :type qty: :py:class:`xarray.DataArray`
    
    :arg index: index name along which to centre
        (generally one of 'gridX', 'gridY', or 'depth')
    :type index: str

    :returns qty: u, v, or w component values at grid cell centres
    :rtype: :py:class:`xarray.DataArray`
    """
    qty = (qty + qty.shift(**{index: 1})) / 2
    
    return qty

# No-canyon files
Grid = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/04_Ast03_No_Cny_Argo/gridGlob.nc' 
Ptracers = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/04_Ast03_No_Cny_Argo/ptracersGlob.nc' 
State = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/04_Ast03_No_Cny_Argo/stateGlob.nc' 
Flux = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/04_Ast03_No_Cny_Argo/FluxTR02Glob.nc'
# Define shelf box indices
# E = east boundary, 
# W = west boundary, 
# N = north boundary (open coast), 
# S = south boundary (virtual wall from surface to shelf break).

zind = 30
Zind = slice(0,zind)

Nx = slice(0,616) # open
Ny = 358           

Sx = slice(0,616) # open
Sy = 131

Wx = 0
Wy = slice(131,359) # periodic

Ex = 616
Ey = slice(131,359) # periodic

Lidx = slice(0,616)
Lidy = slice(131,359)
Lidz = zind
# Grid variables - these are the same for N and S
with Dataset(Grid, 'r') as nbl:
    dxg = nbl.variables['dxG'][Sy,Sx]
    dyg = nbl.variables['dyG'][Ey,Ex]
    drf = nbl.variables['drF'][Zind]
    depth = nbl.variables['Z'][Zind]
    X = nbl.variables['X'][:]
    Y = nbl.variables['Y'][:]
    ra = nbl.variables['rA'][Lidy,Lidx]
    hfacsS = nbl.variables['HFacS'][:zind,Sy,Sx]
    hfacsN = nbl.variables['HFacS'][:zind,Ny,Nx]


dxg_exp = np.expand_dims(dxg,0)
dyg_exp = np.expand_dims(dyg,0)
drf_exp = np.expand_dims(drf,1)
dxg_exp = dxg_exp + np.zeros(np.shape(hfacsS))
dyg_exp = dyg_exp + np.zeros((len(range(zind)),len(range(359-Sy))))
drf_exp = drf_exp + np.zeros(np.shape(hfacsS))
# Calculate transport as u_i*A_i where u_i is the perpendicular velocity to the area A_i for cell i.
# East and west balance each other becasue bcs are periodic so I won't add them.
#
# Do I need to unstagger u,v?
#

# Shelf break
with Dataset(State, 'r') as nbl:
    v_south = nbl.variables['V'][:,Zind,Sy,Sx]
    S_trans = v_south*dxg_exp*hfacsS*drf_exp    
    v_north = nbl.variables['V'][:,Zind,Ny,Nx]
    N_trans = v_north*dxg_exp*hfacsN*drf_exp 
    u_east = nbl.variables['U'][:,Zind,Ey,Ex]
    E_trans = u_east*dyg_exp*drf_exp[:,Ey]
    u_west = nbl.variables['U'][:,Zind,Wy,Wx]
    W_trans = u_west*dyg_exp*drf_exp[:,Ey] 
    
    eta = nbl.variables['Eta'][:,Lidy,Lidx]     # approximation of volume change due to ssh drop
    eta_trans = ((eta[1:]-eta[:-1])*ra)/(3600*12)
    
    w_lid = nbl.variables['W'][:,Lidz,Lidy,Lidx]
    Lid_trans = w_lid*ra


S_trans_ts = np.sum(np.sum(S_trans,axis=1),axis=1)
N_trans_ts = np.sum(np.sum(N_trans,axis=1),axis=1)
E_trans_ts = np.sum(np.sum(E_trans,axis=1),axis=1)
W_trans_ts = np.sum(np.sum(W_trans,axis=1),axis=1)

eta_trans_ts = np.sum(np.sum(eta_trans,axis=1),axis=1)
Lid_trans_ts = np.sum(np.sum(Lid_trans,axis=1),axis=1)

fig,ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(np.arange(19)/2,W_trans_ts-E_trans_ts, label='W-E')
ax.plot(np.arange(1,19)/2,eta_trans_ts, label='SSH dropping')
ax.plot(np.arange(19)/2,-N_trans_ts,label='Open coast')
ax.plot(np.arange(19)/2,S_trans_ts,label='Shelf wall')
ax.plot(np.arange(19)/2,W_trans_ts-E_trans_ts-N_trans_ts+S_trans_ts,'k',
        label='net transport')

ax.set_title('Water transport through shelf box walls')
ax.set_xlabel('days')
ax.set_ylabel('Transport (m$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')<jupyter_output><empty_output><jupyter_text>SSH drop was not added to the net transport because the water is not actually leaving through the surface. I just wanted to see if it was enough to account for the transport through the shelf wall. Note that all the water that leaves the shelf (Shelf wall) is the same as the water transport due to the sea surface height drop.   Tracer balance (Linear tracer)<jupyter_code># Grid variables - these are the same for N and S
# East and west balance each other becasue bcs are periodic so I won't add them.
# I'm not counting diffusive flux

# Shelf break
with Dataset(Flux, 'r') as nbl:
    v_south_tr = nbl.variables['ADVyTr02'][:,Zind,Sy,Sx]
    S_trans_tr = v_south_tr    # I didn't multiply by hfacS in my previous calculations
    v_north_tr = nbl.variables['ADVyTr02'][:,Zind,Ny,Nx]
    N_trans_tr = v_north_tr
    u_east_tr = nbl.variables['ADVxTr02'][:,Zind,Ey,Ex]
    E_trans_tr = u_east_tr
    u_west_tr = nbl.variables['ADVxTr02'][:,Zind,Wy,Wx]
    W_trans_tr = u_west_tr
    
with Dataset(Ptracers,'r') as nbl: 
    conc_surf = nbl.variables['Tr02'][:,0,Lidy,Lidx] # pick surface level
    surf_trans_tr = (((eta[1:]*conc_surf[1:])-(eta[:-1]*conc_surf[:-1]))*ra)/(3600*12)

eta_trans_tr_ts = np.sum(np.sum(surf_trans_tr,axis=1),axis=1)
S_trans_tr_ts = np.sum(np.sum(S_trans_tr,axis=1),axis=1)
N_trans_tr_ts = np.sum(np.sum(N_trans_tr,axis=1),axis=1)
E_trans_tr_ts = np.sum(np.sum(E_trans_tr,axis=1),axis=1)
W_trans_tr_ts = np.sum(np.sum(W_trans_tr,axis=1),axis=1)
fig, ax = plt.subplots(1,2,figsize=(12,4))
mesh = ax[0].pcolor(v_south_tr[7,:,:], cmap=cmo.cm.delta)
plt.colorbar(mesh, ax=ax[0])
mesh2 = ax[1].pcolor(hfacsS, cmap=cmo.cm.thermal)
plt.colorbar(mesh2, ax=ax[1])
fig,ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(np.arange(1,19)/2,W_trans_tr_ts-E_trans_tr_ts, label='W-E')
ax.plot(np.arange(1,19)/2,eta_trans_tr_ts,label='SSH drop')
ax.plot(np.arange(1,19)/2,-N_trans_tr_ts,label='Open coast')
ax.plot(np.arange(1,19)/2,S_trans_tr_ts,label='Shelf wall')
ax.plot(np.arange(1,19)/2,W_trans_tr_ts-E_trans_tr_ts+S_trans_tr_ts-N_trans_tr_ts,'k',label='net transport')

ax.set_title('Tracer 02 transport through shelf box walls')
ax.set_xlabel('days')
ax.set_ylabel('Transport ($\mu$Mm$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')<jupyter_output><empty_output><jupyter_text>Tracer transport is not zero, there is tracer coming onto the shelf throughout the whole simulation (positive transport through shelf wall)  Canyon case<jupyter_code># Canyon files
GridCan = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/03_Ast03_Argo/gridGlob.nc' 
PtracersCan = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/03_Ast03_Argo/ptracersGlob.nc' 
StateCan = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/03_Ast03_Argo/stateGlob.nc' 
FluxCan = '/data/kramosmu/results/TracerExperiments/UPW_10TR_BF2_AST/03_Ast03_Argo/FluxTR02Glob.nc'
# Grid variables - these are the same for N and S
with Dataset(GridCan, 'r') as nbl:
    dxg = nbl.variables['dxG'][Sy,Sx]
    dyg = nbl.variables['dyG'][Ey,Ex]
    drf = nbl.variables['drF'][Zind]
    depth = nbl.variables['Z'][Zind]
    X = nbl.variables['X'][:]
    Y = nbl.variables['Y'][:]
    ra = nbl.variables['rA'][Lidy,Lidx]
    hfacsS = nbl.variables['HFacS'][:zind,Sy,Sx]
    hfacsN = nbl.variables['HFacS'][:zind,Ny,Nx]


dxg_exp = np.expand_dims(dxg,0)
dyg_exp = np.expand_dims(dyg,0)
drf_exp = np.expand_dims(drf,1)
dxg_exp = dxg_exp + np.zeros(np.shape(hfacsS))
dyg_exp = dyg_exp + np.zeros((len(range(zind)),len(range(359-Sy))))
drf_exp = drf_exp + np.zeros(np.shape(hfacsS))
# Calculate transport as u_i*A_i where u_i is the perpendicular velocity to the area A_i for cell i.
# East and west balance each other becasue bcs are periodic so I won't add them.
#
# Do I need to unstagger u,v?
#

# Shelf break
with Dataset(StateCan, 'r') as nbl:
    v_southCan = nbl.variables['V'][:,Zind,Sy,Sx]
    S_transCan = v_southCan*dxg_exp*hfacsS*drf_exp    
    v_northCan = nbl.variables['V'][:,Zind,Ny,Nx]
    N_transCan = v_northCan*dxg_exp*hfacsN*drf_exp 
    u_eastCan = nbl.variables['U'][:,Zind,Ey,Ex]
    E_transCan = u_eastCan*dyg_exp*drf_exp[:,Ey]
    u_westCan = nbl.variables['U'][:,Zind,Wy,Wx]
    W_transCan = u_westCan*dyg_exp*drf_exp[:,Ey] 
    
    etaCan = nbl.variables['Eta'][:,Lidy,Lidx]     # approximation of volume change due to ssh drop
    eta_transCan = ((etaCan[1:]-etaCan[:-1])*ra)/(3600*12)
    
    w_lidCan = nbl.variables['W'][:,Lidz,Lidy,Lidx]
    Lid_transCan = w_lidCan*ra


S_transCan_ts = np.sum(np.sum(S_transCan,axis=1),axis=1)
N_transCan_ts = np.sum(np.sum(N_transCan,axis=1),axis=1)
E_transCan_ts = np.sum(np.sum(E_transCan,axis=1),axis=1)
W_transCan_ts = np.sum(np.sum(W_transCan,axis=1),axis=1)

eta_transCan_ts = np.sum(np.sum(eta_transCan,axis=1),axis=1)
Lid_transCan_ts = np.sum(np.sum(Lid_transCan,axis=1),axis=1)

fig,ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(np.arange(19)/2,W_transCan_ts-E_transCan_ts, label='W-E')
ax.plot(np.arange(1,19)/2,eta_transCan_ts, label='SSH dropping')
ax.plot(np.arange(19)/2,-N_transCan_ts,label='Open coast')
ax.plot(np.arange(19)/2,S_transCan_ts,label='Shelf wall')
ax.plot(np.arange(19)/2,Lid_transCan_ts,label='Vertical')
ax.plot(np.arange(19)/2,W_transCan_ts-E_transCan_ts-N_transCan_ts+S_transCan_ts+Lid_transCan_ts,'k',
        label='net transport')

ax.set_title('Water transport through shelf box walls canyon case')
ax.set_xlabel('days')
ax.set_ylabel('Transport (m$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')
fig,ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(np.arange(19)/2,W_transCan_ts-E_transCan_ts, label='W-E')
ax.plot(np.arange(1,19)/2,eta_transCan_ts, label='SSH dropping')
ax.plot(np.arange(19)/2,N_transCan_ts,label='Open coast')
ax.plot(np.arange(19)/2,Lid_transCan_ts+S_transCan_ts,label='Vertical+Shelf Wall')
ax.plot(np.arange(19)/2,W_transCan_ts-E_transCan_ts-N_transCan_ts+S_transCan_ts+Lid_transCan_ts,'k',
        label='net transport')

ax.set_title('Water transport through shelf box walls canyon case')
ax.set_xlabel('days')
ax.set_ylabel('Transport (m$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')
with Dataset(FluxCan, 'r') as nbl:
    v_southCan_tr = nbl.variables['ADVyTr02'][:,Zind,Sy,Sx]
    S_transCan_tr = v_southCan_tr    # I didn't multiply by hfacS in my previous calculations
    v_northCan_tr = nbl.variables['ADVyTr02'][:,Zind,Ny,Nx]
    N_transCan_tr = v_northCan_tr 
    u_eastCan_tr = nbl.variables['ADVxTr02'][:,Zind,Ey,Ex]
    E_transCan_tr = u_eastCan_tr
    u_westCan_tr = nbl.variables['ADVxTr02'][:,Zind,Wy,Wx]
    W_transCan_tr = u_westCan_tr 
    w_lidCan_tr = nbl.variables['ADVrTr02'][:,Lidz,Lidy,Lidx]
    Lid_transCan_tr = w_lidCan_tr

with Dataset(PtracersCan,'r') as nbl: 
    conc_surfCan = nbl.variables['Tr02'][:,0,Lidy,Lidx] # pick surface level
    surf_transCan_tr = (((etaCan[1:]*conc_surfCan[1:])-(etaCan[:-1]*conc_surfCan[:-1]))*ra)/(3600*12)

eta_transCan_tr_ts = np.sum(np.sum(surf_transCan_tr,axis=1),axis=1)
S_transCan_tr_ts = np.sum(np.sum(S_transCan_tr,axis=1),axis=1)
N_transCan_tr_ts = np.sum(np.sum(N_transCan_tr,axis=1),axis=1)
E_transCan_tr_ts = np.sum(np.sum(E_transCan_tr,axis=1),axis=1)
W_transCan_tr_ts = np.sum(np.sum(W_transCan_tr,axis=1),axis=1)
lid_transCan_tr_ts = np.sum(np.sum(Lid_transCan_tr,axis=1),axis=1)


fig, ax = plt.subplots(1,2,figsize=(12,4))
mesh = ax[0].pcolor(v_south_tr[7,:,:], cmap=cmo.cm.delta)
plt.colorbar(mesh, ax=ax[0])
mesh2 = ax[1].pcolor(hfacsS, cmap=cmo.cm.thermal)
plt.colorbar(mesh2, ax=ax[1])
fig,ax = plt.subplots(1,1,figsize=(6,4))
ax.plot(np.arange(1,19)/2,W_transCan_tr_ts-E_transCan_tr_ts, label='W-E')
ax.plot(np.arange(1,19)/2,eta_transCan_tr_ts,label='SSH drop')
ax.plot(np.arange(1,19)/2,N_transCan_tr_ts,label='Open coast')
ax.plot(np.arange(1,19)/2,S_transCan_tr_ts,label='Shelf wall')
ax.plot(np.arange(1,19)/2,lid_transCan_tr_ts,label='Vertical')

ax.plot(np.arange(1,19)/2,W_transCan_tr_ts-E_transCan_tr_ts+S_transCan_tr_ts-N_transCan_tr_ts+
        lid_transCan_tr_ts,'k',label='net transport')

ax.set_title('Tracer transport through shelf box walls canyon case')
ax.set_xlabel('days')
ax.set_ylabel('Transport ($\mu$Mm$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')
fig,ax = plt.subplots(1,1,figsize=(6,4))

ax.plot(np.arange(1,19)/2,W_transCan_tr_ts-E_transCan_tr_ts+S_transCan_tr_ts-N_transCan_tr_ts+
        lid_transCan_tr_ts,label='net transport canyon')
ax.plot(np.arange(1,19)/2,W_trans_tr_ts-E_trans_tr_ts+S_trans_tr_ts-N_trans_tr_ts
        ,label='net transport no canyon')
ax.plot(np.arange(1,19)/2,(W_transCan_tr_ts-E_transCan_tr_ts+S_transCan_tr_ts-N_transCan_tr_ts+
        lid_transCan_tr_ts)-(W_trans_tr_ts-E_trans_tr_ts+S_trans_tr_ts-N_trans_tr_ts),
        'k',label='net transport anomaly')


ax.set_title('Tracer transport through shelf box walls')
ax.set_xlabel('days')
ax.set_ylabel('Transport ($\mu$Mm$^3$s$^{-1}$)')
ax.legend(loc=0)
ax.axhline(0,color='0.5',linestyle='--')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16153334810547307, 'fraction_numerical': 0.030799911367161533, 'mean_word_length': 5.739671478347437, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '171653', 'license': 'Apache-2.0', 'n_tokens_mistral': 6119, 'n_tokens_neox': 5915, 'n_words': 911, 'path': 'forPaper2/balance_water_on_shelf_Ast_ARGO.ipynb', 'repo_name': 'UBC-MOAD/outputanalysisnotebooks', 'size': '13539'}