Turning ripples into waves

Water-mass mixing percentage in the South Atlantic correction

@ClarkGRichards, one of oce developers, pointed out two mistakes I made in the original Water Mass Mixing post. One was the fact that I did not use conservative temperature from TEOS-10 (or at least potential temperature from EOS-10) for the T$-$S plot. That resulted in an inversion at the end of the T$-$S diagram. The other mistake was the fact that I used units in the salinity label (g kg$^{-1}$), meaning Absolute Salinity, but instead I plotted Practical Salinity (no units).

I am re-posting the original notebook here fixing those mistakes.

In [2]:
import numpy as np
import as ma

def mixing(T, S, inds):
    Compute the water mass mixing percentage using Mamayev's (1975) mixing

    T : Conservative Temperature
    S : Absolute Salinity
    inds :  2x3 array with thermohaline indices
            [T1 T2 T3
            S1 S2 S3]

    m1, m2, m3 : Water mass percentage for masses 1, 2 e 3.


    a = np.r_[inds, np.ones((1, 3))]
    b = np.c_[T.ravel(), S.ravel(), np.ones(T.shape).ravel()].T
    m = np.linalg.solve(a, b)
    m1 = m[0].reshape(T.shape)
    m2 = m[1].reshape(T.shape)
    m3 = m[2].reshape(T.shape)

    # Mask values outside mixing triangle.
    m1 = ma.masked_outside(ma.masked_invalid(m1), 0, 1)
    m2 = ma.masked_outside(ma.masked_invalid(m2), 0, 1)
    m3 = ma.masked_outside(ma.masked_invalid(m3), 0, 1)

    m1 = 100 * m1
    m2 = 100 * m2
    m3 = 100 * m3
    return m1, m2, m3
In [3]:
import os
from oceans.datasets import woa_subset

lon = -25.5
kw = dict(clim_type='annual', resolution='1deg', levels=slice(0, 40),
          llcrnrlat=-70, urcrnrlat=-10, llcrnrlon=lon, urcrnrlon=lon)

if os.path.isfile('./data/data.npz'):
    npz = np.load('./data/data.npz')
    depth = npz['p']
    SP, t, p = npz['SP'], npz['t'], npz['p']
    lon, lat = npz['lon'], npz['lat']
    salinity = woa_subset(var='salinity', **kw)
    salinity = salinity['OA Climatology']['annual'].squeeze()

    temperature = woa_subset(var='temperature', **kw)
    temperature = temperature['OA Climatology']['annual'].squeeze()

    SP = salinity.values.filled(fill_value=np.nan)
    t = temperature.values.filled(fill_value=np.nan)
    p = temperature.columns.values.astype(float)
    lat = temperature.index.values.astype(float)
    lon = np.array([lon]*len(lat))
    np.savez('./data/data.npz', t=t, SP=SP, lon=lon, lat=lat, p=p)

The T$-$S data we just downloaded above are in-situ temperature and practical salinity respectively.

Before we plot a $\Theta-S$ plot we need to compute conservative temperature and Absolute Salinity, if we choose to follow TEOS-10, or we should compute at least Potential Temperature ($\theta$) from EOS-80.

In [4]:
SP, p, lon, lat = np.broadcast_arrays(SP, p[None, :], lon[:, None], lat[:, None])
In [5]:
import gsw
SA = gsw.SA_from_SP(SP, p, lon, lat)
CT = gsw.CT_from_t(SA, t, p)
In [6]:
cores = np.array([[26.42, 5, 4.05], [37.14, 33.90, 35.07]])
sacw, aaiw, nadw = mixing(CT, SA, cores)
In [7]:
s = ma.masked_invalid(SA).mean(axis=0)
t = ma.masked_invalid(CT).mean(axis=0)
Te = np.linspace(t.min(), t.max(), 25)
Se = np.linspace(s.min(), s.max(), 25)

Tg, Sg = np.meshgrid(Te, Se)
sigma_theta = gsw.sigma0(Sg, Tg)
cnt = np.linspace(sigma_theta.min(), sigma_theta.max(), 10)
In [8]:
import brewer2mpl

Reds = brewer2mpl.get_map('Reds', 'Sequential', 9).mpl_colormap
Blues = brewer2mpl.get_map('Blues', 'Sequential', 9).mpl_colormap
Greens = brewer2mpl.get_map('Greens', 'Sequential', 9).mpl_colormap

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [9]:
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(14, 6), facecolor='w')
divider = make_axes_locatable(ax)

# Main figure.
percent = [50, 60, 70, 80, 90, 100]
m1 = ax.contourf(lat, p, sacw, percent, cmap=Reds, zorder=3)
m2 = ax.contourf(lat, p, aaiw, percent, cmap=Greens, zorder=2)
m3 = ax.contourf(lat, p, nadw, percent, cmap=Blues, zorder=1)

m1.set_clim(percent[0], percent[-1])
m2.set_clim(percent[0], percent[-1])
m3.set_clim(percent[0], percent[-1])

# Colorbars.
dy = 0.04
left, bottom, height, width = 0.14, 0.3, 0.13, 0.02

rect1 = [left, bottom, height, width]  # Top.
rect2 = [left, bottom - dy, height, width]  # Center.
rect3 = [left, bottom - 2*dy, height, width]  # Bottom.

cax1 = plt.axes(rect1, axisbg='w')
cax2 = plt.axes(rect2, axisbg='w')
cax3 = plt.axes(rect3, axisbg='w')

kw = dict(orientation='horizontal', extend='min')
cb1 = fig.colorbar(m1, cax=cax1, **kw)
cb2 = fig.colorbar(m2, cax=cax2, **kw)
cb3 = fig.colorbar(m3, cax=cax3, **kw)


new_labels = ['%s%%' % l.get_text() for l in cax3.xaxis.get_ticklabels()]

kw = dict(rotation=0, labelpad=20, y=1,
            verticalalignment='center', horizontalalignment='left')'SACW', **kw)"right")'AAIW', **kw)"right")'NADW', **kw)"right")

# Some tweaking.


new_labels = [r'%s $^\circ$S' % l.get_text()[1:] for l in

new_labels = [r'%s m' % l.get_text() for l in
ax.set_xlabel(u'WOA09 water masses at longitude %3.1f\u00B0' % lon[0, 0])

axin = inset_axes(ax, width="35%", height="35%", loc=4)
inmap = Basemap(projection='ortho', lon_0=lon[0, 0], lat_0=0,
                ax=axin, anchor='NE')
inmap.plot(lon, lat, 'r', latlon=True)

# T-S Diagram.
axTS = divider.append_axes("right", 2, pad=0.1)
cs = axTS.contour(Sg, Tg, sigma_theta, colors='grey', levels=cnt, zorder=1)
kw = dict(color='r', fontsize=14, fontweight='black')
axTS.text(35.00, 11.9, 'SACW', **kw)
axTS.text(34.50, 3.12, 'AAIW', **kw)
axTS.text(34.82, 1.72, 'NADW', **kw)

axTS.plot(s, t, 'k')
axTS.set_xlabel("Absolute Salinity [g kg$^{-1}$]")
axTS.set_ylabel("Conservative Temperature [$^\circ$C]", rotation=-90, labelpad=20)
axTS.set_title("$\Theta$$-$S Diagram")

Thanks to Clark for pointing this out. Don't forget to check the oce project and blog. (Maybe my next post should be with A-R-gh....)

In [10]:

This post was written as an IPython notebook. It is available for download or as a static html.

Creative Commons License
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at