Balanced mass and flow

10 May 2021 by MiniUFO

---

[TOC]

---

1. Introduction

Geopotential model (Yuan et al. 2008), also known as the divergence equation, represents a relation between the mass (geopotential \(\Phi\) in atmospheric context, or sea surface height \(\eta\) in oceanic context) and horizontal flow \((u,v)\) fields. Given the mass field, one can invert the flow field and vice versa. The simplest case is to derive the geostrophic flow given the mass field. Similarly, given the flow field, one can obtain the balanced mass field. This is particularly useful in data assimilation, such as the the technique of normal mode initialization. In this initialization, given a mass (flow) observation, one need to calculate the balanced flow (mass) field so that the high-frequency oscillations are removed. This prevents the model integration from overflow by initial shocks such as fast gravity waves.

In Cartesian coordinates, assuming barotropic flow, applying the horizontal divergence operator \(\nabla_h\cdot()\) to the vector form of the horizontal momentum equation of \(\vec{v_h}\), one obtains the divergence \(D=\nabla_h\cdot\vec{v_h}\) equation as:

\[\begin{align} \frac{\partial D}{\partial t}=-\vec{v_h}\cdot\nabla D - D^2 - \nabla^2 \Phi + f\zeta-u\beta+\nabla\cdot\vec{F}\tag{1} \end{align}\]

The first two terms on the rhs of Eq. (1) is the nonlinear terms, the third to fifth terms are the geostrophic balance terms, and the last term is the divergence of frictional stress. Here we are going to consider a case on the spherical coordinates, possibly include the vertical advection if flow is not barotropic, and keep as many terms as possible for the inversion.

2. Theory

Here the derivation of the geopotential model (or equivalently, the divergence equation) is briefly reviewed. The start point is the horizontal momentum \((u,v)\) equations in spherical-pressure coordinates:

\[\begin{split}\begin{align} \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\omega\frac{\partial u}{\partial p}=-\frac{\partial \Phi}{\partial x}+fv+\frac{uv\tan\phi}{a}+F_{\lambda}\tag{2}\\ \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\omega\frac{\partial v}{\partial p}=-\frac{\partial \Phi}{\partial y}-fu-\frac{u^2\tan\phi}{a}+F_{\phi} \tag{3} \end{align}\end{split}\]

where \(\partial x=a\cos\phi\partial \lambda\) and \(\partial y=a\partial\phi\) are the short-hands. Taking \(\partial/\partial x\) and \(\partial \cos\phi() /\partial y/\cos\phi\) respectively of Eqs. (2) and (3), and suming up the results gives the divergence (\(D\)) equation:

\[\begin{split}\begin{align} \frac{\partial D}{\partial t}=&-\vec{u}\cdot\nabla D-D^2-\frac{\partial}{\partial x}\left(\omega\frac{\partial u}{\partial p}\right)-\frac{1}{\cos \phi}\frac{\partial}{\partial y}\left(\omega\cos\phi\frac{\partial v}{\partial y}\right)\notag\\&-\nabla^2\Phi+f\zeta-u\beta\notag\\ &+\frac{\partial }{\partial x}\left(\frac{uv\tan\phi}{a}\right)-\frac{\partial}{\partial y}\left(\frac{u^2\sin\phi}{a}\right)+\nabla \cdot\vec{F}\tag{4} \end{align}\end{split}\]

2.1 Geopotential model

If the geopotential \(\Phi\) is taken as the unknown and the 3D flow \((u,v,\omega)\) as the knowns, then Eq. (4) is a Poisson equation with respect to \(\Phi\) as:

\[\begin{split}\begin{align} \nabla^2\Phi=&-\frac{\partial D}{\partial t}-\vec{u}\cdot\nabla D-D^2-\frac{\partial}{\partial x}\left(\omega\frac{\partial u}{\partial p}\right)-\frac{1}{\cos \phi}\frac{\partial}{\partial y}\left(\omega\cos\phi\frac{\partial v}{\partial y}\right)\notag\\&+f\zeta-u\beta\notag\\ &+\frac{\partial }{\partial x}\left(\frac{uv\tan\phi}{a}\right)-\frac{1}{\cos \phi}\frac{\partial}{\partial y}\left(\frac{u^2\sin\phi}{a}\right)+\nabla \cdot\vec{F}\tag{5} \end{align}\end{split}\]

One needs to calculate the forcing terms on the rhs of Eq. (4) using 3D flow \((u,v,\omega)\) data and then invert for geopotential \(\Phi\).

2.2 Geostrophic streamfunction model

If the flow fields are unknown and the mass field can be observed (e.g., in oceanic context, the sea surface height can be directly observed by altimetry satellites), Eq. (5) has to be simplied somehow for inversion. First, one need to remove nonlinear terms and the frictional term, which is also known as linearization. Then introducing the streamfunction \(\psi\) so that the horizontal flow is nondivergent. This results in the geostrophic balance as:

\[\begin{align} \nabla \cdot(f\nabla\psi)=\nabla^2\Phi\tag{6} \end{align}\]

Similarly, given the mass field, one can invert the geostrophic streamfunction iteratively. However, near the equator, \(f\rightarrow0\) and Eq. (6) is not elliptic at the equator.

3. Examples

3.1 The geopotential model example

Now we have a global horizontal wind field \((u,v)\), and we want to recover the geopotential field \(\Phi\). All the terms on the rhs of Eq. (5) will be calculated and then inverted using SOR accordingly for each responses. The code is as follows.

### geopotential model example ###
import sys
sys.path.append('../../../')
import xarray as xr

# loading data from a dataset
ds = xr.open_dataset('../../../Data/atmos3D.nc').sel({'LEV':500})

u = ds.U
v = ds.V
h = ds.hgt

Calculate each forcing on the rhs of Eq. (5).

### calculate each forcing ###
import numpy as np
from xinvert import FiniteDiff

Rearth = 6371200.0  # Earth radius
omega = 7.292e-5    # angular speed of the earth rotation

fd = FiniteDiff(dim_mapping={'Y':'lat', 'X':'lon'},
                BCs={'Y':'reflect', 'X':'periodic'}, coords='lat-lon')

vort = fd.curl(u, v)
cosF = np.cos(np.deg2rad(u.lat))
sinF = np.sin(np.deg2rad(u.lat))
beta = 2.0 * omega * cosF / Rearth
fcor = 2.0 * omega * sinF


# force 2
tmp2 = fd.grad(u, 'X', BCs='periodic') * u
frc2 = -fd.grad(tmp2, 'X').load()
frc2 = xr.where(np.isfinite(frc2), frc2, np.nan)

# force 3
tmp3 = fd.grad(v, 'X') * u * cosF
frc3 = -fd.grad(tmp3, 'Y', BCs='extend').load() / cosF
frc3 = xr.where(np.isfinite(frc3), frc3, np.nan)

# force 4
tmp4 = fd.grad(u, 'Y', BCs='extend') * v
frc4 = -fd.grad(tmp4, 'Y').load()
frc4 = xr.where(np.isfinite(frc4), frc4, np.nan)

# force 5
tmp5 = fd.grad(v, 'Y', BCs='extend') * v*cosF
frc5 = -fd.grad(tmp5, 'Y', BCs='extend').load() / cosF
frc5 = xr.where(np.isfinite(frc5), frc5, np.nan)

# force 9
frc9 = (fcor * vort).load()
frc9 = xr.where(np.isfinite(frc9), frc9, np.nan)

# force 10
frc10 = -(u * beta).load()
frc10 = xr.where(np.isfinite(frc10), frc10, np.nan)

# force 11
tmp11 = -u**2 * sinF / Rearth
frc11 = fd.grad(tmp11, 'Y', BCs='extend').load() / cosF
frc11 = xr.where(np.isfinite(frc11), frc11, np.nan)

# force 12
tmp12 = u * v * sinF / cosF / Rearth
frc12 = fd.grad(tmp12, 'X').load()
frc12 = xr.where(np.isfinite(frc12), frc12, np.nan)


# all forces
frcAll = frc2 + frc3 + frc4 + frc5 + frc9 + frc10 + frc11 + frc12
frcAll = xr.where(np.isfinite(frcAll), frcAll, np.nan)
frcAll = frcAll.fillna(0)

Inverting geopotential response to each forcing.

### inverting ###
from xinvert import invert_Poisson

hbc = (h*9.81).load() # boundary values

iParams = {
    'BCs'      : ['fixed', 'periodic'],
    'mxLoop'   : 20000,
    'tolerance': 1e-12,
}

sf = invert_Poisson(frcAll, dims=['lat','lon'], icbc=hbc, iParams=iParams)

{} loops 20000 and tolerance is 5.473987e-09

Now we can plot the forcings and responses.

### plot the sum-ups, observed, and difference ###
import proplot as pplt

array = [
    [1, 1, 2, 2],
    [0, 3, 3, 0],
]

fig, axes = pplt.subplots(array, figsize=(7, 3.8), sharex=3, sharey=3,
                          proj=pplt.Proj('cyl', lon_0=180))

fontsize = 13

ax = axes[0]
p=ax.contourf(h, cmap='jet', levels=np.linspace(4900, 5900,21))
ax.set_title('reanalysis', fontsize=fontsize)
ax.set_ylim([-90, 90])
ax.set_xlim([-180, 180])

ax = axes[1]
p=ax.contourf(sf/9.81, cmap='jet', levels=np.linspace(4900, 5900,21))
ax.set_title('all forcings + boundary', fontsize=fontsize)
ax.set_ylim([-90, 90])
ax.set_xlim([-180, 180])
ax.colorbar(p, loc='r', length=1, label='', ticks=150)

ax = axes[2]
p=ax.contourf(sf/9.81-h, cmap='jet', levels=np.linspace(-30, 30,30), extend='both')
ax.set_title('difference', fontsize=fontsize)
ax.set_ylim([-90, 90])
ax.set_xlim([-180, 180])
ax.colorbar(p, loc='r', length=1, label='', ticks=10)

axes.format(abc='(a)', coast=True,
            lonlines=60, latlines=30, lonlabels='b', latlabels='l',
            grid=False, labels=False)

C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "
C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "
C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "

From Fig. (a) and (b) one can see clearly that the inverted and reanalysis are almost the same. Differences may be caused by neglecting the terms related to vertical motion and temporal change of divergence.

3.2 The geostrophic flow example

In contrast to the above example, if we know the mass field \(\Phi\), we can calculate the geostrophic flow field accordingly. In order to prevent the \(f=0\) at the equator, we shift the grid slightly so that it locates on both sides of the equator.

from xinvert import FiniteDiff

# load data
Phi = (h*9.81).load()

fd = FiniteDiff(dim_mapping={'Z':'LEV', 'Y':'lat', 'X':'lon'},
                BCs={'Z':['extend', 'extend'],
                     'Y':['extend', 'extend'],
                     'X':['periodic','periodic']},
                coords='lat-lon', fill=0)

# calculate rhs forcing of Eq. (6)
force = fd.Laplacian(Phi, dims=['Y', 'X'])

latNew = (h.lat + h.lat.diff('lat')[0]/2.0)[:-1]

ny, nx = force.shape

# shift half latitude so that no data is located at equator (avoiding f=0 exactly)
forceHalf = force.interp_like(xr.DataArray(np.zeros((ny-1, nx)), dims=['lat','lon'], coords={'lat':latNew, 'lon':h.lon})).load()

Inverting for geostrophic streamfunction \(\Phi\):

#%% invert
from xinvert import invert_geostrophic

iParams = {
    'BCs'      : ['fixed', 'periodic'],
    'mxLoop'   : 5000,
    'tolerance': 1e-12,
}

sf = invert_geostrophic(forceHalf, dims=['lat','lon'], iParams=iParams)

{} loops 5000 and tolerance is 8.480886e-10

Results can be plotted as:

#%% plot wind and streamfunction
import proplot as pplt

m = np.hypot(u, v)

lat, lon = xr.broadcast(u.lat, u.lon)
fig, axes = pplt.subplots(nrows=2, ncols=1, figsize=(7, 7), sharex=3, sharey=3,
                          proj=pplt.Proj('cyl', lon_0=180))
fontsize = 13

axes.format(abc='(a)', coast=True,
            lonlines=60, latlines=30, lonlabels='b', latlabels='l',
            grid=False, labels=False)

ax = axes[0]
p = ax.contourf(force*1e7, cmap='seismic',
                levels=np.linspace(-1, 1, 21))
ax.set_title('Laplacian of geopotential', fontsize=fontsize)
ax.colorbar(p, loc='r', label='', ticks=0.2, length=1)

ax = axes[1]
p = ax.contourf(sf, levels=31, cmap='jet')
ax.quiver(lon.values, lat.values, u.values, v.values,
              width=0.0007, headwidth=12., headlength=15.)
              # headwidth=1, headlength=3, width=0.002)
ax.set_title('wind field and inverted streamfunction', fontsize=fontsize)
ax.colorbar(p, loc='r', label='', length=1)

C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "
C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "

<matplotlib.colorbar.Colorbar at 0x28c2670d7b0>

Note that there is a jump of the streamfunction \(\Phi\) at the equator. If we add a constant to the Southern Hemisphere to remove the jump, then the streamfunction field will look much better:

#%% plot wind and streamfunction
import proplot as pplt

sfm = sf.mean('lon')
sf2 = sf.copy()

sf2.loc[:0] += ((sfm[72]-sfm[71])/2.0 - sfm[71]).values

fig, axes = pplt.subplots(nrows=1, ncols=1, figsize=(7, 4.4), sharex=3, sharey=3,
                          proj=pplt.Proj('cyl', lon_0=180))
fontsize = 14

axes.format(abc='(a)', coast=True,
            lonlines=60, latlines=30, lonlabels='b', latlabels='l',
            grid=False, labels=False)

ax = axes[0]
p = ax.contourf(sf2, levels=21, cmap='jet')
ax.set_ylim([-80, 80])
ax.set_xlim([-180, 175])
ax.quiver(lon.values[::2,::2], lat.values[::2,::2], u.values[::2,::2], v.values[::2,::2],
              width=0.0007, headwidth=12., headlength=15.)
              # headwidth=1, headlength=3, width=0.002)
ax.set_title('wind field and inverted streamfunction', fontsize=fontsize)
ax.colorbar(p, loc='b', label='', length=1)

C:\ProgramData\anaconda3\lib\site-packages\cartopy\mpl\geoaxes.py:406: UserWarning: The `map_projection` keyword argument is deprecated, use `projection` to instantiate a GeoAxes instead.
  warnings.warn("The `map_projection` keyword argument is "

<matplotlib.colorbar.Colorbar at 0x28c2b942560>

This is the observed wind field and is aligned with the streamfunction, especially away from the tropics. Even in the equatorial regions, the wind is along the streamfunction.

---

References

1. Yuan Z., J. Wu, X. Cheng, and M. Jian 2008: The derivation of a numerical diagnostic model for the forcing of the geopotential. Q. J. R. Meteorol. Soc., 134, 2067-2078.