.ipynb .pdf

repository open issue

Binder

Contents

Ocean Barotropic Volume Transport

In this lecture, we will use the ECCO State Estimate to explore and describe the ocean circulation. As discussed in Ocean State Estimation, the flow field produced by ECCO is constrained by nearly all relevant and available observations of the ocean and has been carefully validated across multiple dimensions. Nevertheless, we must remember that ECCO is not the real ocean. A particular shortcoming of this solution is its relatively low spatial resolution of approx. 1 degree. This means it does not resolve sharp boundary currents, fronts, or mesoscale variability.

Barotropic (Horizontal) Circulation

As a starting point for our study of ocean circulation, we will examine the Barotropic Circulation. Formally, a barotropic fluid is one where density depends only on pressure. This would be the case for the ocean only if temperature and salinity were uniform (or at least dependent only on depth). The real ocean is not a barotropic fluid; it is baroclinic (the opposite of barotropic).

However, in common usage in physical oceanography, Barotropic is taken to mean “depth-averaged”. By averaging a proprerty (usually the flow field) over depth, we obtain a simplified picture of the flow. As we will see in Vorticity, Sverdrup Transport, and Ocean Gyres, classical ocean circulation theory makes predicts for the barotropic flow.

Defining the Barotropic Flow

The full 3D velociety field of the ocean is the vector field \(\mathbf{u}(\mathbf{x}, t) = (u, v, w)\). Visualizing and understanding a time-varying vector field is very challenging, so we often want to remove one dimension.

As a first step, we define vertically-integrated flow components as

(29)\[\begin{align} U &= \int_{-H}^\eta u dz \\ V &= \int_{-H}^\eta v dz \\ \end{align}\]

Continuity Equation

Now we apply the same integral to the Boussinesq continuity equation: \(\nabla \cdot \mathbf{u} = 0\) to obtain

(30)\[\begin{align} \int_{-H}^\eta \frac{\partial u}{\partial x} dz+ \int_{-H}^\eta \frac{\partial v}{\partial y} dz + \int_{-H}^\eta \frac{\partial w}{\partial z} dz &= 0 \\ \frac{\partial U}{\partial x} + u|_{z=\eta} \frac{\partial \eta}{\partial x} +\frac{\partial V}{\partial y} + v|_{z=\eta} \frac{\partial \eta}{\partial x} +w|_{z=\eta} - w|_{z=-H} &= 0 \end{align}\]

The surface boundary condition is:

\[ \frac{D \eta}{D t} = \frac{\partial \eta}{\partial t} + u|_{z=\eta} \frac{\partial \eta}{\partial x} + v|_{z=\eta} \frac{\partial \eta}{\partial y} = w|_{z=\eta} - (E - P - R) \ . \]

We will assume that the bottom is rigid and flat, in which case the bottom boundary condition is just \(w|_{z=-H} = 0\).

With these boundary conditions, we obtain.

\[ \frac{\partial \eta}{\partial t} = - \frac{\partial U}{\partial x} - \frac{\partial V}{\partial y} - (E - P - R) \]

This equation describes how the sea surface height (\(\eta\)) evolves in time in a Boussinesq, barotropic ocean. It can change either due to a convergence of horizontal volume fluxes or due to a flux of water from the atmosphere.

“Rigid Lid” and Streamfunction

A common and appropriate approximation is to assume that \(\partial \eta / \partial t\) and \(E - P - R\) are small compared to the terms related to the horizontal velocities. With this, we obtain

\[ \frac{\partial U}{\partial x} + \frac{\partial V}{\partial y} = 0 \ . \]

This is call two-dimensional non-divergent flow. We can always represent such flow using a streamfunction \(\Psi\), defined as

\[ U = -\frac{\partial \Psi}{\partial y} \ , \ \ V = \frac{\partial \Psi}{\partial x} \]

Any velocity field that is defined in this way will satisfy the 2D continuity equation by construction.

\[ \frac{\partial U}{\partial x} + \frac{\partial V}{\partial y} = -\frac{\partial^2 \Psi}{\partial x \partial y} + \frac{\partial^2 \Psi}{\partial x \partial y} = 0 \ . \]

Streamfunctions are great because they summarize the vector flow field in a single scalar field. \(\Psi\) as defined here has units of m\(^3 /\) s. Horizontal boundary conditions are auomatically satisfied by setting \(\Psi = \)const. on the boundaries.

import xgcm
import xarray as xr
from matplotlib import pyplot as plt
import hvplot.xarray
import numpy as np

plt.rcParams['figure.figsize'] = (12, 8)

from dask_gateway import Gateway
gway = Gateway()
cluster = gway.new_cluster()
cluster

/srv/conda/envs/notebook/lib/python3.8/site-packages/dask_gateway/client.py:21: FutureWarning: format_bytes is deprecated and will be removed in a future release. Please use dask.utils.format_bytes instead.
  from distributed.utils import LoopRunner, format_bytes

from ecco_utils import load_ecco_v4
ds = load_ecco_v4()
ds

<xarray.Dataset>
Dimensions:    (face: 13, lon_c: 360, lon_g: 360, lat_c: 240, lat_g: 240, Z: 50, Zl: 50, Zp1: 51, Zu: 50, time: 288, time_snp: 287)
Coordinates: (12/42)
  * face       (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
    i          (lon_c) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    i_g        (lon_g) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    j          (lat_c) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    j_g        (lat_g) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    k          (Z) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    dyG        (lat_c, lon_g) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
    dxG        (lat_g, lon_c) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
  * lon_c      (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lon_g      (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c      (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lat_g      (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables: (12/35)
    ADVr_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ADVr_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ...         ...
    UVELMASS   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELMASS   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    UVELSTAR   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELSTAR   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    oceTAUX    (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    oceTAUY    (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>

try:
    ds_bt = xr.open_dataset("tmp_data/ECCO_barotropic_velocity.nc")
    U_bt = ds_bt.U_bt
    V_bt = ds_bt.V_bt
except FileNotFoundError:
    U_bt = (ds.UVELMASS * ds.drF).sum('Z')
    V_bt = (ds.VVELMASS * ds.drF).sum('Z')
    cluster.scale(10)
    with cluster.get_client():
        U_bt.load()
        V_bt.load()    
    cluster.scale(0)
    ds_tmp = xr.merge([U_bt.rename('U_bt'), V_bt.rename('V_bt')])
    ds_tmp.to_netcdf("tmp_data/ECCO_barotropic_velocity.nc")
U_bt

<xarray.DataArray 'U_bt' (time: 288, lat_c: 240, lon_g: 360)>
[24883200 values with dtype=float32]
Coordinates:
    i_g      (lon_g) int64 0 1 2 3 4 5 6 7 8 ... 352 353 354 355 356 357 358 359
    j        (lat_c) int64 30 31 32 33 34 35 36 ... 263 264 265 266 267 268 269
  * time     (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
    iter     (time) int64 732 1428 2172 2892 ... 208164 208908 209628 210360
    rAw      (lat_c, lon_g) float32 7.659e+08 7.671e+08 ... 3.766e+08 1.931e+08
    dxC      (lat_c, lon_g) float32 2.072e+04 2.071e+04 ... 2.25e+04 1.812e+04
    dyG      (lat_c, lon_g) float32 3.697e+04 3.704e+04 ... 1.971e+04 1.597e+04
  * lon_g    (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 318.0 319.0 320.0 321.0
  * lat_c    (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47

grid = xgcm.Grid(ds, periodic=['X'])

rA = ds.rA.compute()
vol_flux_convergence_U = - (grid.diff(U_bt * U_bt.dyG, 'X', boundary='extend') / rA).fillna(0)
vol_flux_convergence_V = - (grid.diff(V_bt * V_bt.dxG, 'Y', boundary='extend') / rA).fillna(0)
vol_flux_convergence = vol_flux_convergence_V + vol_flux_convergence_U
rho0 = 1029.
vol_flux_surface = (ds.oceFWflx.compute()/rho0)

big_plot = (vol_flux_convergence_U.hvplot.quadmesh('lon_c', 'lat_c', clim=(-1e-4, 1e-4), cmap='RdBu_r', label='Zonal_Flux_Convergence') +
 vol_flux_convergence_V.hvplot.quadmesh('lon_c', 'lat_c', clim=(-1e-4, 1e-4), cmap='RdBu_r', label='Meridional_Flux_Convergence') +
 vol_flux_convergence.hvplot.quadmesh('lon_c', 'lat_c', clim=(-1e-7, 1e-7), cmap='RdBu_r', label='Total_Lateral_Flux_Convergence') +
 vol_flux_surface.hvplot.quadmesh('lon_c', 'lat_c', clim=(-1e-7, 1e-7), cmap='RdBu_r', label='Surface_Flux')
).opts(shared_axes=False, width=400).cols(2)

Testing “Rigid Lid” Approximation

Below we examine the two ocean transport terms in depth-integrated continuity equation individually, and then summed together and compared with the surface volume flux.

This demonstrates the high degree of cancelation between the transport terms; their sum, mostly balanced by \(E - P - R\) is more than three order of magnitude smaller than each term individually.

big_plot

Global Barotropic Streamfunction

The plot of the global barotropic circulation is dominated by the ACC.

Psi_bt = grid.cumsum(-U_bt * ds.dyG.compute(), 'Y', boundary='fill')
mask_Z = grid.interp(ds.hFacS[0].compute(), 'X')
Psi_bt = Psi_bt.where(mask_Z) / 1e6
Psi_bt.attrs['long_name'] = 'Barotropic Streamfunction'
Psi_bt.attrs['units'] = 'Sv'

# mv U_bt and V_bt to the Psi grid
U = grid.interp(U_bt, 'Y', boundary='fill').rename('U')
V = grid.interp(V_bt, 'X').rename('V')
UV = xr.merge([U, V])

import cartopy.crs as ccrs

fig = plt.figure(figsize=(14, 6))
ax = plt.axes(projection=ccrs.Robinson(central_longitude=160), facecolor='0.9')
ax.coastlines()

clevels = np.arange(-170, 51, 20)

Psi_bt.mean('time').plot.contourf(levels=clevels, ax=ax, cmap='viridis', transform=ccrs.PlateCarree(), cbar_kwargs={'shrink': 0.5, 'label': 'Sv'})
plt.close()

fig

Antarctic Circumpolar Current

import cartopy.crs as ccrs

fig = plt.figure(figsize=(10, 9))
ax = plt.axes(projection=ccrs.SouthPolarStereo(), facecolor='0.9')
ax.coastlines()

clevels = np.arange(-190, 51, 10)

Psi_to_plot = Psi_bt.sel(lat_g=slice(None, -20)).mean('time')

Psi_to_plot.plot.contourf(levels=clevels, ax=ax, cmap='viridis', transform=ccrs.PlateCarree(), cbar_kwargs={'shrink': 0.5, 'label': 'Sv'})
Psi_to_plot.plot.contour(levels=[-140, 0], ax=ax, colors='white', linewidths=3, transform=ccrs.PlateCarree(), add_colorbar=False)
ax.set_extent((0, 360, -90, -40), crs=ccrs.PlateCarree())

points_to_label = {
    'Weddell Gyre': (-20, -65),
    'Ross Gyre': (-150, -68),
    'Drake Passage': (-85, -55),
    'Kerguelen Plateau': (65, -60)
}


for label, (lon, lat) in points_to_label.items():
    ax.annotate(label, xy=ax.projection.transform_point(lon, lat, src_crs=ccrs.PlateCarree()), color='white')
ax.set_title("Southern Ocean Barotropic Streamfunction")
    
plt.close()

fig

Important Points

• There are two clockwise-spinning Gyres off the Antarctic continent

  • The Weddell Gyre in the Atlantic sector has a strength of about 30 Sv

  • The Ross Gyre in the Pacific sector has a strength of about 10 Sv

• The Antarctic Circumpolar Current (ACC) is a very strong eastward-flowing current that circumnavigates Antarctica

• The ACC has a strong choke-point in Drak Passage and a large split / meander at Kerguelen

Drake Passage Transport

fig, ax = plt.subplots()
Psi_Drake_Passage = Psi_bt.sel(lon_g=290, method='nearest').sel(lat_g=slice(-70, -55))
Psi_Drake_Passage.plot.contourf(x='time', levels=np.arange(-180, 1, 20), cmap='viridis', ax=ax)
plt.close()

fig

fig, ax = plt.subplots()
(-Psi_Drake_Passage[:, -2]).plot(ax=ax)
ax.grid()
ax.set_title('ECCO Drake Passage Transport')
plt.close()

fig

North Pacific Gyre

fig = plt.figure(figsize=(16, 6))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=180), facecolor='0.9')
ax.coastlines()

Psi_npac = Psi_bt.sel(lat_g=slice(10, 62), lon_g=slice(110, 260))
Psi_npac.mean('time').plot.contourf(ax=ax, levels=20, transform=ccrs.PlateCarree())
ax.set_extent((120, 250, 15, 62), crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)

points_to_label = {
    'Kuroshio': (140, 33),
    'Oyashio': (160, 50),
    'Kuroshio Extension / North Pacific Current': (160, 40),
    'Alaska Current': (-150, 53),
    'North Equatorial Current': (150, 20),
    'California Current': (-140, 32),
}

for label, (lon, lat) in points_to_label.items():
    ax.annotate(label, xy=ax.projection.transform_point(lon, lat, src_crs=ccrs.PlateCarree()), color='white')
ax.set_title("North Pacific Barotropic Streamfunction")

plt.close()

fig

Important Points

• Large Subtropical gyre w/ Kuroshio western boundary current

• Smaller Subpolar gyre with Oyashio western boundary current

• North Equatorial Current feeds Kuroshio

Seasonality

The North Pacific Gyre has a strong seasonal cycle.

Psi_npac_mm = Psi_npac.groupby('time.month').mean('time')
Psi_npac_mm_anom = Psi_npac_mm - Psi_npac_mm.mean('month')


subplot_kws = dict(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=190), facecolor='0.9')
#ax.coastlines()

fg = Psi_npac_mm_anom.sel(month=slice(None, None, 3)).plot.contourf(
    levels=20, transform=ccrs.PlateCarree(), col='month', col_wrap=2, subplot_kws=subplot_kws,
    cbar_kwargs={'shrink': 0.5},
    figsize=(18, 8)
)
for ax in fg.axes.ravel():
    ax.set_extent((120, 250, 15, 62), crs=ccrs.PlateCarree())
    ax.coastlines()
fig = plt.gcf()
plt.close()

fig

Indian Ocean

fig = plt.figure(figsize=(16, 6))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=75), facecolor='0.9')
ax.coastlines()


Psi_ind = Psi_bt.sel(lat_g=slice(-40, 30), lon_g=slice(20, 130))
Psi_ind.mean('time').plot.contourf(ax=ax, levels=40, transform=ccrs.PlateCarree())
ax.set_extent((25, 120, -35, 25), crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)
ax.set_title("Indian Ocean Barotropic Streamfunction")

plt.close()

fig

Seasonality

Psi_ind_mm = Psi_ind.groupby('time.month').mean('time')
Psi_ind_mm_anom = Psi_ind_mm - Psi_ind_mm.mean('month')


subplot_kws = dict(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=75), facecolor='0.9')
#ax.coastlines()

fg = Psi_ind_mm_anom.sel(month=slice(None, None, 3)).plot.contourf(
    levels=20, transform=ccrs.PlateCarree(), col='month', col_wrap=2, subplot_kws=subplot_kws,
    cbar_kwargs={'shrink': 0.5},
    figsize=(18, 8)
)
for ax in fg.axes.ravel():
    ax.set_extent((25, 120, -35, 25), crs=ccrs.PlateCarree())
    ax.coastlines()
    #(Psi_ind + 140).mean('time').plot.contour(ax=ax, levels=20, colors='k', transform=ccrs.PlateCarree())

fig = plt.gcf()
plt.close()

fig

Important Points - Indian

• Subtropical Gyre merges with ACC at the basin edge

• Equatorial Jets

• Circulation in Arabian Sea / Bay of Bengal is highly seasonal. (Dominated by Monsoon wind anomalies.)

North Atlantic Ocean

Psi_natl = Psi_bt.roll(lon_g=-50).isel(lon_g=slice(250, None)).sel(lat_g=slice(0, None))
lon_g_fixed = Psi_natl.lon_g.values
lon_g_fixed[lon_g_fixed > 200] -= 360
Psi_natl = Psi_natl.assign_coords({'lon_g': ('lon_g', lon_g_fixed, Psi_natl.lon_g.attrs)})

/srv/conda/envs/notebook/lib/python3.8/site-packages/xarray/core/dataarray.py:3267: FutureWarning: roll_coords will be set to False in the future. Explicitly set roll_coords to silence warning.
  ds = self._to_temp_dataset().roll(

fig = plt.figure(figsize=(16, 6))
ax = plt.axes(projection=ccrs.LambertAzimuthalEqualArea(central_longitude=-40), facecolor='0.9')
ax.coastlines()

Psi_natl.mean('time').plot.contourf(ax=ax, levels=40, transform=ccrs.PlateCarree())
ax.set_extent((-80, 0, 15, 65), crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)

points_to_label = {
    'Gulf Stream': (-79, 30),
    'Gulf Stream Extension / North Atlantic Current': (-60, 41),
    'Subpolar Gyre': (-45, 55),

}

for label, (lon, lat) in points_to_label.items():
    ax.annotate(label, xy=ax.projection.transform_point(lon, lat, src_crs=ccrs.PlateCarree()), color='white')
ax.set_title("North Atlantic Barotropic Streamfunction")

plt.close()

fig

previous

Ocean State Estimation

next

Overturning Circulation