# Spectral Solution for the Ekman Layer

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.

In [2]:
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

In [3]:
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)

In [4]:
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!

In [5]:
HTML(html)

Out[5]:

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