Assignment: Ekman Pumping#

In this exercise we will apply the lesson learned about derivatives in Python to compute the Ekman transport in the ocean, and then compute the associated up- and downwelling.

# Import the tools we are going to need today:
import matplotlib.pyplot as plt  # plotting library
import numpy as np  # numerical library
import xarray as xr  # netCDF library
import cartopy  # Map projections libary
import cartopy.crs as ccrs  # Projections list
# Some defaults:
plt.rcParams['figure.figsize'] = (12, 5)  # Default plot size

Data#

We are going to work again with the ERA5_LowRes_MonthlyAvg_uvslp.nc and the ERA5_LowRes_Invariant.nc files that you already downloaded for a previous unit:

ds = xr.open_dataset('../data/ERA5_LowRes_MonthlyAvg_uvslp.nc')
dsi = xr.open_dataset('../data/ERA5_LowRes_Invariant.nc').isel(time=0)

Surface wind stress on the oceans#

If you remember the lesson (Oceans), the stress exerted by surface winds on the ocean can be reasonably well approximated by the bulk formula (Marshall and Plumb, p198):

\[(\tau_{wind_x}, \tau_{wind_y}) = \rho_{air} c_D ws_{10} (u_{10}, v_{10})\]

Where \(c_D\) is a bulk transfer coefficient for momentum (here we use \(c_D\) = 1.5\(\times\)10\(^{-3}\)), \(\rho_{air}\) is the denisty of air at the surface (here we use \(\rho_{air}\) = 1.225 kg m\(^3\)), \(ws\) is the wind speed at 10 m, and \(u_{10}\) , \(v_{10}\) are the components of the wind vector at 10 m height.

Q: \(\tau\) is expressed in N m\(^{-2}\). verify that \(c_D\) is a dimensionless quantity.

E: compute \(\overline{\tau_x}\) and \(\overline{\tau_y}\). This quantity only makes sense over the oceans. Read the invariant data and keep only the values of \(\overline{\tau_x}\) and \(\overline{\tau_y}\) over the oceans.

# Your answer here

E: compute the total wind stress \(\tau = \sqrt{\tau_x^2 + \tau_y^2}\). Plot it on a map using the colormap “Greens” (so that you can compare with the lecture). Add the (\(\tau_x\), \(\tau_y\)) vectors as arrows on top.

# Your answer here

Ekman transport#

The transport of water integrated over the Ekman layer is written (Marshall and Plumb, p200):

\[M_{ek} = \frac{\tau_{wind} \times \hat{z}}{f}\]

Or, in the component form:

\[m_{ek_{x}} = \frac{\tau_{wind_y}}{f} \]
\[m_{ek_{y}} = -\frac{\tau_{wind_x}}{f} \]

E: first, we need to compute the coriolis parameter \(f\):

f = 2. * 7.2921150e-5 * np.sin(np.deg2rad(ds.latitude))
f = f.where((np.abs(ds.latitude) > 3) & (np.abs(ds.latitude) < 87))  # Mask out the poles and equator regions

E: compute the components mek_x and mek_y. Plot the transport vectors on a map using the quiver() function. Discuss.

# Your answer here

Ekman pumping#

The convergence (or divergence) of the Ekman transport fields generates a vertical motion at the bottom of the Ekman Layer called Ekman pumping. It is written as (Marshall and Plumb, p204):

\[w_{ek} = \frac{1}{\rho_{0}} \left( \frac{\partial}{\partial x} \frac{\tau_{wind_{y}}}{f} - \frac{\partial}{\partial y} \frac{\tau_{wind_{x}}}{f} \right)\]

with \(\rho_{0}\) the density of water at the ocean surface, here assumed to be 1000 kg m\(^{-3}\).

As discussed in the lesson, we shouldn’t forget that we are using data defined in spherical coordinates. The partial derivatives in (x, y) are computed from spherical coordinates as:

\[\frac{\partial}{\partial x} = \frac{1}{R \cos \varphi}\frac{\partial}{\partial\lambda}\]
\[\frac{\partial}{\partial y} = \frac{1}{R}\frac{\partial}{\partial\varphi}\]

E: compute the partial derivatives \(\frac{\partial}{\partial \lambda} \frac{\tau_{wind_{y}}}{f}\) and \(-\frac{\partial}{\partial \varphi} \frac{\tau_{wind_{x}}}{f}\) out of tau_x and tau_y, and store them in the variables dlambda and dphi. Hint: you will need two calls to np.gradient() to do this computation, each call returning one variable that you will need and another that you won’t need. This second useless variable can be called “tmp” in each call for example.

_, dlambda = np.gradient(tau_y / f, -np.deg2rad(0.75), np.deg2rad(0.75))
dphi, _ = np.gradient(tau_x / f, -np.deg2rad(0.75), np.deg2rad(0.75))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_15451/1181974712.py in <cell line: 1>()
----> 1 _, dlambda = np.gradient(tau_y / f, -np.deg2rad(0.75), np.deg2rad(0.75))
      2 dphi, _ = np.gradient(tau_x / f, -np.deg2rad(0.75), np.deg2rad(0.75))

NameError: name 'tau_y' is not defined

Now we will convert these numpy arrays back to xarray DataArrays:

dlambda = mek_x*0 + dlambda
dphi = mek_x*0 + dphi
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_15451/2576805713.py in <cell line: 1>()
----> 1 dlambda = mek_x*0 + dlambda
      2 dphi = mek_x*0 + dphi

NameError: name 'mek_x' is not defined

E: compute the factor fx = \(\frac{1}{R \cos \varphi}\) and fy = \(\frac{1}{R}\), with an Earth radius of 6371000 m.

# Your answer here

E: you now have all elements to compute the Ekmann pumping “wek”. What is its unit? Convert it to m yr\(^{-1}\).

# Your answer here

E: plot wek on a map. Set vmin and vmax to -50 and +50 and use another colormap (for example “BrBG”). Compare your plot with the one computed by Marshall and Plumb (p205). Discuss the major features of the plot. Where is upwelling taking place? Downwelling? Can you explain the location and strength of these features? (click here to see the plot by Marshall and Plumb)

# Your answer here

To go further#

Read the few pages of chapter 10.3 in the book of Marshall and Plumb (available on OLAT). In particular, we are interested in their Fig. 10.21 and equation 10.20. Can you think about a strategy to realize the same plot with our data? Without implementing it (unless you really really want to!), describe the steps needed to come to this result.