Skip to article frontmatterSkip to article content

Phytoplankton nutrient limitation

A centric diatom, another type of phytoplankton. Art credit: Kristen Krumhardt

Overview

In previous notebooks, we explored the distribution of different nutrients in the ocean. Here we examine how the growth of phytoplankton communities is limited by these nutrient distributions.

  1. General setup
  2. Subsetting
  3. Processing - long-term mean
  4. Mapping nutrient limitation at the surface
  5. Mapping biomass-weighted nutrient limitation in the top 100m
  6. Making monthly climatology maps

Prerequisites

ConceptsImportanceNotes
MatplotlibNecessary
Intro to CartopyNecessary
Dask CookbookHelpful
Intro to XarrayHelpful
  • Time to learn: 30 min

Imports

import xarray as xr
import glob
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors
import cartopy
import cartopy.crs as ccrs
import pop_tools
from dask.distributed import LocalCluster
import s3fs
import netCDF4

from module import adjust_pop_grid

General setup (see intro notebooks for explanations)

Connect to cluster

cluster = LocalCluster()
client = cluster.get_client()

Bring in POP grid utilities

ds_grid = pop_tools.get_grid('POP_gx1v7')
lons = ds_grid.TLONG
lats = ds_grid.TLAT
depths = ds_grid.z_t * 0.01

Load the data

jetstream_url = 'https://js2.jetstream-cloud.org:8001/'

s3 = s3fs.S3FileSystem(anon=True, client_kwargs=dict(endpoint_url=jetstream_url))

# Generate a list of all files in CESM folder
s3path = 's3://pythia/ocean-bgc/cesm/g.e22.GOMIPECOIAF_JRA-1p4-2018.TL319_g17.4p2z.002branch/ocn/proc/tseries/month_1/*'
remote_files = s3.glob(s3path)

# Open all files from folder
fileset = [s3.open(file) for file in remote_files]

# Open with xarray
ds = xr.open_mfdataset(fileset, data_vars="minimal", coords='minimal', compat="override", parallel=True,
                       drop_variables=["transport_components", "transport_regions", 'moc_components'], decode_times=True)

ds

Subsetting

variables =['sp_Fe_lim_Cweight_avg_100m','sp_P_lim_Cweight_avg_100m','sp_N_lim_Cweight_avg_100m',
             'diat_Fe_lim_Cweight_avg_100m', 'diat_P_lim_Cweight_avg_100m','diat_N_lim_Cweight_avg_100m',
             'diat_SiO3_lim_Cweight_avg_100m','diaz_P_lim_Cweight_avg_100m',
             'diaz_Fe_lim_Cweight_avg_100m','cocco_Fe_lim_Cweight_avg_100m','cocco_C_lim_Cweight_avg_100m','cocco_N_lim_Cweight_avg_100m',
             'cocco_P_lim_Cweight_avg_100m','sp_Fe_lim_surf','sp_P_lim_surf','sp_N_lim_surf',
             'diat_Fe_lim_surf', 'diat_P_lim_surf','diat_N_lim_surf','diat_SiO3_lim_surf',
             'diaz_P_lim_surf','cocco_Fe_lim_surf','cocco_C_lim_surf','cocco_N_lim_surf',
             'cocco_P_lim_surf','diaz_Fe_lim_surf']
keep_vars=['z_t','z_t_150m','dz','time_bound', 'time', 'TAREA','TLAT','TLONG'] + variables
ds = ds.drop_vars([v for v in ds.variables if v not in keep_vars])

Processing - long-term mean

Pull in the function we defined in the nutrients notebook...

def year_mean(ds):
    """
    Properly convert monthly data to annual means, taking into account month lengths.
    Source: https://ncar.github.io/esds/posts/2021/yearly-averages-xarray/
    """
    
    # Make a DataArray with the number of days in each month, size = len(time)
    month_length = ds.time.dt.days_in_month

    # Calculate the weights by grouping by 'time.season'
    weights = (
        month_length.groupby("time.year") / month_length.groupby("time.year").sum()
    )

    # Test that the sum of the weights for each season is 1.0
    np.testing.assert_allclose(weights.groupby("time.year").sum().values, np.ones((len(ds.groupby("time.year")), )))

    # Calculate the weighted average
    return (ds * weights).groupby("time.year").sum(dim="time")

Take the long-term mean of our data set. We process monthly to annual with our custom function, then use xarray’s built-in .mean() function to process from annual data to a single mean over time, since each year is the same length.

ds_year = year_mean(ds).mean("year")

Mapping nutrient limitation at the surface

Phytoplankton need a specific ratio of nutrients to grow and produce organic matter. In general, this is known as the Redfield ratio, first proposed by Redfield et al., 1963, and is approximately 106 C:16 N:1 P. Micronutrients like silicate and iron are also needed in more variable amounts depending on plankton type. To learn more about nutrient limitation, see Sarmiento and Gruber Chapter 4: Organic Matter Production. Our dataset uses a numerical notation to specify which nutrient is limiting in each area for each phytoplankton functional type, as specified below:

  • 0 = PO4
  • 1 = Fe
  • 2 = NO3 (only for small phytoplankton and diatoms)
  • 3 = Si (only for diatoms)
  • 3 = C (only for coccolithophores)

To turn this information into a single array, we concatenate the limitation terms along the nutrient dimension for each phytoplankton functional type.

limarray_sp=xr.concat((ds_year.sp_P_lim_surf, ds_year.sp_Fe_lim_surf,ds_year.sp_N_lim_surf),dim='nutrient')
limarray_diat=xr.concat((ds_year.diat_P_lim_surf, ds_year.diat_Fe_lim_surf, ds_year.diat_N_lim_surf, ds_year.diat_SiO3_lim_surf),dim='nutrient')
limarray_diaz=xr.concat((ds_year.diaz_P_lim_surf, ds_year.diaz_Fe_lim_surf),dim='nutrient')
limarray_cocco=xr.concat((ds_year.cocco_P_lim_surf, ds_year.cocco_Fe_lim_surf, ds_year.cocco_N_lim_surf, ds_year.cocco_C_lim_surf),dim='nutrient')
                          
most_lim_sp=limarray_sp.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_diat=limarray_diat.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_diaz=limarray_diaz.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_cocco=limarray_cocco.argmin(dim='nutrient', skipna=False).squeeze()
mask = np.isnan(ds_year.sp_N_lim_surf.squeeze())
fig = plt.figure(figsize=(8,13))
colorbar_specs = {'ticks' : np.arange(0,4,1)}

ax = fig.add_subplot(4,1,2, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Diat surface nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_diat.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,1, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('SP surface nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_sp.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,3, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Cocco surface nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_cocco.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,4, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Diaz surface nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_diaz.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
cbar = fig.colorbar(pc, cax=cbar_ax,**colorbar_specs)
cbar.ax.set_yticklabels(['P lim', 'Fe lim', 'N lim','SiO3/C lim']);

Mapping biomass-weighted nutrient limitation in the top 100m

limarray_sp=xr.concat((ds_year.sp_P_lim_Cweight_avg_100m, ds_year.sp_Fe_lim_Cweight_avg_100m,ds_year.sp_N_lim_Cweight_avg_100m),dim='nutrient')
limarray_diat=xr.concat((ds_year.diat_P_lim_Cweight_avg_100m, ds_year.diat_Fe_lim_Cweight_avg_100m, ds_year.diat_N_lim_Cweight_avg_100m, ds_year.diat_SiO3_lim_Cweight_avg_100m),dim='nutrient')
limarray_diaz=xr.concat((ds_year.diaz_P_lim_Cweight_avg_100m, ds_year.diaz_Fe_lim_Cweight_avg_100m),dim='nutrient')
limarray_cocco=xr.concat((ds_year.cocco_P_lim_Cweight_avg_100m, ds_year.cocco_Fe_lim_Cweight_avg_100m, ds_year.cocco_N_lim_Cweight_avg_100m, ds_year.cocco_C_lim_Cweight_avg_100m),dim='nutrient')
                          
most_lim_sp=limarray_sp.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_diat=limarray_diat.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_diaz=limarray_diaz.argmin(dim='nutrient', skipna=False).squeeze()
most_lim_cocco=limarray_cocco.argmin(dim='nutrient', skipna=False).squeeze()
mask = np.isnan(ds_year.sp_N_lim_surf.squeeze())
fig = plt.figure(figsize=(8,13))
colorbar_specs = {'ticks' : np.arange(0,4,1)}

ax = fig.add_subplot(4,1,2, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Diat biomass-weighted nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_diat.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,1, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('SP biomass-weighted nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_sp.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,3, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Cocco biomass-weighted nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_cocco.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax = fig.add_subplot(4,1,4, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('Diaz biomass-weighted nutrient limitation', fontsize=12)
lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_diaz.where(~mask))
pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
cbar = fig.colorbar(pc, cax=cbar_ax,**colorbar_specs)
cbar.ax.set_yticklabels(['P lim', 'Fe lim', 'N lim','SiO3/C lim']);

Making monthly climatology maps

Make a monthly climatology dataset

A monthly climatology is a dataset where data from each month, including over different years, is averaged together. So for our dataset, the groups would include the average of all Januaries 2010-2019, all Februaries 2010-2019, and so on. This would usually be over a longer time period such as 30 or more years, but we use our shorter dataset so that it processes faster. This technique is helpful for looking at seasonal phenomena while averaging out year-to-year fluctuations.

mon_ds = ds.copy()
mon_ds = ds.groupby('time.month').mean('time')
limarray_sp=xr.concat((mon_ds.sp_P_lim_surf, mon_ds.sp_Fe_lim_surf,mon_ds.sp_N_lim_surf),dim='nutrient')
limarray_diat=xr.concat((mon_ds.diat_P_lim_surf, mon_ds.diat_Fe_lim_surf, mon_ds.diat_N_lim_surf, mon_ds.diat_SiO3_lim_surf),dim='nutrient')
limarray_diaz=xr.concat((mon_ds.diaz_P_lim_surf, mon_ds.diaz_Fe_lim_surf),dim='nutrient')
limarray_cocco=xr.concat((mon_ds.cocco_P_lim_surf, mon_ds.cocco_Fe_lim_surf, mon_ds.cocco_N_lim_surf, mon_ds.cocco_C_lim_surf),dim='nutrient')

most_lim_sp=limarray_sp.argmin(dim='nutrient', skipna=False)
most_lim_diat=limarray_diat.argmin(dim='nutrient', skipna=False)
most_lim_diaz=limarray_diaz.argmin(dim='nutrient', skipna=False)
most_lim_cocco=limarray_cocco.argmin(dim='nutrient', skipna=False)
mask = np.isnan(ds_year.sp_N_lim_surf.squeeze())
fig = plt.figure(figsize=(12,23))

month_list = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]

for row in np.arange(1,13):
    ts=row-1
    
    plot = row*3 - 2
    
    #row 1 Jan
    ax = fig.add_subplot(12,3,plot, projection=ccrs.Robinson(central_longitude=305.0))
    ax.set_title('Diat surface nutrient limitation', fontsize=12)
    lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_diat.isel(month=ts).where(~mask))
    pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
    ax.add_feature(land)
    ax.set_ylabel(month_list[ts])
    ax.set_yticks([]) # necessary to get ylabel to show up
    colorbar_specs = {'ticks' : np.arange(0,4,1)}
    
    plot = row*3 - 1
    ax = fig.add_subplot(12,3,plot, projection=ccrs.Robinson(central_longitude=305.0))
    ax.set_title('SP surface nutrient limitation', fontsize=12)
    lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_sp.isel(month=ts).where(~mask))
    pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
    ax.add_feature(land)
    colorbar_specs = {'ticks' : np.arange(0,4,1)}
    
    plot = row*3
    ax = fig.add_subplot(12,3,plot, projection=ccrs.Robinson(central_longitude=305.0))
    ax.set_title('Cocco surface nutrient limitation', fontsize=12)
    lon, lat, field = adjust_pop_grid(lons, lats,  most_lim_cocco.isel(month=ts).where(~mask))
    pc=ax.pcolormesh(lon, lat, field, cmap=matplotlib.colormaps['viridis'].resampled(4),vmin=-0.5,vmax=3.5,transform=ccrs.PlateCarree())
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
    ax.add_feature(land)
    colorbar_specs = {'ticks' : np.arange(0,4,1)}

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
cbar = fig.colorbar(pc, cax=cbar_ax,**colorbar_specs)
cbar.ax.set_yticklabels(['P lim', 'Fe lim', 'N lim','SiO3/C lim']);

And close the Dask cluster we spun up at the beginning.

cluster.close()

Summary

You’ve learned how to evaluate phytoplankton nutrient limitation.

Resources and references

References
  1. (2013). In Ocean Biogeochemical Dynamics (pp. 102–172). Princeton University Press. 10.2307/j.ctt3fgxqx.7
Plankton
Phytoplankton biomass
Plankton
Zooplankton biomass