This post is part of the
Introduction to Geophysical Fluid Dynamics, 2nd Edition
post series on converting the MatlabTM examples to python. Chapter 8 is about
Ekman Dynamics. This chapter contains two m_files ekmanh.m and
spectralekman.m where the former is the exact solution and the latter the
spectral solution for the Ekman layer with constant eddy viscosity.
Starting from the following system of equation:
$$ -fv = -\dfrac{1}{\rho_o} \dfrac{\partial p}{\partial x} +\nu_E \dfrac{\partial^2u}{\partial z^2} $$ $$ +fu = -\dfrac{1}{\rho_o} \dfrac{\partial p}{\partial y} +\nu_E \dfrac{\partial^2v}{\partial z^2} $$ $$ 0 = -\dfrac{1}{\rho_o} \dfrac{\partial p}{\partial z} $$
If we align the x–axis with the direction of $\bar{u}$ we can write the boundary conditions as:
$$ \text{at } z=0 \text{ (bottom): } u = v = 0 $$ $$ \text{at } z>> d \text{ (interior): } u = \bar{u},\quad v = 0,\quad p=\bar{p}(x, y) $$
$$ 0 = -\dfrac{1}{\rho_o} \dfrac{\partial \bar{p}}{\partial x} $$
$$ f\bar{u} = -\dfrac{1}{\rho_o} \dfrac{\partial \bar{p}}{\partial y} = \text{constant} $$
Plugin these in the equations above:
$$ -fv = \nu E \dfrac{d^2u}{d z^2} $$ $$ f(u - \bar{u}) = \nu E \dfrac{d^2v}{d z^2} $$
we will seek a solution in the form of $u = \bar{u} + A\exp(\lambda z)$ and $v = \exp(\lambda z)$, we know that $\lambda$ is $\nu^2\lambda^4+f^2 = 0$,
$$\lambda = \pm(1\pm i)\dfrac{1}{d},$$
where $d$ is,
$$d = \sqrt{\dfrac{2\nu_E}{f}}$$
Below are the python version and the parameters used for this problem. For further information read chapter 8 of the book.
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
def spectralekman(N):
"""Function that calculates the spectral solution of the Ekman layer with constant eddy viscosity.
The parameter N is the number of modes retained. If run with N=5 you can recreate figure 8--12.
For higher values if N you can see the behaviour of Figure 8--13. Two see how a solution behaves
for the wind-stress case, use isol=True."""
# 1st.
j = np.arange(N)
s = 2 * np.sqrt(2 * h) / (np.pi * (2 * (j+1)-1))
phih = np.sqrt(2 / h) * (-1)**(j)
eigen = (2 * (j+1)-1)**2 * np.pi * np.pi * nu / (4 * h**2)
FX = -pressx * s + taux * phih
FY = -pressy * s + tauy * phih
a = (eigen * FX + f * FY) / (eigen**2 + f**2)
b = (eigen * FY - f * FX) / (eigen**2 + f**2)
# 2nd.
U = np.zeros(KM)
V = np.zeros_like(U)
uex = np.zeros_like(U)
vex = np.zeros_like(U)
z = (np.arange(1, KM + 1) - 1) * h / (KM - 1)
for j in np.arange(N)[::-1]:
U += a[j] * np.sqrt(2 / h) * np.sin((2 * (j+1)-1) * np.pi / 2 * z / h)
V += b[j] * np.sqrt(2 / h) * np.sin((2 * (j+1)-1) * np.pi / 2 * z / h)
uex, vex = ekmanh(h / dek, z / dek)
if isol:
uex = uex * taux / nu * dek
vex = vex * taux / nu * dek
else:
uex = -(1 + uex) * pressy / f
vex = -vex * pressy / f
err = 0
for k in range(0, KM):
err += (U[k] - uex[k])**2 + (V[k] - vex[k])**2
err = np.sqrt(err / KM)
aa = np.log(a**2)
xx = np.log(np.arange(1, N+1))
return uex, vex, U, V, z, err, aa, xx
def ekmanh(h, z):
"""Function called by `spectralekman`. It provides the exact solution of the problem."""
if isol: # No pressure gradient, du/dz=1 dv/dz=0.
u = ((np.sin(h) * np.sinh(h) * (np.cosh(z) * np.sin(z) - np.cos(z) * np.sinh(z)) +
np.cos(h) * np.cosh(h) * (np.cosh(z) * np.sin(z) + np.cos(z) * np.sinh(z))) /
(np.cos(2 * h) + np.cosh(2 * h)))
v = ((-np.sin(h) * np.sinh(h) * (np.cosh(z) * np.sin(z) + np.cos(z) * np.sinh(z)) +
np.cos(h) * np.cosh(h) * (np.cosh(z) * np.sin(z) -np.cos(z) * np.sinh(z))) /
(np.cos(2 * h) + np.cosh(2 * h)))
else: # No wind stress, ug = 1, vg=0.
u = (-np.exp(-z) * (np.exp(2 * h) * ( 1 + np.exp(2 * z)) * np.cos(2 * h - z) +
(np.exp(4 * h) - np.exp(2 * z)) * np.cos(z)) /
(1 + np.exp(4 * h) + 2 * np.exp(2 * h) * np.cos(2 * h)))
v = (np.exp(-z) * (np.exp(2 * h) * (-1 + np.exp(2 * z)) * np.sin(2 * h - z) +
(np.exp(4 * h) - np.exp(2 * z)) * np.sin(z)) /
(1 + np.exp(4 * h) + 2 * np.exp(2 * h) * np.cos(2 * h)))
return u, v
isol = False
title = "Wind stress solution." if isol else "Pressure gradient solution."
nu = 0.01
h = 100
f = 0.0001
tauy = 0
pressx = 0
if isol:
taux = 10e-6 * 10 * 10 * 1
pressy = 0
else:
taux = 0
pressy = 9.81 * 1/1000000
KM = 400
dek = np.sqrt(2 * nu / f)
uex, vex, U, V, z, err, aa, xx = spectralekman(N=5)
fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5))
ax0.plot(uex, z, vex, z, U, z, 'r-.', V, z, 'k-.')
ax0.set_xlabel('velocity')
ax0.set_ylabel('z')
ax1.bar(xx, aa + 40, width=0.1)
ax1.set_xlabel('$\log(j)$')
ax1.set_ylabel(r'$\log(a_j^2)$')
_ = fig.suptitle('Ekman Layer solution for %s' % title)
isol = True
title = "Wind stress solution." if isol else "Pressure gradient solution."
nu = 0.01
h = 100
f = 0.0001
tauy = 0
pressx = 0
if isol:
taux = 10e-6 * 10 * 10 * 1
pressy = 0
else:
taux = 0
pressy = 9.81 * 1/1000000
KM = 400
dek = np.sqrt(2 * nu / f)
uex, vex, U, V, z, err, aa, xx = spectralekman(N=5)
fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5))
ax0.plot(uex, z, vex, z, U, z, 'r-.', V, z, 'k-.')
ax0.set_xlabel('velocity')
ax0.set_ylabel('z')
ax1.bar(xx, aa + 40, width=0.1)
ax1.set_xlabel('$\log(j)$')
ax1.set_ylabel(r'$\log(a_j^2)$')
_ = fig.suptitle('Ekman Layer solution for %s' % title)
I am using N=5 for both Wind Stress and Pressure Gradient solutions. Play
with the parameters according to the exercises in the book. Have fun
exploring all the possible solutions!
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/.