Part 2: Thermal Wind¶

Andrew Delman, updated 2023-02-05.

Objectives¶

To use ECCO state estimate output to explore thermal wind and its use in understanding the ocean circulation.

By the end of the tutorial, you will be able to:

Plot transects along lines of latitude in the ECCO native grid
Rotate fields in model axes to true zonal/meridional orientation
Relate the density structure to the direction of oceanic flow (thermal wind balance)
Estimate oceanic velocity using thermal wind and assumption of a level of no motion

Introduction¶

In the previous tutorial we established that geostrophic balance explains most of the oceanic circulation–at least at the spatial and temporal scales of ECCOv4r4 output. Now we consider thermal wind, a property of geostrophic flow that has been very useful to atmospheric and ocean scientists alike. Thermal wind balance is derived directly from geostrophic balance (in horizontal directions)

$v_g = \frac{1}{f\rho}\frac{\partial{p}}{\partial{x}}$

$u_g = -\frac{1}{f\rho}\frac{\partial{p}}{\partial{y}}$

and hydrostatic balance (in the vertical direction)

$\frac{\partial{p}}{\partial{z}} = -{\rho}g$

Therefore by differentiating geostrophic balance along the $z$ axis and assuming hydrostatic conditions, we get thermal wind balance:

$\frac{\partial{v_g}}{\partial{z}} = -\frac{g}{f\rho}\frac{\partial{\rho}}{\partial{x}}$

$\frac{\partial{u_g}}{\partial{z}} = \frac{g}{f\rho}\frac{\partial{\rho}}{\partial{y}}$

(cf. Vallis ch. 2.8, Kundu and Cohen ch. 14.5, Gill ch. 7.7)

The thermal wind name originates from the application of these equations to atmospheric dynamics, where the velocities are the wind and density is mostly a function of temperature. But the balance is just as important to oceanic dynamics. In both the atmosphere and ocean, thermal wind enabled early measurements of temperature (and salinity in the ocean) to be used to construct velocity fields, where velocities are often difficult to measure directly. And even in present times when velocities can be measured directly, the thermal wind relation connects changes in heat and salt content to changes in the ocean circulation.

Computing thermal wind balance¶

Thermal wind balance across Atlantic (26 N)¶

Now let’s compute horizontal density gradients (times $g/f$ ) and compare with vertical velocity gradients—i.e., thermal wind balance. This calculation is complicated by the fact that density, $u$ and $v$ all start out at different positions on each grid cell, and that we take differences in both horizontal and vertical directions. To enable this to be done so that both sides of the equations end up on the same points in the grid, the calculations are done in the following order:

Density at center of grid cells (i,j,k) –horiz difference–>
- Horizontal density gradients at edge of grid cells (i_g,j,k) or (i,j_g,k) –horiz interpolate–>
  - Interpolated horizontal density gradients (i,j,k)
$u$ , $v$ at edge of grid cells (i_g,j,k) or (i,j_g,k) –horiz interpolate–>
- Centered $u$ , $v$ (i,j,k) –vert difference–>
  - Vertical $u$ , $v$ gradients (i,j,k_l) –vert interpolate–>
    - Interpolated $u$ , $v$ gradients (i,j,k)

Then we rotate velocities and gradients from model axes to geographical (zonal and meridional) axes. Finally, we apply the lon_depth_along_lat function used above to look at the thermal wind balance across the Atlantic basin.

[10]:

# get xgcm Grid object for ECCO LLC90 grid
xgcm_grid = ecco.get_llc_grid(ds_grid)

[11]:

g = 9.81
# compute f from latitude of grid cell centers
lat = ds_grid.YC
Omega = (2*np.pi)/86164
lat_rad = (np.pi/180)*lat    # convert latitude from degrees to radians
# expand dimensions so that f broadcasts correctly across DataArrays
f = np.expand_dims(2*Omega*np.sin(lat_rad),axis=(0,1))

# compute derivatives of density in x and y
d_rho_dx = (xgcm_grid.diff(dens,axis="X",boundary='extend'))/ds_grid.dxC
d_rho_dy = (xgcm_grid.diff(dens,axis="Y",boundary='extend'))/ds_grid.dyC
# interpolate (vector) gradient values to center of grid cells
import warnings
with warnings.catch_warnings():
    warnings.simplefilter("ignore")    # use this to ignore future warnings caused by interp_2d_vector function
    rho_grads_interp = xgcm_grid.interp_2d_vector({"X":d_rho_dx,"Y":d_rho_dy},boundary='extend')
d_rho_dx = rho_grads_interp['X']
d_rho_dy = rho_grads_interp['Y']
d_rho_dx.name = 'drho_dx'
d_rho_dy.name = 'drho_dy'

# rotate from model to geographical (lon-lat) axes
CS = np.expand_dims(ds_grid.CS,axis=(0,1))
SN = np.expand_dims(ds_grid.SN,axis=(0,1))
d_rho_zonal = (CS*d_rho_dx) - (SN*d_rho_dy)
d_rho_merid = (SN*d_rho_dx) + (CS*d_rho_dy)

# right-hand side of thermal wind balance equations
therm_wind_RHS_1 = -(g/(f*dens))*d_rho_zonal
therm_wind_RHS_2 = (g/(f*dens))*d_rho_merid

[12]:

# load velocity dataset
download_dir = join(download_root_dir,vel_monthly_shortname)
curr_vel_file = list(glob.glob(download_dir + '/*2000-01*.nc'))
ds_vel = xr.open_dataset(curr_vel_file[0])

# interpolate velocities to center of grid cells
ds_vel.UVEL.values[np.isnan(ds_vel.UVEL.values)] = 0
ds_vel.VVEL.values[np.isnan(ds_vel.VVEL.values)] = 0
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    vel_interp = xgcm_grid.interp_2d_vector({'X':ds_vel.UVEL,'Y':ds_vel.VVEL},boundary='extend')
u = vel_interp['X']
v = vel_interp['Y']

# vertical derivatives of velocity
# (negative sign after = because z increases upward, while k indices decrease upward)
du_dz_gridedge = -(u.diff("k").values)/np.expand_dims(ds_grid.drC[1:-1].values,axis=(0,2,3,4))
dv_dz_gridedge = -(v.diff("k").values)/np.expand_dims(ds_grid.drC[1:-1].values,axis=(0,2,3,4))

# vertical interpolation of du/dz, dv/dz
du_dz = np.empty(u.shape)
du_dz.fill(np.nan)    # fill empty numpy array with NaNs
dv_dz = np.empty(v.shape)
dv_dz.fill(np.nan)
du_dz[:,1:-1,:,:,:] = du_dz_gridedge[:,1:,:,:,:] - (np.diff(du_dz_gridedge,axis=1)/2)
dv_dz[:,1:-1,:,:,:] = dv_dz_gridedge[:,1:,:,:,:] - (np.diff(dv_dz_gridedge,axis=1)/2)
# create xarray DataArrays
du_dz_array = xr.DataArray(
                data=du_dz,
                dims=["time","k","tile","j","i"],
                coords=u.coords,
)
dv_dz_array = xr.DataArray(
                data=dv_dz,
                dims=["time","k","tile","j","i"],
                coords=v.coords,
)

# rotate from model to geographical (lon-lat) axes
du_dz_zonal = (CS*du_dz_array) - (SN*dv_dz_array)
dv_dz_merid = (SN*du_dz_array) + (CS*dv_dz_array)

# left-hand side of thermal wind balance equations
therm_wind_LHS_1 = dv_dz_merid
therm_wind_LHS_2 = du_dz_zonal

[13]:

# generate DataArrays along line of latitude
lat_transect = 26.
Atl_W_bound = -82.
Atl_E_bound = -14.
lon_bnds = [Atl_W_bound,Atl_E_bound]

# use llc_grid_idx_along_lat and data_along_lat functions separately (runs faster)

idx_along_lat,XC_transect,XG_transect = llc_grid_idx_along_lat(lat_transect,lon_bnds)

mask_transect = data_along_lat(land_mask,idx_along_lat,XC_transect,XG_transect)
tw_LHS_1_transect = data_along_lat(therm_wind_LHS_1,idx_along_lat,XC_transect,XG_transect)
tw_LHS_2_transect = data_along_lat(therm_wind_LHS_2,idx_along_lat,XC_transect,XG_transect)
tw_RHS_1_transect = data_along_lat(therm_wind_RHS_1,idx_along_lat,XC_transect,XG_transect)
tw_RHS_2_transect = data_along_lat(therm_wind_RHS_2,idx_along_lat,XC_transect,XG_transect)

[14]:

# make figure
plt.rcParams["font.size"] = 12   # change font size
fig,axs = plt.subplots(4,2,figsize=(10,15))
depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     tw_LHS_1_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[:2,0])
axs[0,0].set_title('dv/dz [s$^{-1}$] along 26 N latitude,\n in lon-lat axes')

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     tw_RHS_1_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[:2,1])
axs[0,1].set_title('$-(g/f\\rho)(d\\rho/dx)$ [s$^{-1}$] along 26 N latitude,\n in lon-lat axes')
plot_objs = [axs[0,0].get_children()[0],axs[1,0].get_children()[0],\
             axs[0,1].get_children()[0],axs[1,1].get_children()[0]]
cmap_zerocent_scale_multiplots(plot_objs,0.2)

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     tw_LHS_2_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[2:,0])
axs[2,0].set_title('du/dz [s$^{-1}$] along 26 N latitude')
axs[3,0].set_xlabel('Longitude')

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     tw_RHS_2_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[2:,1])
axs[2,1].set_title('$(g/f\\rho)(d\\rho/dy)$ [s$^{-1}$] along 26 N latitude')
axs[3,1].set_xlabel('Longitude')
plot_objs = [axs[2,0].get_children()[0],axs[3,0].get_children()[0],\
             axs[2,1].get_children()[0],axs[3,1].get_children()[0]]
cmap_zerocent_scale_multiplots(plot_objs,0.2)

# remove colorbars from left column plots (to declutter figure)
fig.get_children()[-4].remove()
fig.get_children()[-2].remove()

plt.show()

Note that the plots in the left and right columns are very similar at depths below ~100 meters, and very different in the top 100 meters. Why do you think this is?

Velocity reconstruction using thermal wind¶

In early studies of the atmosphere and ocean, measuring in-situ velocities accurately was nearly impossible, but properties such as temperature, humidity, and salinity (that are closely related to density) could be more readily observed. Even in present-day observing systems such as atmospheric radiosondes and oceanic Argo floats, observing the density structure can tell us a lot about the movement of the fluid, because of the application of thermal wind balance.

Given the thermal wind balance equations introduced earlier, integrating vertically between two depth levels $z = z_0$ and $z = z_1$ gives us

$v_g(z_1) - v_g(z_0) = -\frac{g}{f}\int_{z_0}^{z_1} \frac{1}{\rho}\frac{\partial{\rho}}{\partial{x}}\, dz$

$u_g(z_1) - u_g(z_0) = \frac{g}{f}\int_{z_0}^{z_1} \frac{1}{\rho}\frac{\partial{\rho}}{\partial{y}}\, dz$

Since velocities in the deep ocean tend to be much smaller than those near the surface, a common technique used in oceanography has been to assume a level of no motion, i.e., let $u_g(z_0) = v_g(z_0) = 0$ . Then upper ocean geostrophic velocities can be inferred from vertical profiles of density, spaced at regular intervals. Here we can assess the accuracy of this technique using ECCO v4r4 output.

Assessment¶

Thermal wind-based velocity reconstruction¶

Let’s integrate geostrophic velocity from a depth level where we assume the velocity is zero (even though we know from ECCO output that the velocity will not be exactly zero). In this example we’ll take the level of no motion to be 3000 meters depth. Then we consider the result of this velocity reconstruction along the Atlantic 26 N transect we just considered.

[15]:

# set level of no motion
z_0 = -3000.

# integrate (sum) RHS of thermal wind balance upward
dist_z0_upper = ds_grid.Zl - z_0
dist_z0_upper[ds_grid.Zl < z_0] = 0
delta_upper = (-dist_z0_upper.diff("k_l")).values
delta_upper = np.concatenate((delta_upper,np.array([0,])))
delta_upper = np.expand_dims(delta_upper,axis=(0,2,3,4))
v_g_reconstr_upward = (((delta_upper*therm_wind_RHS_1)\
                       .sel(k=slice(None,None,-1)))\
                        .cumsum(dim="k",skipna=True))\
                        .sel(k=slice(None,None,-1))    # slice object "flips" array along k dimension
u_g_reconstr_upward = (((delta_upper*therm_wind_RHS_2)\
                       .sel(k=slice(None,None,-1)))\
                        .cumsum(dim="k",skipna=True))\
                        .sel(k=slice(None,None,-1))    # slice object "flips" array along k dimension
# interpolate DataArrays to center of tracer cells
v_g_reconstr_upward_interp = v_g_reconstr_upward
v_g_reconstr_upward_interp.isel(k=np.arange(0,len(ds_grid.k) - 1)).values = \
                            v_g_reconstr_upward.isel(k=np.arange(0,len(ds_grid.k) - 1))\
                            + ((v_g_reconstr_upward.diff("k").values)/2)
u_g_reconstr_upward_interp = u_g_reconstr_upward
u_g_reconstr_upward_interp.isel(k=np.arange(0,len(ds_grid.k) - 1)).values = \
                            u_g_reconstr_upward.isel(k=np.arange(0,len(ds_grid.k) - 1))\
                            + ((u_g_reconstr_upward.diff("k").values)/2)

# integrate (sum) RHS of thermal wind balance downward
dist_z0_lower = ds_grid.Zu - z_0
dist_z0_lower[ds_grid.Zu > z_0] = 0
delta_lower = (-dist_z0_lower.diff("k_u")).values
delta_lower = np.concatenate((np.array([0,]),delta_lower))
delta_lower = np.expand_dims(delta_lower,axis=(0,2,3,4))
v_g_reconstr_downward = (delta_lower*therm_wind_RHS_1)\
                       .cumsum(dim="k",skipna=True)
u_g_reconstr_downward = (delta_lower*therm_wind_RHS_2)\
                        .cumsum(dim="k",skipna=True)
v_g_reconstr_downward_interp = v_g_reconstr_downward.isel(k=np.arange(0,len(ds_grid.k) - 1))\
                                + ((v_g_reconstr_downward.diff("k").values)/2)
v_g_reconstr_downward_interp = v_g_reconstr_downward_interp.pad(k=(1,0),constant_value=0)
u_g_reconstr_downward_interp = u_g_reconstr_downward.isel(k=np.arange(0,len(ds_grid.k) - 1))\
                                + ((u_g_reconstr_downward.diff("k").values)/2)
u_g_reconstr_downward_interp = u_g_reconstr_downward_interp.pad(k=(1,0),constant_value=0)

# combine results from upward and downward integration
v_g_reconstr = v_g_reconstr_upward_interp
v_g_reconstr.values[:,v_g_reconstr.Z.values < z_0,:,:,:] = \
                                v_g_reconstr_downward_interp.values[:,v_g_reconstr.Z.values < z_0,:,:,:]
u_g_reconstr = u_g_reconstr_upward_interp
u_g_reconstr.values[:,u_g_reconstr.Z.values < z_0,:,:,:] = \
                                u_g_reconstr_downward_interp.values[:,u_g_reconstr.Z.values < z_0,:,:,:]

# rotate (actual) velocities from model to geographical axes
u_zonal = (CS*u) - (SN*v)
v_merid = (SN*u) + (CS*v)

[16]:

# generate DataArrays along line of latitude
lat_transect = 26.
Atl_W_bound = -82.
Atl_E_bound = -14.
lon_bnds = [Atl_W_bound,Atl_E_bound]

# # only need this line uncommented if new transect is specified above
# idx_along_lat,XC_transect,XG_transect = llc_grid_idx_along_lat(lat_transect,lon_bnds)

mask_transect = data_along_lat(land_mask,idx_along_lat,XC_transect,XG_transect)
u_zonal_transect = data_along_lat(u_zonal,idx_along_lat,XC_transect,XG_transect)
v_merid_transect = data_along_lat(v_merid,idx_along_lat,XC_transect,XG_transect)
u_g_reconstr_transect = data_along_lat(u_g_reconstr,idx_along_lat,XC_transect,XG_transect)
v_g_reconstr_transect = data_along_lat(v_g_reconstr,idx_along_lat,XC_transect,XG_transect)

# make figure
plt.rcParams["font.size"] = 12   # change font size
fig,axs = plt.subplots(4,2,figsize=(10,15))
depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     v_merid_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[:2,0])
axs[0,0].set_title('Meridional v [m s$^{-1}$]\n along 26 N latitude')

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     v_g_reconstr_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[:2,1])
axs[0,1].set_title('v reconstruction [m s$^{-1}$]\n from z_0 = ' + str(int(-z_0)) + ' meters depth')
plot_objs = [axs[0,0].get_children()[0],axs[1,0].get_children()[0],\
             axs[0,1].get_children()[0],axs[1,1].get_children()[0]]
cmap_zerocent_scale_multiplots(plot_objs,0.5)

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     u_zonal_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[2:,0])
axs[2,0].set_title('Zonal u [m s$^{-1}$]')
axs[3,0].set_xlabel('Longitude')

depth_two_subplots(XG_transect,-ds_grid.Zl,\
                     u_g_reconstr_transect.squeeze(),\
                     k_split=25,cmap='seismic',mask=mask_transect,\
                      fig=fig,axs=axs[2:,1])
axs[2,1].set_title('u reconstruction [m s$^{-1}$]')
axs[3,1].set_xlabel('Longitude')
plot_objs = [axs[2,0].get_children()[0],axs[3,0].get_children()[0],\
             axs[2,1].get_children()[0],axs[3,1].get_children()[0]]
cmap_zerocent_scale_multiplots(plot_objs,0.5)

# remove colorbars from left column plots (to declutter figure)
fig.get_children()[-4].remove()
fig.get_children()[-2].remove()

plt.show()

Latitude and depth dependence¶

Now as we did with geostrophic balance in the previous tutorial, we can compute the normalized difference over the global domain, sorted into bins by latitude and depth. Note: the ``mean_weighted_binned`` function may take a few minutes to run over the global domain.

[17]:

u_diff = u_zonal - u_g_reconstr
v_diff = v_merid - v_g_reconstr
vel_diff_complex = u_diff + (1j*v_diff)    # in Python, imaginary number i is indicated by 1j
vel_complex = u_zonal + (1j*v_merid)

# normalize magnitude of difference vector by magnitude of actual velocity
vel_diff_abs = np.abs(vel_diff_complex)
vel_abs = np.abs(vel_complex)
vel_diff_norm = vel_diff_abs/vel_abs

# bins of latitude
lat_bin_spacing = 2
lat_bin_bounds = np.c_[-90:90:lat_bin_spacing] + np.array([[0,lat_bin_spacing]])
lat_bin_centers = np.mean(lat_bin_bounds,axis=-1)

# bins of depth
depth_bin_bounds = -ds_grid.Z_bnds.values


# apply "small velocity" mask as in geostrophic balance tutorial
mask_threshold = 0.005    # mask out velocities <0.5 cm s-1
mask_smallvel = (vel_abs < mask_threshold)
curr_weighting = ((~mask_smallvel)*(ds_grid.maskC)*(ds_grid.rA)).values

# create depth array that will broadcast correctly across horizontal dimensions
depth_expand_dims = np.expand_dims(-ds_grid.Z,axis=(-3,-2,-1))

# compute normalized differences
diff_norm_lat = mean_weighted_binned(vel_diff_norm.values,curr_weighting,\
                                    ds_grid.YC.values,lat_bin_bounds)
diff_norm_depth = mean_weighted_binned(vel_diff_norm.values,curr_weighting,\
                                    depth_expand_dims,depth_bin_bounds)

# apply additional masks for depth (100-1000 m) and latitude (>5 deg from equator) ranges

depth_mask = np.logical_and(depth_expand_dims > 100,\
                            depth_expand_dims < 1000)
lat_mask = np.expand_dims((np.abs(ds_grid.YC) > 5),axis=0)


# compute norm diff with additional masks

curr_weighting = (depth_mask*(~mask_smallvel)*(ds_grid.maskC)*(ds_grid.rA)).values
diff_norm_lat_moremask = mean_weighted_binned(vel_diff_norm.values,curr_weighting,\
                                                ds_grid.YC.values,lat_bin_bounds)

curr_weighting = (lat_mask*(~mask_smallvel)*(ds_grid.maskC)*(ds_grid.rA)).values
diff_norm_depth_moremask = mean_weighted_binned(vel_diff_norm.values,curr_weighting,\
                                                depth_expand_dims,depth_bin_bounds)

[18]:

# load norm diff from previous geostrophic balance tutorial
# (if you don't have this file run that tutorial,
# or else comment out the lines containing geostr_file below)
geostr_file = np.load('diff_norm_geostr_bal.npz')

# plot normalized RMS difference, alongside values for geostrophic balance
fig,axs = plt.subplots(2,2,figsize=(15,14))
curr_ax = axs[0,0]
curr_ax.plot(lat_bin_centers,geostr_file['diff_norm_lat'],color='blue',label='Geostrophy')
curr_ax.plot(lat_bin_centers,diff_norm_lat,color='red',label='Thermal wind')
curr_ax.set_ylim(0,1)
curr_ax.axhline(y=0.1,color='black',lw=2)
curr_ax.axvline(x=0,color='black',lw=1)
curr_ax.set_xlabel('Latitude')
curr_ax.set_ylabel('Normalized diff')
curr_ax.set_title('Normalized diff variation with latitude\n z_0 = ' + str(int(-z_0)) + ' meters depth')
curr_ax.legend()
curr_ax = axs[0,1]
curr_ax.plot(geostr_file['diff_norm_depth'],-ds_grid.Z,color='blue',label='Geostrophy')
curr_ax.plot(diff_norm_depth,-ds_grid.Z,color='red',label='Thermal wind')
# Note: in older versions of Matplotlib
# may need to use linthreshy instead of linthresh
curr_ax.set_yscale('symlog',linthresh=100)
curr_ax.invert_yaxis()
curr_ax.axvline(x=0.1,color='black',lw=2)
curr_ax.axhline(y=-z_0,color='black',lw=1)
curr_ax.text(0.02,-0.92*z_0,'z_0')
curr_ax.set_yticks([0,50,100,500,1000,5000])
curr_ax.set_yticklabels(['0','50','100','500','1000','5000'])
curr_ax.set_xlim(0,1)
curr_ax.set_xlabel('Normalized diff')
curr_ax.set_ylabel('Depth [m]')
curr_ax.set_title('Normalized diff variation with depth\n z_0 = ' + str(int(-z_0)) + ' meters depth')
curr_ax.legend()
curr_ax = axs[1,0]
curr_ax.plot(lat_bin_centers,geostr_file['diff_norm_lat_moremask'],color='blue',label='Geostrophy')
curr_ax.plot(lat_bin_centers,diff_norm_lat_moremask,color='red',label='Thermal wind')
curr_ax.set_ylim(0,1)
curr_ax.axhline(y=0.1,color='black',lw=2)
curr_ax.axvline(x=0,color='black',lw=1)
curr_ax.set_xlabel('Latitude')
curr_ax.set_ylabel('Normalized diff')
curr_ax.set_title('In 100-1000 m depth range only')
curr_ax.legend()
curr_ax = axs[1,1]
curr_ax.plot(geostr_file['diff_norm_depth_moremask'],-ds_grid.Z,color='blue',label='Geostrophy')
curr_ax.plot(diff_norm_depth_moremask,-ds_grid.Z,color='red',label='Thermal wind')
# Note: in older versions of Matplotlib
# may need to use linthreshy instead of linthresh
curr_ax.set_yscale('symlog',linthresh=100)
curr_ax.invert_yaxis()
curr_ax.axvline(x=0.1,color='black',lw=2)
curr_ax.axhline(y=-z_0,color='black',lw=1)
curr_ax.text(0.02,-0.92*z_0,'z_0')
curr_ax.set_yticks([0,50,100,500,1000,5000])
curr_ax.set_yticklabels(['0','50','100','500','1000','5000'])
curr_ax.set_xlim(0,1)
curr_ax.set_xlabel('Normalized diff')
curr_ax.set_ylabel('Depth [m]')
curr_ax.set_title('>5 degrees from equator only')
curr_ax.legend()

plt.show()

Comparing the geostrophy and thermal wind normalized differences provides some insight into the challenges of reconstructing velocity fields using thermal wind. Geostrophic balance is a more accurate predictor of velocity than thermal wind, since thermal wind assumes:

geostrophic balance
hydrostatic balance (not a factor here since hydrostatic balance is actually assumed in the model code itself)
a level of no motion (a potentially more problematic assumption, especially close to $z_0$ ).

The plots above also illustrate why satellite observations of sea surface height are so valuable to oceanographers. Prior to reliable measurements of sea surface height, ocean velocity fields were typically reconstructed using density profiles and thermal wind. When sea surface height data are combined with density profiles, we can compute subsurface pressure gradients, and then obtain velocity estimates using only geostrophy, which is much more accurate than having to rely on thermal wind.

Exercises¶

Try out different values of $z_0$ (i.e., level of no motion) to use in the velocity reconstructions. How does this affect the normalized differences in the plots above? Which $z_0$ values permit the most accurate reconstructions (lowest normalized difference) at 100 meters depth? (For bonus points, write some code to determine the optimal $z_0$ value that minimizes normalized difference at 100 m depth, rather than using trial and error.)
Adapt the lon_depth_along_lat function (and the two functions it calls: llc_grid_idx_along_lat and data_along_lat) to plot transects following lines of longitude instead of latitude. Use this to look at the application of thermal wind to reconstruct velocities through Drake Passage (~64 deg W, 55-65 deg S).

References¶

Gill, A.E. (1982). Atmosphere-Ocean Dynamics. Academic Press, Elsevier.

Kundu, P.K. and Cohen, I.M. (2008). Fluid Mechanics (4th ed.). Elsevier.

Vallis, G.K. (2006). Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press.