Data Assimilation with DART–CESM

Run a real forecast–assimilate–update cycle: a 3-member regional MOM6 ensemble in CESM, updated every 24 hours (one day) by the observations you made in Tutorial 1.

A typical DART-CESM experiment consists of the following steps:

Install CESM_DA.
(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. Working with Real Observations · 2. Creating Synthetic Observations · 3. Data Assimilation with DART–CESM

SECTION 1: Install the CESM DA fork¶

DART is integrated as the ESP component of CESM in CROCODILE. It is built and called by CESM like any other component following the setup, build, submit workflow. For this tutorial, we’ll setup the case in CrocoDash, set options using ./xmlchange and user_nl_dart, build the case with ./case.build, and submit the case with ./case.submit.

Clone the fork, check out the DA branch, and populate the components:

git clone https://github.com/CROCODILE-CESM/CESM.git CESM_DA
cd CESM_DA/
git checkout full_regional_cesm_da
./bin/git-fleximod update

You will point cesmroot at this clone in Section 3.

For more detail on the DART interface to CESM, such as what compsets are available with DA, see the DART_interface documentation. For much more detail on DART take a look at the DART documentation or learn about research with DART at dart.ucar.edu.

SECTION 2: The model domain¶

Step 2.1: Set Up Your Experiment Parameters¶

Remember to match these with your observations from Tutorial 1.

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

DA_PROJECT_DIR = Path("<DART_DA_PROJECT_DIR>")        # your DA project working directory on Derecho

START = datetime.datetime(2023, 6, 15)   # must match RUN_STARTDATE in Tutorial 3
END   = datetime.datetime(2023, 6, 18)   # 3 days -> 3 one-day assimilation windows, centered on midnight
FREQ  = datetime.timedelta(hours=24)

# Bounding box:
LAT_MIN, LAT_MAX = 20.0, 25.0
LON_MIN, LON_MAX = -160.0, -155.0

OBS_TYPES = ["ARGO_TEMPERATURE", "ARGO_SALINITY"]

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

Step 2.2: Regenerate the tutorial domain¶

Domain generation is taught in the CrocoDash tutorial. Here we just run the three cells that build the same Hawaii domain, because creating a case needs the live grid, topography, and vertical-grid objects.

from CrocoDash.grid import Grid

grid = Grid(
  resolution = 0.05, # in degrees
  xstart = 200.0, # min longitude in [0, 360]
  lenx = 5.0, # longitude extent in degrees
  ystart = 20.0, # min latitude in [-90, 90]
  leny = 5.0, # latitude extent in degrees
  name = "hawaii",
)

from CrocoDash.topo import Topo

topo = Topo(
    grid=grid,
    min_depth=9.5,  # in meters
)

from pathlib import Path

bathymetry_path = Path("<GEBCO>")

if not bathymetry_path.exists():
    raise FileNotFoundError(
        "Bathymetry file not found, please replace with path to bathymetry file"
    )

topo.set_from_dataset(
    bathymetry_path=bathymetry_path,
    longitude_coordinate_name="lon",
    latitude_coordinate_name="lat",
    vertical_coordinate_name="elevation",
)

topo.depth.plot()

from CrocoDash.vgrid import VGrid

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

SECTION 3: Create a multi-instance CESM case¶

Step 3.1: Set the case details¶

Make sure cesmroot points at the CESM_DA directory.

from pathlib import Path
# CESM case (experiment) name
casename = "hawaii"

# CESM source root - this is your CESM_DA directory from Section 1.
cesmroot = "<CESM_DA>"

# Place where all your input files will be stored
inputdir = Path(DA_PROJECT_DIR) / "input_files" / casename

# CESM case directory
caseroot = Path(DA_PROJECT_DIR) / casename

print(f"Case:")
print(f" casename = {casename}")
print(f" cesmroot = {cesmroot}")
print(f" inputdir = {inputdir}")
print(f" caseroot = {caseroot}")

print(f"Observation dirctories - these will be used once you have created your case")
print(f" REAL_OBS_DIR = {REAL_OBS_DIR}")
print(f" SYNTHETIC_OBS_DIR = {SYNTHETIC_OBS_DIR}")

Step 3.2: Create the case¶

One new concept compared to the CrocoDash tutorial is running multi-instance CESM. An ensemble filter needs an ensemble, that is a group of model forcasts. CESM runs ninst copies (“instances”) of the model. In an assimilation experiment, each ensemble member (instance) starts from slightly different initial conditions. The spread of the ensemble members is what gives us information on model uncertainty.

For more detail on ensemble data assimilation, see the DART introduction to ensemble data assimilation.

Three members is workshop-sized, small enough to build and run in a tutorial session. Real ocean DA experiments use 30–80 members, and spin up the oceans from different initial conditions to get ensemble spread. We will perturb the ensemble members in this tutorial to generate ensemble spread.

from CrocoDash.case import Case

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

Global Alternative: raw create_newcase

If you are not using CrocoDash (for example, on a global grid), the same case can be created directly with CIME. G_JRA_DA is the global DA-enabled ocean compset:

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

Answer

Statistics. 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 available. The ensemble is the model uncertainty.

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: make sure the forcing cover the date range you are interested in.

case.configure_forcings(
    date_range=["2023-06-15 00:00:00", "2023-06-18 00:00:00"],
    boundaries=["north", "south", "east", "west"],
    function_name="get_glorys_data_from_rda",
)

case.process_forcings()

SECTION 4: Configure the case for assimilation¶

We’ll now swap to the terminal to use CESM’s xmlchange command to query and set our experiment options.

Step 4.1: Turn on data assimilation¶

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_DA_PROJECT_DIR>/hawaii/
./case.setup

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

# 3 one-day cycles, centered on midnight:
./xmlchange STOP_OPTION=ndays,STOP_N=1
./xmlchange DATA_ASSIMILATION_CYCLES=3

./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

Check the forcing calendar¶

MOM6 expects the calendar for the input to be named “gregorian” rather than standard.

module load nco
cd inputdir/casename/ocnice
for f in *.nc; do
    ncatted -O -a calendar,time,o,c,"gregorian" "$f" "$f"
done

Using the cpu dev queue¶

./xmlchange JOB_QUEUE=develop
./xmlchange NTASKS=2
+ setting root pes

Step 4.2 Tell CESM when your observations are¶

Unlike regular CESM components, DART also requires a list of observation files, to use in the assimilation. Tell CESM when your observations are, <REAL_OBS_DIR> from the experiment parameters cell in section 2.1

./xmlchange DART_OBS_ROOT="<REAL_OBS_DIR>"

During the CESM build step, the list of observation files DART expects is written to to Buildconf/dart.input_data_list. We can use ./preview_namelists to check the file is created correctly:

./preview_namelists --comp esp

Check Buildconf/dart.input_data_list contains the files you expect.

cat Buildconf/dart.input_data_list

Step 4.3: Build the case¶

qcmd -- ./case.build

The build takes a while. Let’s take a look at how the assimilation will while we wait. The job will execute 3 cycles data assimilation since DATA_ASSIMILATION_CYCLES=3:


               one cycle
            ┌──────────────────────────────────────────────────────────┐
            ▼                                                          │
 1. FORECAST      all 3 MOM6 instances advance 24 hours (1 day)        │
 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 24-hour segment ─────────────┘

Step 4.4: A tour of DART namelists¶

During case.build CESM writes a DART namelist to Buildconf/dartconf/input.nml, which is read-only and regenerated for each run of DART.

To change DART settings, edit user_nl_dart in the case directory, exactly like user_nl_mom for MOM6.

The three key namelists 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 are the default). Only state-vector variables are updated by the filter.

&obs_kind_nml: which observation types are used by filter:

assimilate_these_obs_types - these observations are allowed to impact the model state.
evaluate_these_obs_types - these observations do not impact the model state, but their forward operator is computed and recorded.

After assimilation the DART QC indicates whether an observation was evaluate only or used in the assimilation: DART outgoing quality control.

For this tutorial we’ll assimilate ARGO_SALINITY and ARGO_TEMPERATURE. Add the following to user_nl_dart

&obs_kind_nml
assimilate_these_obs_types = 'ARGO_SALINITY', 'ARGO_TEMPERATURE'
/

The ./preview_namelists --component esp command will create input.nml.ocn with all the input options options for DART. Check the file contains the options you expect:

cat Buildconf/dartconf/input.nml.ocn

For much more detail on DART options take a look at the filter namelist documentation and the MOM6 model_mod documentation.

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 a way to test model performance against independent observations without altering or updating the model’s internal state.

SECTION 5: Submit the Data Assimilation Experiment¶

Submmitting the case will run the DATA_ASSIMILATION_CYCLES in one jobs submission.

./case.submit

Each cycle, filter writes 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 optionally posterior) forward operator results, and DART QC value
`output_mean.nc`	Updated ensemble mean
`output_sd.nc`	Updated ensemble spread

SECTION 6: Examine the assimilation¶

State space diagnostics evaluate model variables and error covariances directly in the system’s physical or model coordinate domain, whereas observation space diagnostics map model estimates to the measurement (observation) domain.

Step 6.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 6.2: State-space diagnostics¶

Look at the model state before and after assimilation.

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 a cutoff allows spurious correlations to update the model state.

Recap¶

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

Where to go from here?¶

Run longer: set END = 2023-09-30 in Tutorial 1, regenerate the observations, and increase DATA_ASSIMILATION_CYCLES.
Run an OSSE: swap in synthetic observations from Tutorial 2.
More members: raise ninst and examine how spread, inflation, and RMSE respond.
More observations: add GLIDER_* or BOTTLE_* types in Tutorial 1. Start them in evaluate mode.
Dive into the DART documentation, explore pyDARTdiags for observation space diagnostics, and try the Crocodile gallery’s mom6-tools page to look the mean ensemble member.