Open In Colab   Open in Kaggle

Tutorial 4: Oceanic Wind-driven Circulation#

Week 1, Day 2: Ocean and Atmospheric Reanalysis

Content creators: Abigail Bodner

Content reviewers: Katrina Dobson, Danika Gupta, Maria Gonzalez, Will Gregory, Nahid Hasan, Sherry Mi, Beatriz Cosenza Muralles, Jenna Pearson, Chi Zhang, Ohad Zivan

Content editors: Brodie Pearson, Jenna Pearson, Chi Zhang, Ohad Zivan

Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan

Our 2023 Sponsors: NASA TOPS and Google DeepMind

Tutorial Objectives#

The ocean’s motion is driven by radiation from the sun, winds, and various sources & sinks of fresh water (precipitation, rivers, melting and freezing ice). In the previous tutorial, you quantified the surface winds around the world. The surface winds drags on the surface of the ocean which results in ocean transport, known as Ekman transport.

In this tutorial you will use the ECCO (Estimating the Circulation and Climate of the Ocean) reanalysis data, to visualize the ocean’s surface currents and to compare these currents against local atmospheric conditions.

At the end of this tutorial you will be able to

  • Access and manipulate ocean reanalysis data

  • Plot annual mean ocean surface currents and atmospheric surface winds

  • Compare oceanic and atmospheric circulation patterns

Setup#

# installations ( uncomment and run this cell ONLY when using google colab or kaggle )

# !pip install pythia_datasets
# !pip install cartopy
# !pip install geoviews
# from intake import open_catalog
import matplotlib.pyplot as plt
import matplotlib
import os
import pooch
import tempfile
import numpy as np
import xarray as xr
import warnings
from cartopy import crs as ccrs, feature as cfeature

#  suppress warnings issued by Cartopy when downloading data files
warnings.filterwarnings("ignore")

Helper functions#

Hide code cell source
# @title Helper functions

def pooch_load(filelocation=None, filename=None, processor=None):
    shared_location = "/home/jovyan/shared/Data/tutorials/W1D2_StateoftheClimateOceanandAtmosphereReanalysis"  # this is different for each day
    user_temp_cache = tempfile.gettempdir()

    if os.path.exists(os.path.join(shared_location, filename)):
        file = os.path.join(shared_location, filename)
    else:
        file = pooch.retrieve(
            filelocation,
            known_hash=None,
            fname=os.path.join(user_temp_cache, filename),
            processor=processor,
        )

    return file

Video 1: Ocean Wind-driven Circulation#

Section 1: Load Ocean & Atmosphere Reanalysis Data#

Here you will load atmospheric near-surface winds (at 10-meter height), and then load the oceanic surface currents from ECCO reanalysis data

Note, each of these variables is a velocity with two components (zonal and meridional). These two velocity components must be loaded separately for each variable, so you will load four datasets.

# load data: atmospheric 10m wind from ECCO
# wind in east/west direction labeled here as 'u'
fname_atm_wind_u = "wind_evel_monthly_2016.nc"
url_atm_wind_u = "https://osf.io/ke9yp/download"
atm_wind_u = xr.open_dataarray(pooch_load(url_atm_wind_u, fname_atm_wind_u))
atm_wind_u
Downloading data from 'https://osf.io/ke9yp/download' to file '/tmp/wind_evel_monthly_2016.nc'.
SHA256 hash of downloaded file: 728abbbcf931d178e95f1fdf6df61da45d8ec21627a9b3cc30d65be55b0e8d6d
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.DataArray 'EXFewind' (time: 12, latitude: 360, longitude: 720)> Size: 25MB
[3110400 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 2016-01-16T12:00:00 ... 2016-12-16T1...
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
Attributes:
    units:      m/s
    long_name:  eastward 10-m wind velocity over open water, >0 increases eVel
    mate:       EXFnwind
# wind in north/south direction labeled here as 'v'
fname_atm_wind_v = "wind_nvel_monthly_2016.nc"
url_atm_wind_v = "https://osf.io/9zkgd/download"
atm_wind_v = xr.open_dataarray(pooch_load(url_atm_wind_v, fname_atm_wind_v))
atm_wind_v
Downloading data from 'https://osf.io/9zkgd/download' to file '/tmp/wind_nvel_monthly_2016.nc'.
SHA256 hash of downloaded file: 0d5b2a21c84483395a48bfa2a55e022d4c4cb08a07e29e64a750b0b08f64413a
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.DataArray 'EXFnwind' (time: 12, latitude: 360, longitude: 720)> Size: 25MB
[3110400 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 2016-01-16T12:00:00 ... 2016-12-16T1...
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
Attributes:
    units:      m/s
    long_name:  northward 10-m wind velocity over open water, >0 increases nVel
    mate:       EXFewind
# load data: oceanic surface current from ECCO
# current in east/west direction labeled here as 'u'
fname_ocn_surface_current_u = "evel_monthly_2016.nc"
url_ocn_surface_current_u = "https://osf.io/ndx5a/download"
ocn_surface_current_u = xr.open_dataarray(
    pooch_load(url_ocn_surface_current_u, fname_ocn_surface_current_u)
)
ocn_surface_current_u
Downloading data from 'https://osf.io/ndx5a/download' to file '/tmp/evel_monthly_2016.nc'.
SHA256 hash of downloaded file: 452bdf8b74e49207c53db24cdb82fd7e712fcc48495339fbe6157dd3547f085a
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.DataArray 'EVEL' (time: 12, latitude: 360, longitude: 720)> Size: 25MB
[3110400 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 2014-01-16T12:00:00 ... 2014-12-16T1...
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
Attributes:
    units:          m/s
    long_name:      Eastward Component of Velocity (m/s)
    standard_name:  eastward_sea_water_velocity
# current in east/west direction labeled here as 'v'
fname_ocn_surface_current_v = "nvel_monthly_2016.nc"
url_ocn_surface_current_v = "https://osf.io/qa9ex/download"
ocn_surface_current_v = xr.open_dataarray(
    pooch_load(url_ocn_surface_current_v, fname_ocn_surface_current_v)
)
ocn_surface_current_v
Downloading data from 'https://osf.io/qa9ex/download' to file '/tmp/nvel_monthly_2016.nc'.
SHA256 hash of downloaded file: de285efede0d4a4b919cdc174a68941449a875f5dcf7f70bb9c2072ee9b5d690
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.DataArray 'NVEL' (time: 12, latitude: 360, longitude: 720)> Size: 25MB
[3110400 values with dtype=float64]
Coordinates:
  * time       (time) datetime64[ns] 96B 2014-01-16T12:00:00 ... 2014-12-16T1...
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
Attributes:
    units:          m/s
    long_name:      Northward Component of Velocity (m/s)
    standard_name:  northward_sea_water_velocity

Section 1.1: Exploring the Reanalysis Data#

Let’s examine the time (or temporal/output) frequency, which descibes the rate at which the reanalysis data is provided, for one of the ECCO variables (atm_wind_u).

Note that all the variables should have the same output frequency.

atm_wind_u.time
<xarray.DataArray 'time' (time: 12)> Size: 96B
array(['2016-01-16T12:00:00.000000000', '2016-02-15T12:00:00.000000000',
       '2016-03-16T12:00:00.000000000', '2016-04-16T00:00:00.000000000',
       '2016-05-16T12:00:00.000000000', '2016-06-16T00:00:00.000000000',
       '2016-07-16T12:00:00.000000000', '2016-08-16T12:00:00.000000000',
       '2016-09-16T00:00:00.000000000', '2016-10-16T12:00:00.000000000',
       '2016-11-16T00:00:00.000000000', '2016-12-16T12:00:00.000000000'],
      dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 96B 2016-01-16T12:00:00 ... 2016-12-16T12:...
Attributes:
    long_name:  center time of averaging period
    bounds:     time_bnds
    axis:       T

Questions 1.1#

  1. Why do you think the atmospheric reanalysis dataset from the previous tutorial (ERA5) has higher output frequency (hourly) than the ECCO ocean reanalysis dataset (daily)?

  2. What can you infer about the role of these two systems at different timescales in the climate system?

Click for solution

Section 2: Plotting the Annual Mean of Global Surface Wind Stress#

In this section you will create global maps displaying the annual mean of atmospheric 10m winds.

First, you should compute the annual mean of the surface wind variables. You can do so by averaging over the time dimension using .mean(dim='time'). Since you have monthly data spanning only one year, .mean(dim='time') will give the annual mean for the year 2016.

# compute the annual mean of atm_wind_u
atm_wind_u_an_mean = atm_wind_u.mean(dim="time")
atm_wind_u_an_mean
<xarray.DataArray 'EXFewind' (latitude: 360, longitude: 720)> Size: 2MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [1.14805873, 1.14805873, 1.14805873, ..., 1.14805873, 1.14805873,
        1.14805873],
       [1.03043525, 1.03043525, 1.03043525, ..., 1.03043525, 1.03043525,
        1.03043525],
       [0.60513089, 0.60513089, 0.60513089, ..., 0.60513089, 0.60513089,
        0.60513089]])
Coordinates:
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
# take the annual mean of atm_wind_stress_v
atm_wind_v_an_mean = atm_wind_v.mean(dim="time")
atm_wind_v_an_mean
<xarray.DataArray 'EXFnwind' (latitude: 360, longitude: 720)> Size: 2MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.29360731, 0.29360731, 0.29360731, ..., 0.29360731, 0.29360731,
        0.29360731],
       [0.25386951, 0.25386951, 0.25386951, ..., 0.25386951, 0.25386951,
        0.25386951],
       [0.64899328, 0.64899328, 0.64899328, ..., 0.64899328, 0.64899328,
        0.64899328]])
Coordinates:
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8

You are now almost ready to plot!

However, you currently have seperate zonal and meridional wind velocity components \((u,v)\). An effective way of visualizing the total surface wind stress is to create a global map of the magnitude and direction of the wind velocity vector. This type of plot is known as a vector field. A vector is a special mathematical quantity that has both magnitude and direction, just like the wind! The velocity components describe the intensity of wind blowing in the zonal (\(u\)) or meridional (\(v\)) directions. Specifically, wind can blow eastward (positive \(u\)) or westward (negative \(u\)), as well as northward (positive \(v\)) or southward (negative \(v\)).

The total velocity vector is the vector sum of these two components and exhibits varying magnitude and direction. The magnitude (\(||u||\)) and direction (\(θ\)) of the total velocity vector can be determined using the following equations:

(3)#\[\begin{align} ||u|| = \sqrt{u^2 + v^2}, \ \ \ \ \theta = tan^{-1}(\frac{v}{u}) \end{align}\]

When plotting a vector field using a computer, it is commonly referred to as a quiver plot. In our case, we will utilize a quiver function created by Ryan Abernathey that calculates the magnitude and direction of the total velocity vector based on the given zonal and meridional components.

We will overlay the quiver plot on top of the annual mean ocean surface temperature (labeled here as theta).

fname_surface_temp = "surface_theta.nc"
url_fname_surface_temp = "https://osf.io/98ksr/download"
ocn_surface_temp = xr.open_dataarray(
    pooch.retrieve(url_fname_surface_temp, known_hash=None, fname=fname_surface_temp)
)
ocn_surface_temp
Downloading data from 'https://osf.io/98ksr/download' to file '/home/runner/.cache/pooch/surface_theta.nc'.
SHA256 hash of downloaded file: 2d8db6d497acec50c4d0753ba69c610231dc02294ab1dc99640e0fce0285d369
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.DataArray 'THETA' (latitude: 360, longitude: 720)> Size: 2MB
[259200 values with dtype=float64]
Coordinates:
    i          (longitude) int64 6kB ...
    k          int64 8B ...
    j          (latitude) int64 3kB ...
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
    Z          float32 4B ...
# longitude ad latitude coordinates for plotting
lon = atm_wind_u_an_mean.longitude
lat = atm_wind_u_an_mean.latitude

# calculate magnitude of total velocity
mag = (atm_wind_u_an_mean**2 + atm_wind_v_an_mean**2) ** 0.5

# coarsen the grid so the arrows are distinguishable by only selecting
# some longitudes and latitudes defined by sampling_x and sampling_y.
slx = slice(None, None, 20)
sly = slice(None, None, 20)
sl2d = (sly, slx)

fig, ax = plt.subplots(
    figsize=(12, 6), subplot_kw={"projection": ccrs.Robinson(central_longitude=180)}
)

c = ax.contourf(lon, lat, ocn_surface_temp, alpha=0.5)

# plot quiver arrows indicating vector direction (winds are in blue, alpha is for opacity)
q = ax.quiver(
    lon[slx],
    lat[sly],
    atm_wind_u_an_mean[sl2d],
    atm_wind_v_an_mean[sl2d],
    color="b",
    alpha=0.5,
)

ax.quiverkey(q, 175, 95, 5, "5 m/s", coordinates="data")

fig.colorbar(c, label="Sea Surface Temperature (degC)")
<matplotlib.colorbar.Colorbar at 0x7f88a8afcfd0>
../../../_images/c0cf70e7edc8056cbe48abd4a701a267e60e95749b389f76e70d845e1f0fcf39.png

Section 3: Comparing Global Maps of Surface Currents and Winds#

In this section you will compute the annual mean of the ocean surface currents, similar to your above analyses of atmospheric winds, and you will create a global map that shows both of these variables simultaneously.

# take the annual mean of ocn_surface_current_u
ocn_surface_current_u_an_mean = ocn_surface_current_u.mean(dim="time")
ocn_surface_current_u_an_mean
<xarray.DataArray 'EVEL' (latitude: 360, longitude: 720)> Size: 2MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.02164101, 0.02164101, 0.02164101, ..., 0.02164101, 0.02164101,
        0.02164101],
       [0.01952917, 0.01952917, 0.01952917, ..., 0.01952917, 0.01952917,
        0.01952917],
       [0.01100874, 0.01100874, 0.01100874, ..., 0.01100874, 0.01100874,
        0.01100874]])
Coordinates:
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
# take the annual mean of ocn_surface_current_v
ocn_surface_current_v_an_mean = ocn_surface_current_v.mean(dim="time")
ocn_surface_current_v_an_mean
<xarray.DataArray 'NVEL' (latitude: 360, longitude: 720)> Size: 2MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.00707205, 0.00707205, 0.00707205, ..., 0.00707205, 0.00707205,
        0.00707205],
       [0.00721859, 0.00721859, 0.00721859, ..., 0.00721859, 0.00721859,
        0.00721859],
       [0.01465764, 0.01465764, 0.01465764, ..., 0.01465764, 0.01465764,
        0.01465764]])
Coordinates:
  * latitude   (latitude) float64 3kB -89.75 -89.25 -88.75 ... 88.75 89.25 89.75
  * longitude  (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8

Let’s add ocean surface currents to the previous plot above, using red quivers. Note the scale of the arrows for the ocean and atmosphere are very different.

# longitude ad latitude coordinates for plotting
lon = atm_wind_u_an_mean.longitude
lat = atm_wind_u_an_mean.latitude

# calculate magnitude of total velocity
mag = (atm_wind_u_an_mean**2 + atm_wind_v_an_mean**2) ** 0.5

# coarsen the grid so the arrows are distinguishable by only selecting
# some longitudes and latitudes defined by sampling_x and sampling_y.
slx = slice(None, None, 20)
sly = slice(None, None, 20)
sl2d = (sly, slx)

# fig, ax = plt.subplots(**kwargs)
fig, ax = plt.subplots(
    figsize=(12, 6), subplot_kw={"projection": ccrs.Robinson(central_longitude=180)}
)

c = ax.contourf(lon, lat, ocn_surface_temp, alpha=0.5)

# plot quiver arrows indicating vector direction (winds are in blue, alpha is for opacity)
q1 = ax.quiver(
    lon[slx],
    lat[sly],
    atm_wind_u_an_mean[sl2d],
    atm_wind_v_an_mean[sl2d],
    color="b",
    alpha=0.5,
)

# plot quiver arrows indicating vector direction (ocean currents are in red, alpha is for opacity)
q2 = ax.quiver(
    lon[slx],
    lat[sly],
    ocn_surface_current_u_an_mean[sl2d],
    ocn_surface_current_v_an_mean[sl2d],
    color="r",
    alpha=0.5,
)

ax.quiverkey(q1, 175, 95, 5, "5 m/s", coordinates="data")

ax.quiverkey(q2, 150, 95, 0.05, "5 cm/s", coordinates="data")

fig.colorbar(c, label="Sea Surface Temperature (degC)")
<matplotlib.colorbar.Colorbar at 0x7f88a792de80>
../../../_images/c9d15d2d0eafacce55d9197f744d048577ac29de3db91922efc4aeeac31e5916.png

Questions 3#

You may notice that the surface currents (red) are typically not aligned with the wind direction (blue). In fact, the surface ocean currents flow at an angle of approximately 45 degrees to the wind direction!

The combination of Coriolis force and the friction between ocean layers causes the wind-driven currents to turn and weaken with depth, eventually dissapearing entirely at depth. The resulting current profile is called the Ekman Spiral, and the depth over which the spiral is present is called the Ekman layer. While the shape of this spiral can vary in time and space, the depth-integrated transport of water within the Ekman layer is called Ekman transport. Ekman transport is always 90 degrees to the right of the wind in the Northern Hemisphere, and 90 degrees to the left of the wind in the Southern Hemisphere. Under certain wind patterns or near coastlines, Ekman transport can lead to Ekman Pumping, where water is upwelled or downwelled depending on the wind direction. This process is particularily important for coastal ecosystems that rely on this upwelling (which is often seasonal) to bring nutrient-rich waters near the surface. These nutrients fuel the growth of tiny marine plants that support the entire food chain, upon which coastal economies often rely. These tiny plants are also responsible for most of the oxygen we breathe!

  1. Do you observe a deflection of the surface currents relative to the wind?

  2. Is this deflection the same in both the Northern and Southern Hemispheres?

Click for solution

Summary#

In this tutorial, you explored how atmospheric winds can shape ocean surface currents.

You calculated the surface wind vector an assessed its influence on surface ocean currents called the wind-driven circulation.

This circulation influences global climate (e.g., moving heat between latitudes) and biogeochemical cycles (e.g., by driving upwelling of nutrient rich water). In the next tutorial you will explore another mechanism for ocean movement, density-driven circulation.

Resources#

Data for this tutorial can be accessed here.