Creating Synthetic Observations with perfect_model_obs

Sample a model state as if it were the real ocean - the first step in an Observing System Simulation Experiment (OSSE)!

A typical workflow for generating synthetic observations consists of the following steps:

Design an observing network: we’ll use the locations from Tutorial 1, plus some random ones.
Run a “truth” or nature run of your model: MOM6 for this tutorial
Use dartobsgen with perfect_model_obs create the synthetic observations
Examine the synthetic observations you have created.

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

Missed Tutorial 1?

You can run this notebook standalone:

Point REAL_OBS_DIR at the staged workshop copy, /glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/obs/real`, in Step 3.1.
or
Skip the harvest entirely. Step 3.1 falls back to a regular-grid network if it finds no Tutorial 1 output.

SECTION 1: Why synthetic observations?¶

In an OSSE you treat a model run as the “truth”. You sample your model run where your instruments would sample the real ocean, and add realistic errors such as instrument noise and representativeness error. Because you know the truth exactly, you can use OSSEs to estimate how new instruments and observing networks will impact your data assimilation.

OSSEs answer questions like:

Is my DA plumbing working at all? (If you can’t recover a known truth, real obs won’t help.)
Where should new instruments go to constrain the circulation feature I care about?
How does observation accuracy trade off against observation count?

perfect_model_obs is DART’s tool for the sampling step: it reads a model state, applies the same forward operators the assimilation would use, and perturbs each value with noise drawn from the observation error variance you assign.

SECTION 2: Set Up Your Experiment Parameters¶

dartobsgen’s PerfectModelSource drives perfect_model_obs for you, one assimilation window at a time. It needs three things:

the compiled perfect_model_obs executable,
an input.nml containing a &perfect_model_obs_nml block
a MOM6 state to sample, the “truth”.

For this tutorial, we’ll use a MOM6 run of the Hawaii domain as model truth, and the perfect_model_obs from /glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/pmo.

# --- 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/crocodile2026
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 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      = 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)

OCN_OBS_SEQ_DIR = "ocn_obs_seq"

# Shared, read-only template provided for this tutorial: the compiled
# perfect_model_obs, its input.nml, and the static/geometry/template files
# input.nml names. PerfectModelSource needs to create and remove scratch
# files inside its run directory, so you can't point PMO_RUN_DIR at this
# shared copy directly -- the next cell makes you a private, writable copy.
PMO_RUN_DIR_SOURCE = Path("/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/pmo")
PMO_RUN_DIR = WORKSPACE / "pmo_runs" / "pmo"   # your private writable copy

# MOM6 "truth" run output, one file per assimilation window (or a glob /
# multi-timeslice file). Must contain the variables named by
# model_state_variables in PMO_RUN_DIR/input.nml.
MOM6_HISTORY_GLOB = "/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/hawaii_pmo/ocn/hist/hawaii_pmo.mom6.h.daily.*.nc"

import shutil

# Make your own writable copy of the shared pmo_run_dir template. Skipped if
# you already have one, so re-running this cell won't overwrite local edits
# (e.g. if you've customized input.nml). Delete PMO_RUN_DIR and rerun this
# cell to pick up a fresh copy from PMO_RUN_DIR_SOURCE.
if PMO_RUN_DIR.exists():
    print(f"{PMO_RUN_DIR} already exists; skipping copy.")
else:
    shutil.copytree(
        PMO_RUN_DIR_SOURCE, PMO_RUN_DIR,
        ignore=shutil.ignore_patterns("windows"),  # scratch dir, not part of the template
    )
    print(f"Copied {PMO_RUN_DIR_SOURCE} -> {PMO_RUN_DIR}")

For each window, PerfectModelSource writes a template obs_seq.in holding your network (locations, types, times, error variances), sets the DART namelist for each window, runs the executable in an isolated subdirectory under <PMO_RUN_DIR>/windows/, and collects the resulting obs_seq.out.

SECTION 3: Design the observation network¶

Step 3.1: Harvest the Tutorial 1 locations¶

Often in an OSSE we want to sample the truth where instruments actually were. Let’s use the locations of the real observations from Tutorial 1. We read the first Tutorial 1 window with pyDARTdiags and turn its unique locations into ObsNetworkEntry objects, one per location per observation type.

One thing changes compared to Tutorial 1: the error variance is now yours to choose. Real converters carry the instrument error with the data; in an OSSE the “instrument” is imaginary, so you decide how good it is. We use (0.2 °C)² for temperature and (0.1 PSU)² for salinity, typical Argo-like values.

import numpy as np
import pydartdiags.obs_sequence.obs_sequence as obsq
from dartobsgen import ObsNetworkEntry

# DART's MOM6 model_mod forward operator returns salinity in kg/kg, not PSU,
# so the error variance must be in kg/kg^2 even though we're thinking in PSU.
OBS_ERR_VAR = {
    "ARGO_TEMPERATURE": 0.04,    # (0.2 degC)^2
    "ARGO_SALINITY":    1.0e-8,  # (0.0001 kg/kg = 0.1 psu)^2
}
# Argo-like profile depths (m), used as a fallback if no Tutorial 1 obs
PROFILE_DEPTHS = [10.0, 50.0, 100.0, 200.0, 500.0, 1000.0]  

nb1_files = sorted((REAL_OBS_DIR / OCN_OBS_SEQ_DIR).glob("obs_seq.*.out"))
# Missed Tutorial 1? Use the staged workshop copy instead:
# nb1_files = sorted((Path("/glade/campaign/cgd/oce/projects/CROCODILE/workshops/2026/dart/obs/real") / OCN_OBS_SEQ_DIR).glob("obs_seq.*.out"))

network = []
if nb1_files:
    real = obsq.ObsSequence(str(nb1_files[0]))
    locations = (real.df[["longitude", "latitude", "vertical"]]
                 .drop_duplicates()
                 .reset_index(drop=True))
    print(f"Harvested {len(locations)} unique locations from {nb1_files[0].name}")
    for _, row in locations.iterrows():
        lon = row.longitude if row.longitude <= 180 else row.longitude - 360
        for obs_type in OBS_TYPES:
            network.append(ObsNetworkEntry(
                obs_type=obs_type,
                lat=float(row.latitude),
                lon=float(lon),
                vertical=float(row.vertical),
                vert_unit="height (m)",
                obs_err_var=OBS_ERR_VAR[obs_type],
            ))
else:
    # Fallback: no Tutorial 1 output found -- build a coarse regular-grid network.
    print("No Tutorial 1 output found; building a 1-degree grid network instead.")
    for lat in np.arange(LAT_MIN + 0.5, LAT_MAX, 1.0):
        for lon in np.arange(LON_MIN + 0.5, LON_MAX, 1.0):
            for depth in PROFILE_DEPTHS:
                for obs_type in OBS_TYPES:
                    network.append(ObsNetworkEntry(
                        obs_type=obs_type,
                        lat=float(lat), lon=float(lon),
                        vertical=depth, vert_unit="height (m)",
                        obs_err_var=OBS_ERR_VAR[obs_type],
                    ))

n_harvested = len(network)
print(f"Network so far: {n_harvested} observations")

Step 3.2: Add random locations¶

Now the part you can’t do with real data: invent instruments. We add 20 random profile locations inside the domain, each sampling the standard depths. The seeded random generator means everyone in the workshop gets the same “random” network. You can change the seed to generate a different network.

rng = np.random.default_rng(seed=42)   # fixed seed: reproducible "random" network
N_RANDOM = 20

random_lats = rng.uniform(LAT_MIN, LAT_MAX, N_RANDOM)
random_lons = rng.uniform(LON_MIN, LON_MAX, N_RANDOM)

for lat, lon in zip(random_lats, random_lons):
    for depth in PROFILE_DEPTHS:
        for obs_type in OBS_TYPES:
            network.append(ObsNetworkEntry(
                obs_type=obs_type,
                lat=float(lat), lon=float(lon),
                vertical=depth, vert_unit="height (m)",
                obs_err_var=OBS_ERR_VAR[obs_type],
            ))

print(f"Added {len(network) - n_harvested} obs at {N_RANDOM} random profile locations "
      f"({len(network)} total)")

Step 3.3: Map the network¶

Plot the two networks, the locations harvested from Tutorial 1 and the random additions, together.

import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

net_df = pd.DataFrame([{
    "obs_type": e.obs_type, "lat": e.lat, "lon": e.lon,
    "vertical": e.vertical, "obs_err_var": e.obs_err_var,
} for e in network])

harvested = net_df.iloc[:n_harvested]
random_part = net_df.iloc[n_harvested:]

fig, ax = plt.subplots(figsize=(6, 5), subplot_kw={"projection": ccrs.PlateCarree()})
ax.scatter(harvested["lon"], harvested["lat"], s=18, color="tab:blue",
           label=f"Tutorial 1 locations ({len(harvested)})", transform=ccrs.PlateCarree())
ax.scatter(random_part["lon"], random_part["lat"], s=30, color="tab:orange",
           marker="^", label=f"random additions ({len(random_part)})", transform=ccrs.PlateCarree())

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.gridlines(draw_labels=True)

ax.set_title("Synthetic observing network")
ax.legend()
plt.show()

SECTION 4: Generate The Synthetic Observations¶

Step 4.1: Configure `dartobsgen` For Our New Observation Network¶

ObsGenConfig is almost the same as Tutorial 1. We want the same time period, same bounding box, same observation types, same assimilation frequency, but we want our output written to a different directory: SYNTHETIC_OBS_DIR. The source for “observations” is now a MOM6 run rather than CrocoLake, so we use PerfectModelSource rather than CrocoLake. Because the “truth” changes from day to day, we also give it a MOM6StateProvider built over the MOM6 history output, so each window samples the model state valid at that day rather.

from dartobsgen import (
    MOM6StateProvider,
    ObsGenConfig,
    PerfectModelSource,
    generate_obs_sequences,
    state_vars_from_nml,
)

config = ObsGenConfig(
    start=START, end=END,
    lat_min=LAT_MIN, lat_max=LAT_MAX,
    lon_min=LON_MIN, lon_max=LON_MAX,
    obs_types=OBS_TYPES,
    assimilation_frequency=FREQ,
    output_dir=SYNTHETIC_OBS_DIR / OCN_OBS_SEQ_DIR,
)

# Validates the daily history files against exactly the variables
# PMO_RUN_DIR/input.nml expects; generate_obs_sequences below automatically
# runs check_coverage(), which reports whether the files' times line up
# with the assimilation windows before anything runs.
provider = MOM6StateProvider(
    MOM6_HISTORY_GLOB,
    cache_dir=str(WORKSPACE / "pmo_runs" / "state_cache"),
    required_vars=state_vars_from_nml(str(PMO_RUN_DIR / "input.nml")),
)

source = PerfectModelSource(
    pmo_run_dir=str(PMO_RUN_DIR),
    obs_network=network,
    state_provider=provider,
)

(SYNTHETIC_OBS_DIR / OCN_OBS_SEQ_DIR).mkdir(parents=True, exist_ok=True)

# max_workers=1 runs windows sequentially; set to None to run them in parallel.
written = generate_obs_sequences(config, source, max_workers=1)

print(f"{len(written)} obs_seq file(s) written to {SYNTHETIC_OBS_DIR.name}/{OCN_OBS_SEQ_DIR}/")
for p in written:
    print("  ", Path(p).name)

Step 4.2: Look inside¶

perfect_model_obs has sampled the MOM6 truth state at each location with the DART forward operator and added noise drawn from your obs_err_var. Before moving on, let’s verify that the following is true in your observation sequence files:

ARGO_TEMPERATURE should read as plausible °C for the Hawaii domain (roughly 20-28 °C near the surface), varying smoothly with depth.
ARGO_SALINITY is in kg/kg, not PSU — DART’s MOM6 forward operator returns salinity in kg/kg regardless of what units you’re used to thinking in. A realistic value is ~0.034-0.035, not ~35. If you see values near zero, or negative, salinity isn’t real — it means noise is swamping the signal, almost always because obs_err_var was set in PSU² instead of kg/kg².
obs_err_var for temperature and salinity will look very different in magnitude (e.g. 0.04 vs 1e-8) — that’s expected, since they’re different quantities in different units. What matters is that sqrt(obs_err_var) is small relative to the actual observed value, not that the two numbers look similar to each other.

print("Let's check min/max/mean by type (catches unit/sign errors head() might miss):\n")

for w in written:
    syn = obsq.ObsSequence(str(Path(w)))
    print(f"{Path(w).name}: {len(syn.df)} observations")
    print(syn.df.groupby("type")[["observation", "obs_err_var"]].agg(["min", "max", "mean"]))
    print()

syn.df[["type", "longitude", "latitude", "vertical", "observation", "obs_err_var"]].head(8)

Recap¶

Your takeaway artifact: A directory obs/synthetic/ocn_obs_seq containing synthetic observations under your WORKSPACE.

Where to go from here?¶

Tutorial 3: Cycling DART–CESM: assimilate the real observations, then swap in obs/synthetic/ocn_obs_seq to run the experiment as an OSSE.
Take a look at the script version of this workflow: mom6_perfect_model.py in the dartobsgen repository.
Dive deeper into the DART perfect_model_obs documentation.