Ocean macronutrients

Overview

The availability of several macronutrients controls production in most of the ocean: nitrate, phosphate, and silicate. Here we take a look at maps and depth profiles of these nutrients, and compare them to an observational dataset.

General setup
Subsetting
Processing - means in time and space
Compare to World Ocean Atlas data
Make depth profiles

Prerequisites

Concepts	Importance	Notes
Matplotlib	Necessary
Intro to Cartopy	Necessary
Dask Cookbook	Helpful
Intro to Xarray	Helpful

Time to learn: 30 min

Imports

import xarray as xr
import glob
import numpy as np
import matplotlib.pyplot as plt
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()

cluster.scale(20)

client

Client

Client-e8cd09a5-550a-11ef-8a5b-6045bdbbe7f0

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

LocalCluster

4e1a915c

Dashboard: http://127.0.0.1:8787/status	Workers: 4
Total threads: 4	Total memory: 15.61 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-b6e94331-dba2-4241-80e9-aa36105e09ea

Comm: tcp://127.0.0.1:37105	Workers: 4
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 15.61 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:37653	Total threads: 1
Dashboard: http://127.0.0.1:37001/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:38187
Local directory: /tmp/dask-scratch-space/worker-rbw_hvph

Worker: 1

Comm: tcp://127.0.0.1:43765	Total threads: 1
Dashboard: http://127.0.0.1:44847/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:43741
Local directory: /tmp/dask-scratch-space/worker-3tiygbrr

Worker: 2

Comm: tcp://127.0.0.1:37639	Total threads: 1
Dashboard: http://127.0.0.1:34117/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:33443
Local directory: /tmp/dask-scratch-space/worker-c0zl8su_

Worker: 3

Comm: tcp://127.0.0.1:45899	Total threads: 1
Dashboard: http://127.0.0.1:45051/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:43045
Local directory: /tmp/dask-scratch-space/worker-2aurs0s7

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)

Subsetting

Make our dataset smaller so it has just a couple of macronutrient variables we’re interested in.

variables =['PO4','NO3','SiO3']

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])

Let’s take a quick look at nitrate to make sure that things look okay…

ds.NO3.isel(time=0,z_t=0).plot(cmap="viridis")

/home/runner/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/distributed/client.py:3362: UserWarning: Sending large graph of size 10.01 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

<matplotlib.collections.QuadMesh at 0x7f9e7c0d5070>

../_images/d605f0caefc082176d59817b8b9c993b3ba98508bcbbd80897eac0d20070453e.png

Transforming from monthly to annual data

We can’t just use xarray’s regular mean() function because months have different numbers of days in them, so we have to weight by that to ensure the annual mean is accurate. See this ESDS blog post for a more detailed explanation with examples!

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.year'
    weights = (
        month_length.groupby("time.year") / month_length.groupby("time.year").sum()
    )

    # Test that the sum of the year 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")

ds_annual = year_mean(ds)
ds_annual

<xarray.Dataset> Size: 2GB
Dimensions:   (year: 10, z_t: 60, nlat: 384, nlon: 320, z_t_150m: 15)
Coordinates:
    TLAT      (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    TLONG     (nlat, nlon) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
  * z_t       (z_t) float32 240B 500.0 1.5e+03 2.5e+03 ... 5.125e+05 5.375e+05
  * z_t_150m  (z_t_150m) float32 60B 500.0 1.5e+03 2.5e+03 ... 1.35e+04 1.45e+04
  * year      (year) int64 80B 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019
Dimensions without coordinates: nlat, nlon
Data variables:
    NO3       (year, z_t, nlat, nlon) float64 590MB dask.array<chunksize=(1, 15, 96, 80), meta=np.ndarray>
    PO4       (year, z_t, nlat, nlon) float64 590MB dask.array<chunksize=(1, 15, 96, 80), meta=np.ndarray>
    SiO3      (year, z_t, nlat, nlon) float64 590MB dask.array<chunksize=(1, 15, 96, 80), meta=np.ndarray>

Note that our time coordinate is now called year instead, and has only years now. We can select specific years to plot:

ds_annual['NO3'].sel(year=2010).isel(z_t=0).plot()

/home/runner/miniconda3/envs/ocean-bgc-cookbook-dev/lib/python3.12/site-packages/distributed/client.py:3362: UserWarning: Sending large graph of size 10.01 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

<matplotlib.collections.QuadMesh at 0x7f9e8c5e7950>

../_images/32e7e9fea1a15e6ac66e354e4b77e9b12e0cdbf38addf719d49859e5fa76bcca.png

Let’s make a nicer-looking map

fig = plt.figure(figsize=(8,6))

ax = fig.add_subplot(1,1,1, projection=ccrs.Robinson(central_longitude=305.0))
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

ax.set_title('CESM surface NO$_3$', fontsize=10)
lon, lat, field = adjust_pop_grid(lons, lats, ds_annual.NO3.sel(year=2010).isel(z_t=0))
pc1=ax.pcolormesh(lon, lat,field, vmin=0, vmax=20, cmap='Greens',
                 transform=ccrs.PlateCarree())

cbar1 = fig.colorbar(pc1, ax=ax,extend='max',label='NO$_3$ (mmol m$^{-3}$)')

../_images/9a7f244a3d5822c0265b014a7fe1555877266961dd31c1957f544d69360c9772.png

Global mean macronutrient profiles

Let’s write a function to take a global mean of the variables we’re interested in, so that we can look at some depth profiles rather than maps. Also remember that we already took a mean over the whole time range (and the WOA dataset already had this mean taken), so this is a mean in time as well. Like the above maps, we also plot a bias panel to directly compare the difference between the datasets.

def global_mean(ds, ds_grid, compute_vars, include_ms=False):
    """
    Compute the global mean on a POP dataset. 
    Return computed quantity in conventional units.
    """
    dso = xr.Dataset({v: ds_grid[v] for v in ['z_t']})
    
    for var in compute_vars:
        
        area_depth = np.full([384,320,60],np.nan)
        var_profile = np.full([60],np.nan)

        for z in np.arange(0,60,1):
            
            if include_ms: # marginal seas
                area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.KMT > 0).where(ds[var].isel(z_t=z) > 0)
            
            else: 
                area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)  
            
        area_depth = xr.DataArray(area_depth,dims=('nlat','nlon','z_t'))
            
        for z in np.arange(0,60,1):
            
            var_profile[z] = (ds[var].isel(z_t=z) * area_depth.isel(z_t=z)).sum(dim=('nlon','nlat')) / area_depth.isel(z_t=z).sum(dim=('nlon','nlat'))
    
        dso[var] = var_profile
                
    return dso

ds_glb = global_mean(ds_mean, ds_grid, ['NO3','PO4','SiO3']).compute()

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

ds_glb_woa = global_mean(ds_woa, ds_grid, ['NO3','PO4','SiO3']).compute()

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

/tmp/ipykernel_2651/680694341.py:19: DeprecationWarning: __array__ implementation doesn't accept a copy keyword, so passing copy=False failed. __array__ must implement 'dtype' and 'copy' keyword arguments.
  area_depth[:,:,z] = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).where(ds[var].isel(z_t=z) > 0)

Rather than setting up a dictionary of parameters, here we choose to make the plots inline since there aren’t as many.

fig = plt.figure(figsize=(6,10))

plt.suptitle('Global mean macronutrient profiles', fontsize=14)

### Row 1  - NO3

ax = fig.add_subplot(3,2,1)
ax.set_title('Global mean NO$_3$')
ax.plot(ds_glb_woa['NO3'].values, depths, label='WOA', linewidth=3, color='lightgreen')
ax.plot(ds_glb['NO3'].values, depths, label='CESM', linewidth=3, color='green')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='NO$_3$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()

# Bias plot

ax = fig.add_subplot(3,2,2)
ax.plot(ds_glb['NO3'].values - ds_glb_woa['NO3'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='NO$_3$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()

### Row 2  - PO4

ax = fig.add_subplot(3,2,3)
ax.set_title('Global mean PO$_4$')
ax.plot(ds_glb_woa['PO4'].values, depths, label='WOA', linewidth=3, color='peachpuff')
ax.plot(ds_glb['PO4'].values, depths, label='CESM', linewidth=3, color='orange')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='PO$_4$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()

# Bias plot

ax = fig.add_subplot(3,2,4)
ax.plot(ds_glb['PO4'].values - ds_glb_woa['PO4'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='PO$_4$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()

### Row 3  - SiO3

ax = fig.add_subplot(3,2,5)
ax.set_title('Global mean SiO$_3$')
ax.plot(ds_glb_woa['SiO3'].values, depths, label='WOA', linewidth=3, color='lightblue')
ax.plot(ds_glb['SiO3'].values, depths, label='CESM', linewidth=3, color='blue')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='SiO$_3$ (mmol m$^{-3}$)')
plt.gca().invert_yaxis()

# Bias plot

ax = fig.add_subplot(3,2,6)
ax.plot(ds_glb['SiO3'].values - ds_glb_woa['SiO3'].values, depths, label='bias', linewidth=3, color='black')
ax.legend()
ax.set(ylabel='depth (m)',xlabel='SiO$_3$ bias (mmol m$^{-3}$)')
ax.axvline(x=0,color='black',linestyle='--',linewidth=0.5)
plt.gca().invert_yaxis()

fig.tight_layout()

../_images/62694116e54a73063f550beb6f06311ec3afd5032cfa79a79a0b0bdbd4975eb1.png

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

cluster.close()

2024-08-07 22:20:01,054 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,055 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,055 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,056 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,057 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,058 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,059 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,061 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,062 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,065 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,066 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,073 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,077 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,080 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,081 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,085 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,093 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,093 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,101 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

2024-08-07 22:20:01,112 - distributed.nanny - WARNING - Worker process still alive after 4.0 seconds, killing

Summary

You’ve learned how to plot and evaluate the distribution of some key ocean nutrients in CESM output.