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/romstools-test/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
_images/cdr_analysis_8_1.png

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)
_images/cdr_analysis_12_0.png

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.

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
_images/cdr_analysis_15_1.png

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)

Ensembles of CDR efficiency curves#

Let’s now look at ensembles. For this we have run a second CDR experiment. Let’s open its output!

[13]:

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

The release location is the same as in the first CDR experiment.

[14]:

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

[########################################] | 100% Completed | 101.21 ms
_images/cdr_analysis_27_1.png

However, plotting the total alkalinity source flux indicates that the release occurred over a short period in early April 2000, whereas in the first experiment it took place in early January 2000.

[15]:

fig = plot_alk_flux_series(another_roms_output.ds)
_images/cdr_analysis_29_0.png

Note also that the second experiment only ran for three months, while the first ran for six months.

[16]:

another_roms_output.cdr_metrics()

2026-04-23 12:22:52 - INFO - Max absolute difference between flux-based and DIC-based uptake efficiency: 1.088e-02
_images/cdr_analysis_31_1.png

We can now construct a small ensemble consisting of two members: our two experiments. Note that ensemble members do not need to share the same start or end dates.

[17]:

from roms_tools import Ensemble

The members variable is a dictionary whose values can be either:

  • file paths (strings) pointing to the CDR metrics, or

  • xarray.Dataset objects containing the metrics directly.

Let’s use a mix of both!

[18]:

ensemble = Ensemble(
    members={
        "January_release": cdr_filepath,
        "April_release": another_roms_output.ds_cdr
    }
)

The .ds attribute contains an xarray.Dataset that aggregates all releases, aligns them relative to the release start date, and computes the ensemble mean and standard deviation for both efficiency and CO\(_2\) uptake.

[19]:

ensemble.ds

[19]:

<xarray.Dataset> Size: 13kB
Dimensions:                     (time: 182)
Coordinates:
  * time                        (time) timedelta64[s] 1kB -7 days ... 174 days
Data variables:
    January_release             (time) float64 1kB nan nan nan ... 0.7043 0.705
    April_release               (time) float64 1kB nan nan nan ... nan nan nan
    ensemble_efficiency_mean    (time) float64 1kB nan nan nan ... 0.7043 0.705
    ensemble_efficiency_std     (time) float64 1kB nan nan nan ... 0.0 0.0 0.0
    January_release_co2_uptake  (time) float64 1kB 0.0 0.0 0.0 ... 4.82 4.825
    April_release_co2_uptake    (time) float64 1kB 0.0 0.0 0.0 ... nan nan nan
    ensemble_uptake_mean        (time) float64 1kB 0.0 0.0 0.0 ... 4.82 4.825
    ensemble_uptake_std         (time) float64 1kB 0.0 0.0 0.0 ... 0.0 0.0 0.0

We can make a plot and save it as follows.

[ ]:

ensemble.plot(save_path="figures/ensemble.png")
_images/cdr_analysis_39_0.png 
[ ]: