Ocean carbon fluxes - Ocean Biogeochemistry Cookbook

Overview¶

The carbon cycle is a key part of ocean biogeochemistry and, more broadly, Earth’s climate system. Here we learn how to make maps of some key variables modeled by CESM related to the marine carbon cycle.

General setup
Subsetting
Processing data
Making maps

Prerequisites¶

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

Time to learn: 15 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 dask
import distributed
import s3fs

from module import adjust_pop_grid

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 7
      3 import numpy as np
      4 import matplotlib.pyplot as plt
      5 import cartopy
      6 import cartopy.crs as ccrs
----> 7 import pop_tools
      8 from dask.distributed import LocalCluster
      9 import dask
     10 import distributed

File ~/micromamba/envs/ocean-bgc-cookbook-dev/lib/python3.14/site-packages/pop_tools/__init__.py:4
      1 """Top-level module for pop_tools"""
----> 4 from pkg_resources import DistributionNotFound, get_distribution
      6 from .calc import cfc11sol, cfc12sol
      7 from .datasets import DATASETS

ModuleNotFoundError: No module named 'pkg_resources'

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

ds_grid

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)
s3.invalidate_cache()
# 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 =['FG_CO2','photoC_TOT_zint','POC_FLUX_100m']
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 - means in time and space¶

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.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")

We also define a new function to take global mean in space.

def global_mean(ds, ds_grid, compute_vars, normalize=True, include_ms=False):
    """
    Compute the global mean on a POP dataset. 
    Return computed quantity in conventional units.
    """

    other_vars = list(set(ds.variables) - set(compute_vars))

    # note TAREA is in cm^2, which affects units
    
    if include_ms: # marginal seas!
        surface_mask = ds_grid.TAREA.where(ds_grid.KMT > 0).fillna(0.)
    else:
        surface_mask = ds_grid.TAREA.where(ds_grid.REGION_MASK > 0).fillna(0.)        
    
    masked_area = {
        v: surface_mask.where(ds[v].notnull()).fillna(0.) 
        for v in compute_vars
    }
    
    with xr.set_options(keep_attrs=True):
        
        dso = xr.Dataset({
            v: (ds[v] * masked_area[v]).sum(['nlat', 'nlon'])
            for v in compute_vars
        })
        
        if normalize:
            dso = xr.Dataset({
                v: dso[v] / masked_area[v].sum(['nlat', 'nlon'])
                for v in compute_vars
            })            
                
    return dso

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_mean(ds).mean("year")

Do some global integrals, to check if our values look reasonable¶

ds_glb = global_mean(ds, ds_grid, variables,normalize=False).compute()

# convert from nmol C/s to Pg C/yr
nmols_to_PgCyr = 1e-9 * 12. * 1e-15 * 365. * 86400.

for v in variables:
    ds_glb[v] = ds_glb[v] * nmols_to_PgCyr        
    ds_glb[v].attrs['units'] = 'Pg C yr$^{-1}$'
    
ds_glb

We can compare these values to some observationally derived values. Each of these is calculated in a different way with combinations of data and models--please reference each linked paper for detailed discussion. Takahashi et al., 2002 estimate global air-sea CO $_2$ flux to be 2.2 (+22% or −19%) Pg C yr $^{−1}$ . Our value (shown above as FG_CO2) is 2.779 Pg C yr $^{−1}$ . This value is outside of these bounds, but still on the same order of magnitude. We note that these values are calculated over different time periods, so we also don’t expect them to be an exact comparison. photoC_TOT_zint represents global vertically-integrated NPP; Behrenfeld and Falkowski, 1997 estimate this value to be 43.5 Pg C yr $^{−1}$ . Our value is 53.26 Pg C yr $^{−1}$ , which is within 22% of the observationally derived value. POC_FLUX_100m represents the particulate organic carbon flux at 100 m depth. DeVries and Weber, 2017 calculated this flux integrated over the entire euphotic zone to be 9.1 ± 0.2 Pg C yr $^{−1}$ . Since the depth ranges are different, this isn’t an exact comparison, but the orders of magnitude are similar. This first-pass analysis tells us that CESM is on the right track for these values.

Make some maps¶

First, convert from mmol/m3 cm/s to mmol/m2/day.

for var in variables:
    ds[var] = ds[var] * 0.01 * 86400.

Then, make a few maps of key carbon-related variables.

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

ax = fig.add_subplot(3,1,1, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('a) Air-sea CO$_2$ flux', fontsize=12,loc='left')
lon, lat, field = adjust_pop_grid(lons, lats,  ds.FG_CO2)
pc=ax.pcolormesh(lon, lat, field, cmap='bwr',vmin=-5,vmax=5,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

cbar1 = fig.colorbar(pc, ax=ax,extend='both',label='mmol m$^{-2}$ d$^{-1}$')


ax = fig.add_subplot(3,1,2, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('b) NPP', fontsize=12,loc='left')
lon, lat, field = adjust_pop_grid(lons, lats,  ds.photoC_TOT_zint)
pc=ax.pcolormesh(lon, lat, field, cmap='Greens',vmin=0,vmax=100,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

cbar1 = fig.colorbar(pc, ax=ax,extend='max',label='mmol m$^{-2}$ d$^{-1}$')

ax = fig.add_subplot(3,1,3, projection=ccrs.Robinson(central_longitude=305.0))
ax.set_title('c) POC flux at 100m', fontsize=12,loc='left')
lon, lat, field = adjust_pop_grid(lons, lats,  ds.POC_FLUX_100m)
pc=ax.pcolormesh(lon, lat, field, cmap='Oranges',vmin=0,vmax=10,transform=ccrs.PlateCarree())
land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale='110m', edgecolor='k', facecolor='white', linewidth=0.5)
ax.add_feature(land)

cbar1 = fig.colorbar(pc, ax=ax,extend='max',label='mmol m$^{-2}$ d$^{-1}$');

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

cluster.close()

Summary¶

You’ve learned how to make maps of some key quantities related to oceanic carbon.

Resources and references¶

References¶

Takahashi, T., Sutherland, S. C., Sweeney, C., Poisson, A., Metzl, N., Tilbrook, B., Bates, N., Wanninkhof, R., Feely, R. A., Sabine, C., Olafsson, J., & Nojiri, Y. (2002). Global sea–air CO2 flux based on climatological surface ocean pCO2, and seasonal biological and temperature effects. Deep Sea Research Part II: Topical Studies in Oceanography, 49(9–10), 1601–1622. 10.1016/s0967-0645(02)00003-6
Behrenfeld, M. J., & Falkowski, P. G. (1997). Photosynthetic rates derived from satellite‐based chlorophyll concentration. Limnology and Oceanography, 42(1), 1–20. 10.4319/lo.1997.42.1.0001
DeVries, T., & Weber, T. (2017). The export and fate of organic matter in the ocean: New constraints from combining satellite and oceanographic tracer observations. Global Biogeochemical Cycles, 31(3), 535–555. 10.1002/2016gb005551
(2013). In Ocean Biogeochemical Dynamics (pp. 318–358). Princeton University Press. 10.2307/j.ctt3fgxqx.12