Data Assimilation with DART–CESM

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

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

Generate the model domain.
Create a multi-instance (ensemble) CESM case.
Configure the case for assimilation.
Stage your observations
Build and submit the case.
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

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 1: The model domain¶

Step 1.1: Set Up Your Experiment Parameters¶

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

# Change this WORKSPACE path if you installed CROCODILEworkspace somewhere other then /glade/work/$USER
WORKSPACE = Path(os.path.expandvars("/glade/work/$USER/crocodile2026/workspace")) # CROCODILE workspace directory 

START = datetime.datetime(2023, 6, 15)   # must match RUN_STARTDATE in Tutorial 1,2
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 = "/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/obs/real"
SYNTHETIC_OBS_DIR = "/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/obs/synthetic"
# When you have run tutorial 1 and 2 you can use your own generated observations:
#REAL_OBS_DIR      = WORKSPACE / "obs" / "real"        # DART_OBS_ROOT for real obs (Tutorial 1)
#SYNTHETIC_OBS_DIR = WORKSPACE / "obs" / "synthetic"   # DART_OBS_ROOT for synthetic obs (Tutorial 2)

# CESM's DART interface builds its observation file list from DART_OBS_ROOT/{comp}_obs_seq;
# "ocn" is the component name for this ocean-only tutorial series, so the obs_seq files
# Tutorials 1 and 2 write, and Tutorial 3 reads, live in <REAL_OBS_DIR>/ocn_obs_seq and
# <SYNTHETIC_OBS_DIR>/ocn_obs_seq.
OCN_OBS_SEQ_DIR = "ocn_obs_seq"

Step 2.2: Generate the tutorial domain¶

Domain generation is taught in the CrocoDash tutorial. Here run the three cells that build a 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("/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/CrocoDash/data/gebco_2026/GEBCO_2026_coarse_x8.nc")

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 
# Change this WORKSPACE path if you installed ROCODILEworkspace somewhere other then /glade/work/$USER
cesmroot = Path(os.path.expandvars("/glade/work/$USER/crocodile2026/CESM_DA"))

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

# CESM case directory
caseroot = Path(WORKSPACE) / 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="UCGD0009",
    override=True,
    machine="derecho",
    compset="CR_JRA_DA", 
)

Global Alternative: raw create_newcase

If you are not using CrocoDash (for example, running CESM on a global grid), the case can be created directly with CIME. G_JRA_DA is a 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 covers 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.

We’ll build and run our case on Derecho, so from Casper connect to Derecho with:

ssh casper
ssh derecho

Step 4.1: Turn on data assimilation¶

XML settings turn a regional ocean case into a cycling DA experiment. We run xmlchange in the case directory, so in the terminal cd to your case directory:

cd /glade/work/$USER/crocodile2026/workspace/hawaii/

We’re going to set several XML variables to:

Turn on data assimilation for the ocean.
Stop the model each day to assimilate observations - we’re going to run 3 one-day cycles, centered on midnight.
Turn off short term archiving so our output files stay in the run directory - for this tutorial we only have 3 days of output.
Tell CESM-DART where to find our observations.

./xmlchange DATA_ASSIMILATION_OCN=TRUE
./xmlchange STOP_OPTION=ndays,STOP_N=1
./xmlchange DATA_ASSIMILATION_CYCLES=3
./xmlchange DOUT_S=FALSE
./xmlchange DART_OBS_ROOT=/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/obs/real

Notes on xmlquery¶

You can see what CESM variables are set to with ./xmlquery.

./xmlquery CALENDAR

For data assimilation we use the GREGORIAN calendar, rather than the default CESM NO_LEAP calendar.

You can query multiple variables at once:

./xmlquery STOP_OPTION STOP_N DATA_ASSIMILATION_CYCLES DOUT_S JOB_QUEUE NTASKS ROOTPE_OCN DART_OBS_ROOT DATA_ASSIMILATION

Or search for matches to a string:

./xmlquery --partial DATA_ASSIM

Notes on observation files¶

Step 4.2 Set the DART namelist Options¶

The namelist options for dart are set in user_nl_dart. Note user_nl_dart applies to all active DA components. If you are running DA with multiple components you may want to set component specific namelist options. Since we are just running ocean assimilation, we’l use user_nl_dart.

The are several key namelists for ocean DA: &model_nml, &obs_kind_nml, &assim_tools_nml and &filter_nml.

Open user_nl_dart in your case directory, using “open from path”

/work/crocodile2026/workspace/hawaii/user_nl_dart

Copy and past the following namelist options into user_nl_dart and remember to click save!

&model_nml
use_pseudo_depth = .true. 
layer_name = 'Layer'
model_state_variables        = 'Salt ', 'QTY_SALINITY             ', 'NA', 'NA', 'UPDATE',
                               'Temp ', 'QTY_POTENTIAL_TEMPERATURE', 'NA', 'NA', 'UPDATE',
                               'u    ', 'QTY_U_CURRENT_COMPONENT  ', 'NA', 'NA', 'UPDATE',
                               'v    ', 'QTY_V_CURRENT_COMPONENT  ', 'NA', 'NA', 'UPDATE',
/

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

&assim_tools_nml
cutoff = 0.008
/

&filter_nml
perturb_from_single_instance = .true.
perturbation_method = 'model'
stages_to_write = 'output', 'preassim'
num_output_state_members = 3
output_members               = .true.
/

During the CESM build step, the namelist file DART will read is written to to Buildconf/dartconf/input.nml.ocn. We can use ./preview_namelists in the case directory to check input.nml.ocn is created correctly:

./preview_namelists --comp esp

./preview_namelists will also create list of observation files DART use in Buildconf/dart.input_data_list.

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

What the Namelist Options Do¶

&model_nml:

Sets model specific options - which MOM6 variables are in the state vector and what vertical coordinate we want to use.

We’ve set Salt, Temperature, U and V current in our model state - only theses state-vector variables are updated by filter.
We’ve set pseudo depth as the vertical.

&obs_kind_nml:

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

For this tutorial we’ll assimilate ARGO_SALINITY and ARGO_TEMPERATURE. After assimilation the DART QC indicates whether an observation was evaluate only or used in the assimilation: DART outgoing quality control.

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.

&filter_nml and &assim_tools_nml are where the main assimilation options are set.

For much more detail on DART options take a look at the filter namelist documentation. The two main options for ensemble DA and inflation and localization.

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.2 rad ≈ 11.5° ≈ 1,275 km at the equator

We’re not using inflation in this tutorial, but it is a great topic for the project time. The default cutoff is 0.2 radians which would be larger than our entire 5°×5° Hawaii domain, so we’ve set cutoff much smaller at 0.008 radians.

Generating Spread¶

If three ensemble members have the same initial conditions and the same forcing data, they will all output the same results. This is no good for data assimilation, which requires some ensemble spread.

Generating ensemble spread is an interesting scientific topic. Typically small random perturbations are added to each ensemble member, then the models are integrated until system can stabilize, reach physical balance, and provide a realistic spread of model states for data assimilation. Alternative approaches in include running a model for a long time and creating an ensemble by using the state from the same day each year. For example, running for 40 years, and taking the 40 January 1st as though they occurred on the same year.

For this tutorial we are going to generate ensemble spread by perturbing the temperature of each ensemble member. Note a realistic data assimilation experiment would take time to spin up the model.

Step 4.3 Using the cpudev Queue on Derecho¶

Since we have a small domain, we’ll reduce the number of tasks per instance.

./xmlchange JOB_QUEUE=develop
./xmlchange NTASKS=16
./xmlchange ROOTPE_OCN=0

CESM will tell us that for our changes to take effect we have to run:

./case.setup --reset

Step 4.4 Atmospheric forcing¶

The JRA atmosphere can take a long time to initialize so we are going to restrict its search to a single year by running this script in the case directory:

/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/set_atm_forcings.sh

Step 4.5: Build the case¶

qcmd -A UCGD0009 -- ./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 ─────────────┘

SECTION 5: Submit the Data Assimilation Experiment¶

Submitting 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 include the case name, dart, component, 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

The files in your run directory will be called something like hawaii.dart.ocn_obs_seq_final.2023-06-18-00000 which is {case}.dart.{component}_obs_seq_final.{timestamp}.

Since we are also outputting state members there will be files for the individual ensemble members named *member_0001, *member_0002, *member_0003.

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.

We’ll look at used vs. rejected observations, and the prior and posterior statistics for an observation type.

import pydartdiags.obs_sequence.obs_sequence as obsq
from pydartdiags.matplots import matplots as mp

# Your case run directory: cd caseroot && ./xmlquery --value RUNDIR
RUN_DIR = Path(os.path.expandvars("/glade/derecho/scratch/$USER/hawaii/run"))

# Which assimilation cycle to examine. There are DATA_ASSIMILATION_CYCLES=3 of these
# (2023-06-15/16/17/18-00000) - pick any one, not necessarily the last.
TIMESTAMP = "2023-06-18-00000"

def plot_obs_space(timestamp, obs_type, levels):
    # File name convention: {case}.dart.{component}_obs_seq_final.{timestamp}
    data_file = RUN_DIR / f"{casename}.dart.ocn_obs_seq_final.{timestamp}"
    obs_seq = obsq.ObsSequence(data_file)
    obs_seq.possible_vs_used()
    mp.plot_profile(obs_seq, levels, obs_type, bias=True, rmse=True, totalspread=True, depth=True)
    return obs_seq

Let’s take a look at a profile plot, and what this is telling us about the assimilation. We’ll take a look at the Argo temperature observations, binned at the depths below.

Let’s also print a summary of the the number of observations that were available (possible) and the number of observations that were used by filter, with obs_seq.possible_vs_used()

If you have observations that were not used, take a deeper dive into the DART QC value.

levels = [0.0, 100.0, 150.0, 200.0, 250.0, 300.0, 400.0, 500.0, 700, 850, 925, 1000]  # depth in m
obs_seq = plot_obs_space(TIMESTAMP, "ARGO_TEMPERATURE", levels)
obs_seq.possible_vs_used()

Interpreting the assimilation diagnostics¶

When you look at your plot, consider the following questions:

Did the assimilation change the bias? Compare the prior and posterior bias curves. Where do they differ, and at what depths? A change in bias indicates that the assimilation has moved the ensemble mean relative to the observations. Ideally, assimilation should reduce the magnitude of the bias, although this may not occur everywhere.
Did the assimilation change the RMSE? Compare the prior and posterior RMSE curves. Where does the RMSE increase or decrease? A decrease indicates better agreement between the ensemble and the observations at those depths. Also look at the overall RMSE reported in the statistics box.
Where does the assimilation have the largest effect? Look for depths where the prior and posterior curves are most separated. These are regions where the observations are having the greatest impact on the analysis.
Does the impact occur where observations are available? Compare the locations of the observations with changes in the prior and posterior curves. Do the largest changes occur near observed depths? If not, the assimilation may be spreading observational information to other depths through the ensemble covariance and localization.
What happened to the ensemble spread? Compare the prior and posterior total spread. Does the spread decrease after assimilation? A decrease is expected because the observations provide additional information and constrain the ensemble.
Did the spread decrease everywhere, or only at particular depths? Look at the shape of the prior and posterior spread curves. A localized reduction in spread can indicate where the observations are constraining the ensemble most strongly.
Did the reduction in spread correspond to an improvement in RMSE or bias? This is an important question. A smaller posterior spread means the ensemble is more tightly clustered, but it does not necessarily mean that the analysis is more accurate. Compare the changes in spread with the changes in bias and RMSE.
Has the ensemble become too tightly constrained? If the posterior spread is much smaller than the prior spread but the RMSE has not improved, this may indicate that the ensemble has become underdispersive. In other words, the ensemble is more confident without necessarily being more accurate.
What do the overall statistics tell you? Compare the prior and posterior values in the Grand statistics box. Do the changes in the overall bias, RMSE, and total spread agree with what you see in the depth profiles?

Step 6.2: State-space diagnostics¶

Look at the model state before and after assimilation.

filter writes the ensemble mean state for each stage in stages_to_write. Comparing the mean before assimilation (preassim_mean.nc) to the mean after assimilation (output_mean.nc) gives the increment:

increment = posterior mean - prior mean

Increments show where the observations moved the model state. plot_state_space(TIMESTAMP) plots the surface Temp field for the cycle given by TIMESTAMP (the same cycle you looked at in Step 6.1), pre- and post-assimilation, the increment between them, and the locations of the observations that were actually assimilated during that cycle.

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

def lon180(lon):
    """Convert 0-360 longitudes (DART convention) to -180..180 for plotting."""
    lon = np.asarray(lon, dtype=float)
    return np.where(lon > 180, lon - 360, lon)

# MOM6 is a C-grid: T-point fields (Temp, Salt) live on (lath, lonh), but u lives on
# (lath, lonq) and v on (latq, lonh). ocean_geometry.nc has a matching geolon/geolat
# pair for each point type - pick the right one from the field's own dimensions.
_GEOM_SUFFIX = {
    ("lath", "lonh"): "",   # T-point
    ("lath", "lonq"): "u",  # zonal velocity (Cu) point
    ("latq", "lonh"): "v",  # meridional velocity (Cv) point
    ("latq", "lonq"): "b",  # corner (Bu) point
}

def geo_coords(da, geom):
    dims = da.dims
    if dims not in _GEOM_SUFFIX:
        raise ValueError(f"No matching geolon/geolat in ocean_geometry.nc for dims {dims}")
    suffix = _GEOM_SUFFIX[dims]
    geolon = (dims, lon180(geom[f"geolon{suffix}"].values))
    geolat = (dims, geom[f"geolat{suffix}"].values)
    return geolon, geolat

def plot_state_space(timestamp, var="Temp", level=0, obs_type="ARGO_TEMPERATURE"):
    # File name convention: {case}.dart.{component}_{preassim,output}_mean.{timestamp}.nc
    preassim_file = RUN_DIR / f"{casename}.dart.ocn_preassim_mean.{timestamp}.nc"
    output_file   = RUN_DIR / f"{casename}.dart.ocn_output_mean.{timestamp}.nc"
    obs_file      = RUN_DIR / f"{casename}.dart.ocn_obs_seq_final.{timestamp}"

    preassim  = xr.open_dataset(preassim_file)[var].isel(Layer=level).squeeze()
    posterior = xr.open_dataset(output_file)[var].isel(Layer=level).squeeze()

    # Land points are filled with exactly 0 (not NaN). The gray land feature added
    # below hides them visually, but an unmasked 0 still pulls vmin down to 0 when
    # autoscaling the preassim/posterior color limits, squeezing the real ocean range
    # into a sliver of the colormap. Mask them out before plotting.
    preassim  = preassim.where(preassim != 0)
    posterior = posterior.where(posterior != 0)

    increment = posterior - preassim

    # preassim_mean.nc / output_mean.nc carry no lon/lat coordinates, only grid-index
    # dimensions. ocean_geometry.nc is MOM6's static grid file and has the actual
    # geolon/geolat for each point type (0-360, converted to -180..180 to match LON_MIN/MAX).
    geom_file = (sorted(RUN_DIR.glob("*ocean_geometry*")) or sorted(Path(".").glob("*ocean_geometry*")))[-1]
    geom = xr.open_dataset(geom_file)
    geolon, geolat = geo_coords(preassim, geom)

    preassim  = preassim.assign_coords(geolon=geolon, geolat=geolat)
    posterior = posterior.assign_coords(geolon=geolon, geolat=geolat)
    increment = increment.assign_coords(geolon=geolon, geolat=geolat)

    # Observations from the same cycle, so the overlay matches the state fields above.
    obs_seq = obsq.ObsSequence(obs_file)
    used_obs = obs_seq.select_used_qcs()
    used_obs = used_obs[used_obs["type"] == obs_type]

    fig, axes = plt.subplots(1, 3, figsize=(16, 5), subplot_kw={"projection": ccrs.PlateCarree()})
    fields = [preassim, posterior, increment]
    titles = ["Preassim (prior) mean", "Output (posterior) mean", "Increment (posterior - prior)"]
    vmax = float(np.abs(increment).max())

    for ax, da, title in zip(axes, fields, titles):
        is_increment = da is increment
        # geolon/geolat are the coordinates matching this field's point type, degrees.
        da.plot(
            ax=ax, x="geolon", y="geolat", transform=ccrs.PlateCarree(),
            cmap="RdBu_r" if is_increment else "viridis",
            vmin=-vmax if is_increment else None,
            vmax=vmax if is_increment else None,
            cbar_kwargs={"shrink": 0.7, "label": var},
        )
        ax.coastlines(resolution="10m")
        ax.add_feature(cfeature.LAND, facecolor="lightgray")
        ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
        ax.set_title(title)

    # Overlay the observation locations that were actually assimilated.
    axes[2].scatter(
        lon180(used_obs["longitude"]), used_obs["latitude"],
        s=20, edgecolor="k", facecolor="none", transform=ccrs.PlateCarree(),
        label=f"assimilated {obs_type}",
    )
    axes[2].legend(loc="lower left")

    plt.tight_layout()
    plt.show()

    print(f"Preassim file: {preassim_file.name}")
    print(f"Output file:   {output_file.name}")
    print(f"Obs file:      {obs_file.name}")
    print(f"Max |increment|: {vmax:.4f}")

#plot_state_space(TIMESTAMP, var="Temp", level=30, obs_type="ARGO_TEMPERATURE")
plot_state_space(TIMESTAMP)

Try it: connect increment to observation

Find the largest temperature increment on the map. Overlay the assimilated observation locations from used_obs (watch the longitude convention. The model grid may be 0–360). Is there an observation at the bull’s-eye?
Repeat the plot for salinity: plot_state_space(TIMESTAMP, var="Salt", obs_type="ARGO_SALINITY"). Did temperature observations move the salt field? Why is that possible?
Localization: if you halved cutoff in user_nl_dart and reran, how would the footprint of each increment change?
Compare cycles: call plot_obs_space and plot_state_space with a different TIMESTAMP (e.g. "2023-06-18-00000"). Since both functions take the cycle they should examine as an argument, you’re always looking at obs-space and state-space diagnostics from the same cycle.

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.
If you reduce the cutoff, increment footprints shrink, 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.
Spin up an ensemble for a case you are interested in. How will you create spread?
Dive into the DART documentation, review the CESM DART_interface documentation, explore pyDARTdiags for observation space diagnostics, and try the Crocodile gallery’s mom6-tools page to look the mean ensemble member.