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)}$$
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.
Input
-----
s : salinity
t : temperature [°C]
p : pressure [bars]
Return
------
rho density : [kg/m3]
Original Author
---------------
Olivier Le Calvé
E-mail : lecalve@isitv.univ-tln.fr
http://lecalve.univ-tln.fr/index.html
"""
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,
9.1697e-10]
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)}$$
Where:
- $\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}$
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
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_xlabel('Salinity')
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_xlabel('Salinity')
ax1.set_ylabel('Temperature [$^{\circ}$C]')
Hopefully that will motive students to inspect code and create some discussion on linear and non-linear version of the Equation of State.
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.
- UNESCO (1985) The international system of units (SI) in oceanography. UNESCO Technical Papers No.45, IAPSO Pub. Sci. No. 32, Paris, France.
- 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