Finite divergence for bolus velocity in the first two verical layer near land

[MaceKuailv]

Hi!

I am using the ECCO v4r4 product and realized that the bolus velocity is not purely divergent-free. I think this is the repo that hosts the code that is responsible for converting GMPsiX and GMPsiY to UVELSTAR, VVELSTAR, WVELSTAR (https://github.com/MITgcm/gcmfaces/blob/master/gcmfaces_calc/calc_bolus.m#L1).

The finite divergence appears at grids that have exactly 2 layers in depth around land boundaries, and it only appears in the first two layers. This is what the convergence of bolus velocity looks like:

Here is my code to go from downloading to producing the above plot. (Excuse me for not able to write in MATLAB)

# Download one daily mean data
# !wget https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/ECCO_L4_BOLUS_LLC0090GRID_DAILY_V4R4/OCEAN_BOLUS_VELOCITY_day_mean_1992-01-01_ECCO_V4r4_native_llc0090.nc

import xarray as xr
import numpy as np
import xgcm
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

grd = xr.open_dataset('~/ECCO-grid/ECCO-GRID.nc')
bolus = xr.open_dataset('OCEAN_BOLUS_VELOCITY_day_mean_1992-01-01_ECCO_V4r4_native_llc0090.nc')

ds = xr.merge([bolus,grd])

# Calculate the transport
# NaN means no transport across rocks, fill them as 0
ds['utrans'] = (ds['UVELSTAR']*ds['drF']*ds['dyG']).fillna(0)
ds['vtrans'] = (ds['VVELSTAR']*ds['drF']*ds['dxG']).fillna(0)
ds['wtrans'] = (ds['WVELSTAR']*ds['rA']).fillna(0)

# Creating the xgcm object for the divergence calculation
face_connections = {
    "tile": {
        0: {"i": ((12, "j", False), (3, "i", False)), "j": (None, (1, "j", False))},
        1: {"i": ((11, "j", False), (4, "i", False)),"j": ((0, "j", False), (2, "j", False))},
        2: {"i": ((10, "j", False), (5, "i", False)),"j": ((1, "j", False), (6, "i", False))},
        3: {"i": ((0, "i", False), (9, "j", False)), "j": (None, (4, "j", False))},
        4: {"i": ((1, "i", False), (8, "j", False)),"j": ((3, "j", False), (5, "j", False))},
        5: {"i": ((2, "i", False), (7, "j", False)),"j": ((4, "j", False), (6, "j", False))},
        6: {"i": ((2, "j", False), (7, "i", False)),"j": ((5, "j", False), (10, "i", False)),},
        7: {"i": ((6, "i", False), (8, "i", False)),"j": ((5, "i", False), (10, "j", False)),},
        8: {"i": ((7, "i", False), (9, "i", False)),"j": ((4, "i", False), (11, "j", False)),},
        9: {"i": ((8, "i", False), None), "j": ((3, "i", False), (12, "j", False))},
        10: {"i": ((6, "j", False), (11, "i", False)),"j": ((7, "j", False), (2, "i", False))},
        11: {"i": ((10, "i", False), (12, "i", False)),"j": ((8, "j", False), (1, "i", False))},
        12: {"i": ((11, "i", False), None),"j": ((9, "j", False), (0, "i", False))},
    }
}

xgcmgrd = xgcm.Grid(
    ds,
    periodic=False,
    face_connections=face_connections,
    coords={
        "i": {"center": "i", "left": "i_g"},
        "j": {"center": "j", "left": "j_g"},
        "k": {"center": "k", "left": "k_l"},
    },
)

# Calculating the derivative (actually convergence)
xy_diff = xgcmgrd.diff_2d_vector(
    {"i": ds['utrans'], "j": ds['vtrans']}, boundary="fill", fill_value=0.0
)
x_diff = xy_diff["i"]
y_diff = xy_diff["j"]
hConv = -(x_diff + y_diff) 

vConv = xgcmgrd.diff(ds['wtrans'], "k", boundary="fill", fill_value=0.0) 

# Combine the convergences and get rid of the time dimension
conv = (hConv+vConv).squeeze()

# Plot the convergence in face 10 of the second level (k = 1)
ax = plt.axes(projection = ccrs.PlateCarree())
ct = ax.pcolormesh(ds.XC[10],ds.YC[10],conv[1,10])
ax.coastlines()
plt.colorbar(ct)

For comparison, the divergence should be many orders of magnitude smaller. Here is what it looks like in the third layer.

The code for reproducing this is:

# Plot the convergence for the level beneath it. 
ax = plt.axes(projection = ccrs.PlateCarree())
ct = ax.pcolormesh(ds.XC[10],ds.YC[10],conv[2,10])
ax.coastlines()
plt.colorbar(ct)

Although I don't know how to write MATLAB, please let me know how I can help.

[MaceKuailv]

BTW, it does not make a difference if one includes hFacS and hFacW

ds['utrans'] = (ds['UVELSTAR']*ds['drF']*ds['dyG']*ds['hFacW']).fillna(0)
ds['vtrans'] = (ds['VVELSTAR']*ds['drF']*ds['dxG']*ds['hFacS']).fillna(0)

[MaceKuailv]

I tried to use xgcm for the same calculation,

strm = xr.open_dataset('OCEAN_BOLUS_STREAMFUNCTION_day_mean_1992-01-01_ECCO_V4r4_native_llc0090.nc')

strmx = strm.GM_PsiX.fillna(0)
strmy = strm.GM_PsiY.fillna(0)

u = xgcmgrd.diff(strmx,"k", boundary = 'fill', fill_value = 0.0)/ds.drF
v = xgcmgrd.diff(strmy,"k", boundary = 'fill', fill_value = 0.0)/ds.drF

vstrmx = strmx*ds.dyG
vstrmy = strmy*ds.dxG

xy_diff = xgcmgrd.diff_2d_vector(
    {"i": vstrmx, "j": vstrmy}, boundary="fill", fill_value=0.0
)
x_diff = xy_diff["i"]
y_diff = xy_diff["j"]
hDiv = (x_diff + y_diff) 

w = hDiv/ds.rA

When I plug this into the divergence calculation above, the divergence is very small. The horizontal components are the same but the vertical velocity at the second level are different.

diffu = (u - ds.UVELSTAR).fillna(0)
diffv = (v - ds.VVELSTAR).fillna(0)
diffw = (w - ds.WVELSTAR).fillna(0)
plt.imshow(diffw[0,1,10].T)
plt.colorbar()

which yields:

[gaelforget]

@ifenty @owang1 can you help with this? It relates to python toolboxes and ECCOv4r4 which you know better than me

[MaceKuailv]

Thank you! Let me know if there is anything that I could help.

[owang01]

Thanks @gaelforget

@MaceKuailv Thank you for the detailed description and the Python code regarding the issue. I will look into it and get back to you.

[ifenty]

Thank you! Let me know if there is anything that I could help.

@MaceKuailv sorry that you didn't get a reply sooner. From what I understand, the issue likely comes from either how the vertical bolus velocity is calculated or how it is prepared/distributed.

You mentioned that the issue only appears ... "_ at grids that have exactly 2 layers in depth around land boundaries, and it only appears in the first two layers.
_" which is a very important clue.

@owang01 does the model calculate these bolus terms or do we calculate it offline?

[MaceKuailv]

I think MITgcm does not have uvelstar as an output variable, The ones that are available are GMpsiX and GMpsiY. Here is a python script of mine to calculate the bolus velocity with xgcm. This one does not create a divergent velocity field.

https://github.com/MaceKuailv/seaduck/blob/948557aac8423393b8254e2439c214538e7c5983/seaduck/eulerian_budget.py#L173-L193

[MaceKuailv]

The code is translated from the code in this repo, so I think the MATLAB tool box probably does not have the same bug.

[ifenty]

I think MITgcm does not have uvelstar as an output variable, The ones that are available are GMpsiX and GMpsiY. Here is a python script of mine to calculate the bolus velocity with xgcm. This one does not create a divergent velocity field.

https://github.com/MaceKuailv/seaduck/blob/948557aac8423393b8254e2439c214538e7c5983/seaduck/eulerian_budget.py#L173-L193

So one can calculate divergence-free bolus velocities from the GMxx output fields from the model, but it seems the code we use has a bug that appears in the case of exactly two depth levels next to land. Is that your understanding?

[MaceKuailv]

Yes, exactly.

[owang01]

We calculated these bolus terms offline. It is likely that there was a bug in calculating these bolus terms. I was able to reproduce the issue that @MaceKuailv found but did not find the root cause of the the bug at that time. I will take one more crack at it.

[ifenty]

We calculated these bolus terms offline. It is likely that there was a bug in calculating these bolus terms. I was able to reproduce the issue that @MaceKuailv found but did not find the root cause of the the bug at that time. I will take one more crack at it.

@owang01 Mace seems to calculate the w that yields a zero divergence. I don't fully understand the calculation (should hFac* be involved?) because it seems zero horizontal divergence yields a zero w. Shouldn't it be zero horizontal divergence yields zero vertical divergence, w_k = w_{k+1}?