Analyzing Carbon Dioxide Removal (CDR) Metrics

Analyzing Carbon Dioxide Removal (CDR) Metrics#

This notebook performs a Carbon Dioxide Removal (CDR) analysis using C-Star model output and grid data. We

compute CDR efficiency and CO\(_2\)-uptake curves
visualize the results
and save a dedicated CDR dataset to NetCDF for further analysis or verification.

[1]:

from roms_tools import Grid, ROMSOutput

/global/homes/u/uheede/.conda/envs/roms-tools-env/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

First, we read the grid data from a file. For more details, see this page on creating a grid from an existing NetCDF (or YAML) file.

[2]:

grid = Grid(filename="/global/cfs/projectdirs/m4746/Datasets/ROMSOutput/gom12km/grid_12km.nc")

Next, we read the C-Star model output. For more details, see this page on reading ROMS output data.

[3]:

roms_output = ROMSOutput(
    grid=grid,
    path=[
        "/global/cfs/projectdirs/m4746/Datasets/ROMSOutput/gom12km/atlas_cdr1_cstar.nc"
    ],
    use_dask=True,
)

[4]:

roms_output.ds

[4]:

<xarray.Dataset> Size: 32GB
Dimensions:         (time: 182, eta_rho: 130, xi_rho: 210, s_rho: 50)
Coordinates:
  * time            (time) datetime64[ns] 1kB 2000-01-02 ... 2000-07-01
    lat_rho         (eta_rho, xi_rho) float64 218kB ...
    lon_rho         (eta_rho, xi_rho) float64 218kB ...
Dimensions without coordinates: eta_rho, xi_rho, s_rho
Data variables: (12/23)
    ocean_time      (time) float64 1kB dask.array<chunksize=(1,), meta=np.ndarray>
    zeta            (time, eta_rho, xi_rho) float64 40MB dask.array<chunksize=(1, 130, 210), meta=np.ndarray>
    temp            (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    salt            (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    ALK             (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    DIC             (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    ...              ...
    ALK_ALT_source  (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    DIC_source      (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    DIC_ALT_source  (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    avg_begin_time  (time) float64 1kB dask.array<chunksize=(1,), meta=np.ndarray>
    avg_end_time    (time) float64 1kB dask.array<chunksize=(1,), meta=np.ndarray>
    mask_rho        (eta_rho, xi_rho) int32 109kB ...
Attributes: (12/44)
    global_x:           208
    global_y:           128
    title:              GoM Atlas CDR runs
    grid_file:          ../input_files/partitioned_files/grid_12km.000.nc
    forcing_files:       ../input_files/partitioned_files/surf_frc_phys_20000...
    forcing_info:
    ...                 ...
    pipe_frc_options:   OFF
    particle_options:   OFF
    git_version:        3e02f2d1a2e10b3fb848e5197ea9b360269fb44a
    averaging:          All variables are averaged in time
    history:            Tue Oct  7 16:32:36 2025: ncrcat ../atlas_cdr1_cstar....
    NCO:                netCDF Operators version 4.9.3 (Homepage = http://nco...

In this experiment, an alkalinity source has been added into the Gulf of Mexico off the coast of the Yucatán Peninsula. This can be visualized by plotting the uppermost layer of ALK_source, a variable from the C-Star output. For more details on visualizing ROMS output fields, see this page.

[5]:

roms_output.plot("ALK_source", s=-1, time=7)

[########################################] | 100% Completed | 101.07 ms

In the next figure, we plot the time series of the total alkalinity source flux.

[6]:

import matplotlib.pyplot as plt

[7]:

def plot_alk_flux_series(ds):
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(ds.time, ds.ALK_source.sum(dim=["eta_rho", "xi_rho", "s_rho"]))
    ax.set(title="Total alkalinity source flux", xlabel="time", ylabel=ds.ALK_source.attrs["units"])
    ax.grid()

[8]:

fig = plot_alk_flux_series(roms_output.ds)

The figure indicates that the release has taken place over a short period in early January 2000.

CO₂ uptake and CDR uptake efficiency#

CDR metrics estimate how much extra CO₂ the ocean takes up because of the intervention, and how efficient that uptake is. The CO₂ uptake is reported as mass in tonnes of CO₂; the efficiency is a dimensionless ratio of moles of DIC transfered from the atmosphere to the ocean per moles of added alkalinity. The same physical story is expressed in two parallel ways (see also here).

Flux-based CO₂ uptake. The first method uses the difference between the air–sea CO₂ flux in the intervention run and in a counter-factual (e.g. alternative CO₂) run. That difference is integrated over the model sea surface and accumulated in time. The implied change in ocean carbon is converted to tonnes of CO₂ using the molar mass of CO₂.

Inventory-based CO₂ uptake. The second method uses the difference in dissolved inorganic carbon between the same two simulations, after accounting for layer thickness and cell geometry so the field represents an inventory per column or per volume element. That difference is integrated over the model domain at each time. Again the result is expressed as tonnes of CO₂ using the same mass conversion.

Both integrals cover the interior of the grid. The outermost row and column of rho points on each side are boundary points, whose values the boundary condition sets, so they are left out.

CDR uptake efficiency. Each uptake construction is normalised by the cumulative source term from the intervention—the added alkalinity (written below as alkalinity source). That yields two efficiency curves over time.

Efficiency from CO₂ flux differences: \begin{equation} \eta_\text{flux}(t_f) = \frac{ \displaystyle \int_{t_0}^{t_f} \int \big( \Phi_\text{intervention} - \Phi_\text{baseline} \big) \, dA \, dt }{ \displaystyle \int_{t_0}^{t_f} \int \text{ALK}_\text{source} \, dV \, dt } \end{equation}

Efficiency from DIC differences: \begin{equation} \eta_\text{DIC}(t_f) = \frac{ \displaystyle \int \big( \text{DIC}_\text{intervention} - \text{DIC}_\text{baseline} \big) \, dV }{ \displaystyle \int_{t_0}^{t_f} \int \text{ALK}_\text{source} \, dV \, dt } \end{equation}

Symbol	Meaning	Units
\(t_0\), \(t_f\)	Start and end times	s
\(\Phi_\text{intervention}\), \(\Phi_\mathrm{baseline}\)	Air–sea CO₂ flux, intervention vs baseline	mmol/m²/s
\(\text{DIC}_\text{intervention}\), \(\text{DIC}_\text{baseline}\)	DIC fields for the two runs evaluated at \(t_f\)	mmol/m³
\(\text{ALK}_\text{source}\)	Added alkalinity from the intervention	meq/s
\(dA\), \(dV\)	Surface area and volume elements	m², m³

Together, the two methods give two estimates of CO₂ uptake (tonnes CO₂) and two estimates of CDR uptake efficiency; comparing them checks consistency between the flux and inventory views of the same CDR signal. We can perform a CDR analysis, including computing the CDR uptake efficiency and CO\(_2\) uptake in both ways, by calling the method .cdr_metrics(). This method will automatically compute the required diagnostics (if not already present), save them in the attribute .ds_cdr, and produce a plot of the uptake efficiency over time.

[ ]:

roms_output.cdr_metrics()

2026-04-23 12:20:47 - INFO - Max absolute difference between flux-based and DIC-based uptake efficiency: 1.375e-02

A note on discrete integrals, and differences between methods:

Because the numerical model outputs data at discrete intervals, the analytical continuous integrals must be evaluated using finite summations, which can cause discrepancies between the flux-based and DIC inventory based methods.

Flux in the above equations represents the time-averaged gas flux output across the interval \((t_{i-1}, t_i]\), where DIC represents the instantaneous, 3D dissolved inorganic carbon concentrations at the final evaluation time-step \(t_f\). Because the flux-based method involves time averaging from \(t_{i-1}\) to \(t_i\), it lags the DIC inventory method by half a time-step.

CDR-relevant metrics#

After running .cdr_metrics(), a new attribute ds_cdr is created, containing all CDR-relevant outputs.

[10]:

roms_output.ds_cdr

[10]:

<xarray.Dataset> Size: 4GB
Dimensions:                                (time: 182, eta_rho: 130,
                                            xi_rho: 210, s_rho: 50)
Coordinates:
  * time                                   (time) datetime64[ns] 1kB 2000-01-...
    lat_rho                                (eta_rho, xi_rho) float64 218kB 17...
    lon_rho                                (eta_rho, xi_rho) float64 218kB 26...
Dimensions without coordinates: eta_rho, xi_rho, s_rho
Data variables:
    FG_CO2                                 (time, eta_rho, xi_rho) float64 40MB dask.array<chunksize=(1, 130, 210), meta=np.ndarray>
    FG_ALT_CO2                             (time, eta_rho, xi_rho) float64 40MB dask.array<chunksize=(1, 130, 210), meta=np.ndarray>
    hDIC                                   (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    hDIC_ALT_CO2                           (time, s_rho, eta_rho, xi_rho) float64 2GB dask.array<chunksize=(1, 50, 130, 210), meta=np.ndarray>
    area                                   (eta_rho, xi_rho) float64 218kB 1....
    window_length                          (time) float64 1kB dask.array<chunksize=(1,), meta=np.ndarray>
    cdr_efficiency_from_flux               (time) float64 1kB nan nan ... 0.705
    cdr_efficiency_from_DIC_difference     (time) float64 1kB nan nan ... 0.7045
    cdr_carbon_uptake_from_flux            (time) float64 1kB 0.0 0.0 ... 4.825
    cdr_carbon_uptake_from_DIC_difference  (time) float64 1kB 0.0 0.0 ... 4.822

We can save the CDR-relevant outputs (stored in .ds_cdr) to a NetCDF file using xarray’s .to_netcdf() method.

[11]:

cdr_filepath = "/pscratch/sd/u/uheede/output/cdr_metrics.nc"

[12]:

roms_output.ds_cdr.to_netcdf(cdr_filepath)

Analyzing Carbon Dioxide Removal (CDR) Metrics

Contents

Analyzing Carbon Dioxide Removal (CDR) Metrics#

CO₂ uptake and CDR uptake efficiency#

CDR-relevant metrics#

Ensembles of CDR efficiency curves#

Nested (upscaled) experiments#