Cycling DART–CESM Data Assimilation with Regional MOM6

Run a real forecast–assimilate–update cycle: a 3-member regional MOM6 ensemble in CESM, corrected every six hours by the observations you made in Tutorial 1.

The workflow consists of six main steps:

Install the CESM DA fork.
(Re)generate the model domain.
Create a multi-instance (ensemble) CESM case.
Configure the case for assimilation.
Stage your observations and submit.
Examine what the assimilation did.

This is Part 3 of the DART tutorial series: 1. Creating Observations · 2. Synthetic Observations · 3. Cycling DART–CESM

SECTION 1: Install the CESM DA fork¶

DART is integrated as the ESP component of CESM in the CROCODILE fork: it is built and called by CESM like any other component, so a data-assimilation experiment is submitted with plain ./case.submit. Clone the fork, check out the DA branch, and populate the sub-components:

git clone https://github.com/hkershaw-brown/CESM.git CESM_DA
cd CESM_DA/
git checkout dart-cesm3.0-alphabranch
./bin/git-fleximod update

You will point cesmroot at this clone in Section 3. Your CrocoDash conda environment from the CrocoDash tutorial is the only other software you need.

SECTION 2: The model domain¶

Step 2.1: Parameters¶

The shared series parameters — identical to Tutorials 1 and 2.

# --- CROCODILE DART tutorial series parameters (same cell in all 3 notebooks) ---
from pathlib import Path
import datetime

WORKDIR = Path("<DART_WORKDIR>")        # your scratch working directory on Derecho

START = datetime.datetime(2013, 4, 1)   # must match RUN_STARTDATE in Tutorial 3
END   = datetime.datetime(2013, 4, 4)   # 3 days -> 12 six-hour assimilation windows
FREQ  = datetime.timedelta(hours=6)

# Bounding box: the panama1 domain (lon 278-281E, lat 7-10N) padded by ~2 degrees
LAT_MIN, LAT_MAX = 5.0, 12.0
LON_MIN, LON_MAX = -84.0, -77.0

OBS_TYPES = ["ARGO_TEMPERATURE", "ARGO_SALINITY"]

REAL_OBS_DIR      = WORKDIR / "obs" / "real"       # Tutorial 1 output
SYNTHETIC_OBS_DIR = WORKDIR / "obs" / "synthetic"  # Tutorial 2 output

Step 2.2: Regenerate the panama1 domain¶

Domain generation is taught in the CrocoDash tutorial; here we just re-run the three cells that rebuild the same panama1 domain, because creating a case needs the live grid, topography, and vertical-grid objects.

from CrocoDash.grid import Grid

grid = Grid(
    resolution=0.05,  # degrees
    xstart=278.0,     # min longitude, in [0, 360]
    lenx=3.0,         # longitude extent, degrees
    ystart=7.0,       # min latitude
    leny=3.0,         # latitude extent, degrees
    name="panama1",
)

from CrocoDash.topo import Topo

topo = Topo(grid=grid, min_depth=9.5)

bathymetry_path = Path("<GEBCO>")
topo.set_from_dataset(
    bathymetry_path=bathymetry_path,
    longitude_coordinate_name="lon",
    latitude_coordinate_name="lat",
    vertical_coordinate_name="elevation",
)

from CrocoDash.vgrid import VGrid

vgrid = VGrid.hyperbolic(
    nk=75,                  # number of vertical levels
    depth=topo.max_depth,
    ratio=20.0,             # ratio of top to bottom layer thicknesses
)

SECTION 3: Create a multi-instance CESM case¶

Step 3.1: Case name and directories¶

One new idea compared to the CrocoDash tutorial: ninst. An ensemble filter needs an ensemble — CESM runs ninst copies (“instances”) of the ocean, each from slightly different initial conditions, and DART updates all of them at every cycle.

casename = "panama-da"

cesmroot = "<DART_CESMROOT>"                 # your CESM_DA clone from Section 1
inputdir = WORKDIR / "croc_input" / casename # where input files are written
caseroot = WORKDIR / "croc_cases" / casename # the CESM case directory

Step 3.2: Create the case¶

Three members is workshop-sized — small enough to build and run in a session. Real ocean DA experiments use 30–80 members; with only 3, the sampling noise in the ensemble covariances is large, which is exactly why inflation and localization (Section 4) exist.

from CrocoDash.case import Case

case = Case(
    cesmroot=cesmroot,
    caseroot=caseroot,
    inputdir=inputdir,
    ocn_grid=grid,
    ocn_vgrid=vgrid,
    ocn_topo=topo,
    project="<PROJECT_CODE>",
    override=True,
    machine="derecho",
    compset="GR_JRA",
    ninst=3,             # ensemble size: 3 instances of MOM6
)

Alternative: raw create_newcase

If you are not using CrocoDash (for example, on a pre-existing global grid), the same case can be created directly with CIME. G_JRA_DA is the DA-enabled ocean compset; --multi-driver runs all members in a single job:

./cime/scripts/create_newcase \
    --run-unsupported \
    --res TL319_t232 \
    --compset G_JRA_DA \
    --case $casedir \
    --ninst 3 \
    --multi-driver \
    --project <PROJECT_CODE>

Answer

Everything. The ensemble spread is the filter’s estimate of forecast uncertainty, and the ensemble covariance between an observed quantity and the model state is what turns an observation of temperature at one point into corrections of salinity, currents, and temperature nearby. With one member there is no spread and no covariance — no update. The ensemble is the error model.

3 is a very small ensemble, so the covariance estimates are noisy. In practice, 30–80 members is common, and the ensemble size is a key control on the quality of the assimilation.

Step 3.3: Prepare forcing data¶

Exactly as in the CrocoDash tutorial: an initial condition plus one time-dependent segment per open boundary, cut from the GLORYS reanalysis, covering the experiment period.

case.configure_forcings(
    date_range=[START.strftime("%Y-%m-%d 00:00:00"), END.strftime("%Y-%m-%d 00:00:00")],
    boundaries=["south", "west"],   # the open (non-land) boundaries of panama1
    function_name="get_glorys_data_from_rda",
)

case.process_forcings()

SECTION 4: Configure the case for assimilation¶

Step 4.1: Turn on data assimilation¶

Three XML settings turn a regional ocean case into a cycling DA experiment. Run these in a terminal in your case directory (caseroot above):

cd <DART_WORKDIR>/croc_cases/panama-da
./case.setup

./xmlchange CALENDAR=GREGORIAN
./xmlchange DATA_ASSIMILATION_OCN=TRUE
./xmlchange RUN_STARTDATE=2013-04-01

# 12 six-hour cycles: 2013-04-01 00:00 -> 2013-04-04 00:00
./xmlchange STOP_OPTION=nhours,STOP_N=6
./xmlchange DATA_ASSIMILATION_CYCLES=12

./case.setup --reset

To confirm DART is wired in as the ESP component:

./xmlquery --partial DATA_ASS   # DATA_ASSIMILATION_* flags
./xmlquery --partial ESP        # ESP component should be DART

Step 4.2: Build — and understand the cycle while you wait¶

qcmd -- ./case.build

The build takes a while. Perfect time to walk the loop your job will execute 12 times:

            ┌──────────────────────────────────────────────────────────┐
            ▼                                                          │
 1. FORECAST      all 3 MOM6 instances advance 6 hours                 │
 2. FILTER        CESM's ESP layer calls DART filter:                  │
                    reads obs_seq.<window>.out + all 3 model states    │
                    computes the ensemble update                       │
                    writes diagnostics (obs_seq.final, *assim_mean.nc) │
 3. UPDATE        updated restart files replace the forecast restarts  │
 4. ADVANCE       CESM resubmits the next 6-hour segment ──────────────┘

Step 4.3: A tour of the DART namelists¶

After setup, CESM writes a fully-resolved DART namelist to Buildconf/dartconf/input.nml (read-only — regenerated every build and submission), and the list of observation files DART expects to Buildconf/dart.input_data_list. To change DART settings, edit user_nl_dart in the case directory, exactly like user_nl_mom for MOM6.

The three groups that matter most for ocean DA:

&filter_nml — the ensemble filter itself:

Setting	Meaning
`inf_flavor`, `inf_initial`	inflation — grows ensemble spread to counter the overconfidence of small ensembles
`cutoff`	localization half-width in radians — limits how far one observation reaches. 0.02 rad ≈ 127 km at the equator

&model_nml — which MOM6 variables are in the state vector (temperature, salinity, SSH, velocities). Only state-vector variables are updated by the filter.

&obs_kind_nml — which observation types are assimilate_these_obs_types (they change the state) vs. evaluate_these_obs_types (forward operator computed and recorded, but no influence — a free out-of-sample check).

Answer

Salinity observations still get forward-operator values recorded in obs_seq.final — you can still compute their RMSE — but they no longer update the state; only temperature does. Evaluate mode is the safe way to vet a new observation type: check its statistics against the model for a few cycles, look for biases or gross errors, then promote it to assimilate.

SECTION 5: Stage the observations¶

Step 5.1: Put the obs where DART looks¶

Copy your Tutorial 1 files into the case’s run directory:

cd <DART_WORKDIR>/croc_cases/panama-da
RUNDIR=$(./xmlquery --value RUNDIR)
cp <DART_WORKDIR>/obs/real/obs_seq.*.out $RUNDIR/

Check Buildconf/dart.input_data_list if you want to see exactly which file names DART expects — they follow the obs_seq.YYYY-MM-DD-SSSSS.out convention from Tutorial 1.

A quick completeness check — 12 cycles need 12 windows of observations (windows with no observations are allowed; DART simply has nothing to assimilate that cycle):

obs_files = sorted(REAL_OBS_DIR.glob("obs_seq.*.out"))
print(f"{len(obs_files)} obs_seq file(s) staged for 12 cycles:")
for f in obs_files:
    print("  ", f.name)

SECTION 6: Submit and cycle¶

Step 6.1: Submit¶

./case.submit

That’s it — CESM runs the forecast–filter–update loop 12 times. To watch it:

qstat -u $USER                  # queue status
tail CaseStatus                 # case-level progress, cycle by cycle
ls $RUNDIR/*.log.*              # component logs
cat $RUNDIR/dart_log.out        # DART's own log: obs counts, rejections, inflation

Step 6.2: What appears while it runs¶

Each cycle, filter drops diagnostics into the run directory (file names carry the case name and the cycle timestamp):

File	Contents
`obs_seq.final`	Every observation with its prior (and posterior) ensemble estimates, and whether it was used
`preassim_mean.nc`	Ensemble-mean ocean state after inflation but before assimilation
`postassim_mean.nc`	Ensemble-mean ocean state after assimilation
`output_mean.nc`	Updated ensemble mean, restart form
`output_sd.nc`	Updated ensemble spread, restart form

SECTION 7: Examine the assimilation¶

Step 7.1: Observation-space diagnostics¶

obs_seq.final contains information about which observations were used in assimilation, and why any observations were not used. It also contains the prior and optionally the posterior forward operator values for each observation and their mean and spread.

pyDARTdiags can be used to calculate and plot observation-space diagnostics from obs_seq.final.

Look at used vs. rejected observations, and the prior and posterior statistics for each observation type.

import pydartdiags.obs_sequence.obs_sequence as obsq
from pydartdiags.stats import stats

# Your case run directory: cd caseroot && ./xmlquery --value RUNDIR
RUN_DIR = Path("<RUN_DIR>")

# File names carry the case name and cycle timestamp: ls $RUNDIR/*obs_seq*
obs_seq_final = obsq.ObsSequence(str(RUN_DIR / "obs_seq.final"))

used_obs = obs_seq_final.select_used_qcs()
stats.diag_stats(used_obs)
stats.grand_statistics(used_obs)

Step 7.2: State-space increments¶

The increment — posterior minus prior ensemble — is the map of what the observations changed.

import xarray as xr
import matplotlib.pyplot as plt

preassim  = xr.open_dataset(RUN_DIR / "preassim_mean.nc")
postassim = xr.open_dataset(RUN_DIR / "postassim_mean.nc")

# Sea-surface temperature increment (top model level)
increment = postassim["Temp"].isel(Time=0, zl=0) - preassim["Temp"].isel(Time=0, zl=0)

fig, ax = plt.subplots(figsize=(6, 5))
increment.plot(ax=ax, cmap="RdBu_r", robust=True)
ax.set_title("SST increment (posterior − prior)")
plt.show()

Answer

There should almost always be an observation at or near a strong increment maximum.
Yes — the ensemble covariance between temperature and salinity carries the update across variables. That cross-variable transfer is the whole power (and risk) of ensemble DA with a small ensemble.
Increment footprints shrink toward the observation locations; corrections far from any observation disappear. Too small a cutoff wastes information, too large lets 3-member sampling noise correct the far field with garbage.

Recap¶

Your takeaway: a DA-enabled case you can rerun and reconfigure, plus obs_seq.final and increment maps from your own 12-cycle experiment.

Where to go from here?¶

Run longer: set END = 2013-04-08 in Tutorial 1, regenerate the observations, and bump DATA_ASSIMILATION_CYCLES — spin-up effects fade after the first day.
Run an OSSE: swap in obs/synthetic/ from Tutorial 2 and measure recovery of a known truth.
More members: raise ninst and ens_size, and watch spread, inflation, and RMSE respond.
More observations: add GLIDER_* or BOTTLE_* types in Tutorial 1 — start them in evaluate mode.
Explore pyDARTdiags, the DART documentation, and the gallery’s diagnostics section.