Turning ripples into waves

Cushman-Roisin & Beckers chapter 3: International Equation of State

This post is part of the Introduction to Geophysical Fluid Dynamics, 2nd Edition post series on converting the MatlabTM examples to python. Chapter 3 contains just one mfile with an implementation for the Equation of State-80 (Chapter 3.3 -- Equation of State.)

But before we go any further it is worth mentioning that the Equation of State$^1$ (EOS-80) is superseded by the Thermodynamic Equation of State$^2$ (TEOS-10). Also, both EOS-80 and TEOS-10 have been implemented for FORTRAN, MatlabTM and of-course Python (gsw and seawater)

Now we can start the post! I usually ask my students to open and inspect all code they will use. Never using code in a "black box" manner.

For that goal, EOS-80 (eqs. 1), is a great exercise: easy to understand, quick to implement and with great accompanying documentation$^1$.

$$ \rho(S, t, p) = \dfrac{\rho(S, t, 0)}{\left[ \dfrac{1 - p}{K(S, t, p)} \right]} \quad\quad \text{(1)}$$

In [2]:
import numpy as np

def ies80(s, t, p=0):
    International Equation of State of sea water:
    rho = ies80(s, t, p)

    Function that calculates seawater density from temperature, salinity and pressure.
    Be aware of units.

    s : salinity
    t : temperature [°C]
    p : pressure [bars]
    rho density : [kg/m3]
    Original Author
    Olivier Le Calvé 
    E-mail :
    s, t, p = np.broadcast_arrays(s, t, p)

    r0_coef = [999.842594, 6.793952e-2, -9.09529e-3, 1.001685e-4, -1.120083e-6,
               6.536332e-9, 8.24493e-1, -4.0899e-3, 7.6438e-5, -8.2467e-7,
               5.3875e-9, -5.72466e-3, 1.0227e-4, -1.6546e-6, 4.8314e-4]
    rho_0 = (np.polyval(r0_coef[5::-1], t) +
             np.polyval(r0_coef[10:5:-1], t) * s +
             np.polyval(r0_coef[13:10:-1], t) * s**1.5 +
             r0_coef[14] * s**2)
    K_coef = [19652.21, 148.4206, -2.327105, 1.360447e-2, -5.155288e-5,
              3.239908, 1.43713e-3, 1.16092e-4, -5.77905e-7, 8.50935e-5,
              -6.12293e-6, 5.2787e-8, 54.6746, -0.603459, 1.09987e-2,
              -6.1670e-5, 7.944e-2, 1.6483e-2, -5.3009e-4, 2.2838e-3,
              -1.0981e-5, -1.6078e-6, 1.91075e-4, -9.9348e-7, 2.0816e-8,
    K = (np.polyval(K_coef[4::-1], t) +
         np.polyval(K_coef[8:4:-1], t) * p +
         np.polyval(K_coef[11:8:-1], t) * p**2 +
         np.polyval(K_coef[15:11:-1], t) * s +
         np.polyval(K_coef[18:15:-1], t) * s**1.5 +
         np.polyval(K_coef[21:18:-1], t) * p * s +
         K_coef[22] * p * s**1.5 +
         np.polyval(K_coef[25:22:-1], t) * p**2 * s)
    rho = rho_0 / (1-p / K)

    return rho

The book asks the students to compare the calculated values from the non-linear version with a simpler linear (eq. 2) version of the Equation of State. Also, tell students to estimate $\alpha$ and $\beta$ (Thermal expansion and Salinity contraction coefficients).

I will leave that for my class, but here I'll plot a T-S diagram with contours from both linear and non-linear versions of EOS.

$$\rho = \rho_o \left[ 1-\alpha \left( T - T_o \right) + \beta\left( S - S_o \right) \right] \quad\quad \text{(2)}$$


  • $\rho_o = 1028$ kg m$^{-3}$
  • $T_o = 10^{\circ}$C = 283 K
  • $S_o = 35$
  • $\alpha = 1.7 \times 10^{-4}$ K$^{-1}$
  • $\beta = 7.6 \times 10^{-4}$
In [3]:
def linear_eos(S, T):
    S, T = np.broadcast_arrays(S, T)
    rho_o, To, So = 1028, 10, 35
    alpha, beta = 1.7e-4, 7.6e-4
    rho = rho_o * (1 - alpha * (T-To) + beta * (S-So))
    return rho
In [4]:
import matplotlib.pyplot as plt

S = np.linspace(0, 42, 100)
t = np.linspace(-2, 40, 100)
S, t = np.meshgrid(S, t)

sigma_full = ies80(S, t, p=0) - 1000
sigma_linear = linear_eos(S, t) - 1000

cnt = np.arange(-7, 35, 5)

fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(8, 4))
cs0 = ax0.contour(S, t, sigma_full, colors='blue', levels=cnt)
cs1 = ax0.contour(S, t, sigma_linear, colors='red', alpha=0.5, levels=cnt)
ax0.clabel(cs0, fontsize=9, inline=1, fmt='%2i')
ax0.clabel(cs1, fontsize=9, inline=1, fmt='%2i')
ax0.set_ylabel('Temperature [$^{\circ}$C]')

ax1.axis([32, 34, 20, 30.0])
cs0 = ax1.contour(S, t, sigma_full, colors='blue', levels=[21, 23, 24])
cs1 = ax1.contour(S, t, sigma_linear, colors='red', alpha=0.5, levels=[21, 23, 24])
ax1.clabel(cs0, fontsize=9, inline=1, fmt='%2i')
ax1.clabel(cs1, fontsize=9, inline=1, fmt='%2i')
ax1.set_ylabel('Temperature [$^{\circ}$C]')
<matplotlib.text.Text at 0x7f91c3f1ae50>

Hopefully that will motive students to inspect code and create some discussion on linear and non-linear version of the Equation of State.

In [5]:

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
  1. UNESCO (1985) The international system of units (SI) in oceanography. UNESCO Technical Papers No.45, IAPSO Pub. Sci. No. 32, Paris, France.
  2. McDougall, T.J., R. Feistel, F. J. Millero, D. R. Jackett, D. G. Wright, B. A. King, G. M. Marion, C. T. A. Chen and P. Spitzer, 2009: Calculation of the Thermophysical Properties of Seawater, Global Shipbased Repeat Hydrography Manual, IOCCP Report No. 14, ICPO Publication Series no. 134