diff --git a/.gitignore b/.gitignore
index 50b2ff7d..833a323f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -14,7 +14,7 @@ _dev
 .vscode
 .idea
 .snakemake
-output
+output*
 _sandbox/
 experiment_report.html
 
diff --git a/config/varda-rapid-0.2.yaml b/config/varda-rapid-0.2.yaml
new file mode 100644
index 00000000..7461d434
--- /dev/null
+++ b/config/varda-rapid-0.2.yaml
@@ -0,0 +1,158 @@
+# yaml-language-server: $schema=../workflow/tools/config.schema.json
+description: |
+  Evaluate skill of Varda-rapid against ground observations.
+
+config_label: varda-rapid-0.2
+
+dates:
+  # start: 2025-01-01T06:00
+  # end: 2025-12-31T06:00
+  # frequency: 294h
+  - 2025-03-01T00:00
+  - 2025-04-02T06:00
+  - 2025-07-02T12:00
+  - 2025-11-02T18:00
+
+runs:
+  - temporal_downscaler:
+      checkpoint: https://service.meteoswiss.ch/mlstore#/models/sruc-m-2-interpolator/versions/3
+      label: Varda-rapid-0.2
+      steps: 0/12/1
+      config: resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi_with_nudging.yaml
+      extra_requirements:
+        - anemoi-datasets==0.5.35
+        - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@0e69b5916166a7feb104f892cdf2a091bd45bd56
+        - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
+      forecaster:
+        checkpoint: https://service.meteoswiss.ch/mlstore#/models/sruc-m-1-forecaster/versions/4
+        config: resources/inference/configs/sgm-multidataset-forecaster-global-ich1-oper.yaml
+        steps: 0/12/6
+        extra_requirements:
+          - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@0e69b5916166a7feb104f892cdf2a091bd45bd56
+          - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
+
+  - temporal_downscaler:
+      checkpoint: https://service.meteoswiss.ch/mlstore#/models/sruc-m-2-interpolator/versions/3
+      label: Varda-single-1.0
+      steps: 0/12/1
+      config: resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi.yaml
+      extra_requirements:
+        - anemoi-datasets==0.5.35
+        - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@cc116a80c86346cb8de4449c1091bb2c6299a302
+        - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
+      forecaster:
+        checkpoint: https://service.meteoswiss.ch/mlstore#/models/sruc-m-1-forecaster/versions/4
+        config: resources/inference/configs/sgm-multidataset-forecaster-global-ich1-oper.yaml
+        steps: 0/12/6
+        extra_requirements:
+          - git+https://github.com/MeteoSwiss/anemoi-plugins-meteoswiss@cc116a80c86346cb8de4449c1091bb2c6299a302
+          - git+https://github.com/ecmwf/anemoi-inference@88c1ec13632c570f246cc2580ce324e45f36a03b
+
+  - baseline:
+      label: INCA
+      root: /store_new/mch/msclim/INCA
+      steps: 0/6/1
+  - baseline:
+      label: ICON-CH2-CTRL
+      root: /store_new/mch/msopr/osm/ICON-CH2-EPS
+      steps: 0/12/1
+  - baseline:
+      label: ICON-CH1-CTRL
+      root: /store_new/mch/msopr/osm/ICON-CH1-EPS
+      steps: 0/12/1
+
+truth:
+  label: SwissMetNet
+  root: jretrievedwh:bbox=45.7,48.0,5.8,10.8;use_limitation=40
+
+experiment:
+  params:
+    - T_2M
+    - TD_2M
+    - U_10M
+    - V_10M
+    - TOT_PREC
+    - PMSL
+    - PS
+  stratification:
+    regions:
+      - icon
+      - switzerland
+      - mittelland
+      - berge
+      - alpennordseite
+      - alpensuedseite
+      - jura
+    root: /store_new/mch/msopr/ml/regions/Prognoseregionen_LV95_20220517
+  thresholds:
+    TOT_PREC:
+      gt: [0.0, 1, 5]
+    U_10M:
+      gt: [2.5, 5.0, 10.0]
+    V_10M:
+      gt: [2.5, 5.0, 10.0]
+    T_2M:
+      lt: [273.15]
+      gt: [288.15, 298.15]
+  dashboard:
+    stratification:
+      - season
+      # - region
+  station_holdout:
+    - CHM
+    - FRE
+    - RUE
+    - BAS
+    - PLF
+    - CGI
+    - HOE
+    - PUY
+    - RAG
+    - EIN
+    - ALT
+    - GOS
+    - COV
+    - ULR
+    - BIN
+    - GUE
+    - ROE
+    - BIA
+    - CEV
+    - OTL
+  scoremaps:
+    enabled: false
+
+showcase:
+  params:
+    - T_2M
+    - TD_2M
+    - U_10M
+    - V_10M
+  meteograms:
+    enabled: false
+    stations: [JUN, KLO, LUG, GVE] #, COV, GOR, WFJ, SAE, SAM, DAV, ZER, ANT, VSBAS, BRT, LTB, GOS, CEV, BIA]
+  animations:
+    enabled: true
+    frames_per_second: 0.5
+    domains:
+      # - globe
+      # - europe
+      - alps
+      - icon-ch
+      - switzerland
+
+locations:
+  output_root: output/
+
+profile:
+  executor: slurm
+  global_resources:
+    gpus: 16
+  default_resources:
+    slurm_partition: "normal"
+    cpus_per_task: 1
+    mem_mb_per_cpu: 1800
+    runtime: "1h"
+    gpus: 0
+    slurm_account: s83
+  jobs: 50
diff --git a/resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi_with_nudging.yaml b/resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi_with_nudging.yaml
new file mode 100644
index 00000000..1936dd15
--- /dev/null
+++ b/resources/inference/configs/sgm-temporal-downscaler-global_trimedge_multi_with_nudging.yaml
@@ -0,0 +1,176 @@
+runner: temporal_downscaler
+
+input:
+  cutout:
+    - lam_0:
+        grib:
+          path: forecaster/20*
+          pre_processors:
+            - forward_transform_filter:
+                # does not have an effect if the temporal downscaler does not use tp as prognostic variable
+                rescale:
+                  scale: 0.001 # convert from kg m-2 to m
+                  offset: 0
+                  param: tp
+            # nudge the initial condition toward station observations.
+            - forward_transform_filter:
+                retrieve_observation:
+                  obs_path: observation/nudging_station_obs_raw.parquet
+                  jretrieve_src_path: /scratch/mch/llanzila/sruc/evalml/src
+                  retrieval_bbox: [40.5, 53.0, 0.0, 17.5]
+                  # station_group: "1,2"
+                  variables: [T_2M, TD_2M, U_10M, V_10M, TOT_PREC, PMSL, PS]
+                  use_limitation: 40
+                  run_mode: "devt"
+                  station_filter_mode: domain
+                  trim_bbox: [45.7, 48.0, 5.8, 10.8]
+            - forward_transform_filter:
+                clean_observation:
+                  obs_path_in: observation/nudging_station_obs_raw.parquet
+                  obs_path_out: observation/nudging_station_obs.parquet
+            - forward_transform_filter:
+                nudge_toward_observation:
+                  obs_path: observation/nudging_station_obs.parquet
+                  nudge_variables: [T_2M, TD_2M, U_10M, V_10M, PMSL, PS]
+                  icon_grid_file: /scratch/mch/llanzila/sruc/aux_files/icon_grid_0001_R19B08_mch.nc
+                  icon_orog_file: /scratch/mch/icontest/testing-input-data/c2sm/icon-1/external_parameter_icon_grid_0001_R19B08_mch.nc
+                  d_eff_file: /scratch/mch/llanzila/sruc/aux_files/d_eff_cache_domain_maxdist50km_nbar50x3_bw1500m_elev50_elevdiff100_nsta1939.nc
+                  weight_power: 2.0
+                  max_dist: 50000.0
+                  use_reliability_check: true
+                  number_of_std: 5.0
+                  reliability_min_dist_frac: 0.03
+                  write_diagnostics: true
+                  diagnostics_dir: nudging_diagnostics
+                  run_mode: "devt"
+          namer: &namer
+            rules:
+              - - shortName: T
+                - t_{level}
+              - - shortName: U
+                - u_{level}
+              - - shortName: V
+                - v_{level}
+              - - shortName: W
+                - w_{level}
+              - - shortName: QV
+                - q_{level}
+              - - shortName: FI
+                - z_{level}
+              - - shortName: PMSL
+                - msl
+              - - shortName: FIS
+                - z
+              - - shortName: PS
+                - sp
+              - - shortName: T_2M
+                - 2t
+              - - shortName: TD_2M
+                - 2d
+              - - shortName: T_G
+                - skt
+              - - shortName: U_10M
+                - 10u
+              - - shortName: V_10M
+                - 10v
+              - - shortName: FR_LAND
+                - lsm
+              - - shortName: TOT_PREC
+                - tp
+    - global:
+        grib:
+          path: forecaster/ifs*
+          namer: *namer
+
+constant_forcings:
+  test:
+    use_original_paths: true
+
+patch_metadata: resources/sgm-temporal-downscaler-ich1-oper-patch.yaml
+
+post_processors:
+  - accumulate_from_start_of_forecast: # accumulate tp from start of forecast
+      accumulations:
+        - tp
+  - forward_transform_filter:
+      rescale:
+        scale: 1000 # convert units from m to kg m-2
+        offset: 0
+        param: tp
+
+output:
+  tee:
+    - grib:
+        path: grib/{dateTime}_{step:03}.grib
+        encoding:
+          typeOfGeneratingProcess: 2
+        templates:
+            samples: resources/templates_index_icon.yaml
+        post_processors:
+          - extract_mask: # removes global points
+              mask: "lam_0/cutout_mask"
+              as_slice: true
+          # here, the trimedge mask can be specified when available
+          # at overlap steps (multiples of the forecaster stride) replace the
+          # re-predicted prognostics with the forecaster's boundary values
+          - forward_transform_filter:
+              copy-prognostic-from-forecaster:
+                forecaster_path: forecaster/20*
+                common_leadtime: 6h
+                params_to_keep: [tp] # NOTE, if more diagnostics are added, update accordingly.
+                namer: *namer
+
+    - grib:
+        path: grib/ifs-{dateTime}_{step:03}.grib
+        encoding:
+          typeOfGeneratingProcess: 2
+        templates:
+          samples: resources/templates_index_ifs.yaml
+        post_processors:
+          - extract_mask: # removes lam points
+              mask: "lam_0/cutout_mask"
+              as_slice: true
+              inverse: true
+          - assign_mask: # fill local/global overlapping points with nan
+              mask: "global/cutout_mask"
+          - forward_transform_filter:
+              copy-prognostic-from-forecaster:
+                forecaster_path: forecaster/ifs*
+                common_leadtime: 6h
+                params_to_keep: [tp] # NOTE, if more diagnostics are added, update accordingly.
+                namer: *namer
+        modifiers:
+          - patches:
+            - variable:
+                U_10M: {"param": 165, "shortName": "10u"}
+                10u: {"param": 165, "shortName": "10u"}
+                V_10M: {"param": 166, "shortName": "10v"}
+                10v: {"param": 166, "shortName": "10v"}
+                TD_2M: {"param": 168, "shortName": "2d"}
+                2d: {"param": 168, "shortName": "2d"}
+                T_2M: {"param": 167, "shortName": "2t"}
+                2t: {"param": 167, "shortName": "2t"}
+                FR_LAND: {"param": 172, "shortName": "lsm"}
+                lsm: {"param": 172, "shortName": "lsm"}
+                PMSL: {"param": 151, "shortName": "msl"}
+                msl: {"param": 151, "shortName": "msl"}
+                PS: {"param": 134, "shortName": "sp"}
+                sp: {"param": 134, "shortName": "sp"}
+                SSO_SIGMA: {"param": 163, "shortName": "slor"}
+                slor: {"param": 163, "shortName": "slor"}
+                SSO_STDH: {"param": 160, "shortName": "sdor"}
+                sdor: {"param": 160, "shortName": "sdor"}
+                TOT_PREC: {"param": 228, "shortName": "tp"}
+                tp: {"param": 228, "shortName": "tp"}
+                z: {"param": 129, "shortName": "z"}
+                "^q_(\\d+)$": {"param": 133, "shortName": "q"}
+                "^t_(\\d+)$": {"param": 130, "shortName": "t"}
+                "^u_(\\d+)$": {"param": 131, "shortName": "u"}
+                "^v_(\\d+)$": {"param": 132, "shortName": "v"}
+                "^w_(\\d+)$": {"param": 135, "shortName": "w"}
+                "^z_(\\d+)$": {"param": 129, "shortName": "z"}
+
+# silenced due to bug in anemoi-inference for multi-step temporal downscalers, can be removed when fixed
+verbosity: 0
+allow_nans: true
+output_frequency: "1h"
diff --git a/resources/report/dashboard/script.js b/resources/report/dashboard/script.js
index 830324a0..d89af495 100644
--- a/resources/report/dashboard/script.js
+++ b/resources/report/dashboard/script.js
@@ -34,6 +34,7 @@ function initChoices(id) {
 initChoices("region-select");
 initChoices("season-select");
 initChoices("init-select");
+initChoices("station-group-select");
 initChoices("source-select");
 initChoices("metric-select");
 initChoices("param-select");
@@ -52,8 +53,9 @@ function getSelected(id) {
   window.DATA = raw.data.map(row => {
     const obj = {};
     for (let i = 0; i < cols.length; i++) obj[cols[i]] = row[i];
+    const sgPart = obj.station_group !== "all" ? ", Group: " + obj.station_group : "";
     obj.region_season_init =
-      "Region: " + obj.region + ", Season: " + obj.season + ", Init: " + obj.init_hour;
+      "Region: " + obj.region + ", Season: " + obj.season + ", Init: " + obj.init_hour + sgPart;
     return obj;
   });
 })();
@@ -154,7 +156,7 @@ async function renderLegend(filteredData) {
             legend: {
               orient: "bottom",
               direction: "horizontal",
-              title: "Region / Season / Init",
+              title: "Region / Season / Init / Group",
               labelLimit: 400,
               symbolType: "circle", symbolSize: 120,
             },
@@ -216,13 +218,16 @@ async function updateChart() {
   const selMetrics = getSelected("metric-select");
   const selParams  = getSelected("param-select");
 
-  // Filter data by region / season / init / source
+  const selStationGroups = getSelected("station-group-select");
+
+  // Filter data by region / season / init / source / station_group
   // (metric and param are handled per cell)
   let filtered = DATA;
-  if (selRegions.length) filtered = filtered.filter(d => selRegions.includes(d.region));
-  if (selSeasons.length) filtered = filtered.filter(d => selSeasons.includes(d.season));
-  if (selInits.length)   filtered = filtered.filter(d => selInits.includes(d.init_hour));
-  if (selSources.length) filtered = filtered.filter(d => selSources.includes(d.source));
+  if (selRegions.length)       filtered = filtered.filter(d => selRegions.includes(d.region));
+  if (selSeasons.length)       filtered = filtered.filter(d => selSeasons.includes(d.season));
+  if (selInits.length)         filtered = filtered.filter(d => selInits.includes(d.init_hour));
+  if (selSources.length)       filtered = filtered.filter(d => selSources.includes(d.source));
+  if (selStationGroups.length) filtered = filtered.filter(d => selStationGroups.includes(d.station_group));
 
   // Show / hide table columns (params)
   document.querySelectorAll("#chart-table thead th[data-param]").forEach(th => {
@@ -327,6 +332,7 @@ function resizeChartScroll() {
       getSelected("source-select"),
       getSelected("metric-select"),
       getSelected("param-select"),
+      getSelected("station-group-select"),
     ].flatMap(v => v);
     summary.textContent = parts.join(", ");
   }
@@ -342,7 +348,7 @@ function resizeChartScroll() {
 
   // Keep summary current when selections change (guard: some selects may be absent)
   ["region-select", "season-select", "init-select",
-   "source-select", "metric-select", "param-select"].forEach(id => {
+   "source-select", "metric-select", "param-select", "station-group-select"].forEach(id => {
     const el = document.getElementById(id);
     if (el) el.addEventListener("change", () => {
       if (panel.classList.contains("collapsed")) updateSummary();
diff --git a/resources/report/dashboard/template.html.jinja2 b/resources/report/dashboard/template.html.jinja2
index a6fe1f4a..0da748f0 100644
--- a/resources/report/dashboard/template.html.jinja2
+++ b/resources/report/dashboard/template.html.jinja2
@@ -154,6 +154,16 @@
             </select>
         </div>
         {% endif %}
+        {% if station_holdout %}
+        <div class="control-group">
+            <label>Station group</label>
+            <select id="station-group-select" multiple>
+                {% for sg in station_groups %}
+                <option value="{{ sg | e }}"{% if sg == 'all' %} selected{% endif %}>{{ sg | e }}</option>
+                {% endfor %}
+            </select>
+        </div>
+        {% endif %}
         <div class="control-group">
             <label>Source(s)</label>
             <select id="source-select" multiple>
diff --git a/src/data_input/__init__.py b/src/data_input/__init__.py
index 2af07f7d..c9915628 100644
--- a/src/data_input/__init__.py
+++ b/src/data_input/__init__.py
@@ -686,6 +686,56 @@ def _jretrieve_df_to_xarray(df, short_names, catalog) -> xr.Dataset:
     return xr.Dataset(data_vars=data_vars, coords=coords)
 
 
+def _trim_stations_xr(
+    ds: xr.Dataset, filter_mode: str, domain_bbox: list | None
+) -> xr.Dataset:
+    """Trim the "values" (station) dimension to *filter_mode* — same logic as
+    RetrieveObservation._trim_stations / notebooks/d_eff_generator.ipynb's
+    station-trim cell, adapted for an xarray Dataset with latitude/longitude
+    coords along "values" instead of a pandas DataFrame."""
+    lat = ds["latitude"].values
+    lon = ds["longitude"].values
+
+    if filter_mode == "domain":
+        lat_min, lat_max, lon_min, lon_max = domain_bbox
+        mask = (lat >= lat_min) & (lat <= lat_max) & (lon >= lon_min) & (lon <= lon_max)
+        desc = f"domain bbox {domain_bbox}"
+
+    elif filter_mode == "switzerland":
+        import cartopy.io.shapereader as shpreader
+        from shapely.geometry import Point
+
+        shp_path = shpreader.natural_earth(
+            resolution="10m", category="cultural", name="admin_0_countries"
+        )
+        ch_country = next(
+            r
+            for r in shpreader.Reader(shp_path).records()
+            if r.attributes["ADM0_A3"] == "CHE"
+        )
+        swiss_geom = ch_country.geometry
+
+        def _in_switzerland(lat_, lon_):
+            if pd.isna(lat_) or pd.isna(lon_):
+                return False
+            return swiss_geom.contains(Point(lon_, lat_))
+
+        mask = np.array([_in_switzerland(la, lo) for la, lo in zip(lat, lon)])
+        desc = "Swiss national border (Natural Earth)"
+
+    else:
+        raise ValueError(
+            f"Unknown filter_mode: {filter_mode!r} (expected 'domain' or 'switzerland')"
+        )
+
+    n_before = ds.sizes["values"]
+    ds = ds.isel(values=mask)
+    LOG.info(
+        "Station filter [%s]: %d -> %d stations", desc, n_before, ds.sizes["values"]
+    )
+    return ds
+
+
 def load_obs_data_from_jretrieve(
     root, reftime: datetime, steps: list[int], params: list[str]
 ) -> xr.Dataset:
@@ -721,7 +771,9 @@ def load_obs_data_from_jretrieve(
 
     from data_input import jretrieve as jr
 
-    stations, stage, seq_type = jr.parse_selection(root)
+    stations, stage, seq_type, use_limitation, filter_mode, domain_bbox = (
+        jr.parse_selection(root)
+    )
     jr.check_prerequisites(stage)
 
     want_uv = "U_10M" in params or "V_10M" in params
@@ -751,9 +803,13 @@ def load_obs_data_from_jretrieve(
         increment_minutes=step_hours * 60,
         seq_type=seq_type,
         stage=stage,
+        use_limitation=use_limitation,
     )
     raw = _jretrieve_df_to_xarray(df, short_names, catalog)
 
+    if filter_mode is not None:
+        raw = _trim_stations_xr(raw, filter_mode, domain_bbox)
+
     out = xr.Dataset(coords=raw.coords)
     for icon, short in DWH_PARAM_MAP.items():
         if icon in params and short in raw:
@@ -773,7 +829,18 @@ def load_obs_data_from_jretrieve(
 
     out = out.dropna("values", how="all")
     times = np.datetime64(reftime) + np.asarray(steps, dtype="timedelta64[h]")
-    return _select_valid_times(out, times, strict=True)
+    result = _select_valid_times(out, times, strict=True)
+
+    # Per-variable station-coverage log, so ground-truth coverage can be compared
+    # against what was actually available to nudge.
+    _logged_params = ("T_2M", "TD_2M", "U_10M", "V_10M", "PMSL", "TOT_PREC", "VMAX_10M")
+    n_total = result.sizes["values"]
+    for icon in _logged_params:
+        if icon in result.data_vars:
+            n_valid = int(result[icon].notnull().any("time").sum())
+            LOG.info("Stations with valid %s: %d / %d stations", icon, n_valid, n_total)
+
+    return result
 
 
 def load_truth_data(
diff --git a/src/data_input/jretrieve.py b/src/data_input/jretrieve.py
index ada7de32..8e72f7c2 100644
--- a/src/data_input/jretrieve.py
+++ b/src/data_input/jretrieve.py
@@ -168,20 +168,34 @@ def _stations_to_argv(stations: dict[str, Any]) -> list[str]:
     raise AssertionError("unreachable")
 
 
-def parse_selection(root: Any) -> tuple[dict[str, Any], str, str]:
-    """Parse a truth-root marker into (stations, stage, seq_type).
+def parse_selection(
+    root: Any,
+) -> tuple[dict[str, Any], str, str, int | None, str | None, list | None]:
+    """Parse a truth-root marker into
+    (stations, stage, seq_type, use_limitation, filter_mode, domain_bbox).
 
     Examples (slash-free so they survive ``Path()`` normalisation):
       ``jretrievedwh:SwissMetNet``                       -> group
       ``jretrievedwh:group=SwissMetNet;stage=devt``
       ``jretrievedwh:locations=ARO,KLO``
       ``jretrievedwh:bbox=45.8,47.8,5.9,10.5``
+      ``jretrievedwh:bbox=40.5,53.0,0.0,17.5;use_limitation=40;filter_mode=switzerland``
+        -> retrieve over the given bbox, then trim to stations within the real
+           Swiss national border (see load_obs_data_from_jretrieve).
+      ``jretrievedwh:bbox=40.5,53.0,0.0,17.5;filter_mode=domain;domain_bbox=45.7,48.0,5.8,10.8``
+        -> retrieve over one bbox, then trim to a second, independent bbox.
+
+    Defaults when not given: use_limitation=40, filter_mode=None and
+    domain_bbox=None (no extra trim).
     """
     _, _, rest = str(root).partition(":")
     rest = rest.strip()
     stations: dict[str, Any] = {}
     stage = "prod"
     seq_type = "surface"
+    use_limitation: int | None = 40
+    filter_mode: str | None = None
+    domain_bbox: list | None = None
     for i, part in enumerate([p for p in rest.split(";") if p]):
         if "=" not in part:
             if i == 0:
@@ -196,11 +210,25 @@ def parse_selection(root: Any) -> tuple[dict[str, Any], str, str]:
             stage = value
         elif key == "seq_type":
             seq_type = value
+        elif key == "use_limitation":
+            use_limitation = int(value)
+        elif key == "filter_mode":
+            if value not in ("domain", "switzerland"):
+                raise ValueError(
+                    f"filter_mode must be 'domain' or 'switzerland', got {value!r}"
+                )
+            filter_mode = value
+        elif key == "domain_bbox":
+            domain_bbox = [float(v) for v in value.split(",") if v]
+            if len(domain_bbox) != 4:
+                raise ValueError("domain_bbox must be lat_min,lat_max,lon_min,lon_max.")
         else:
             raise ValueError(f"Unknown jretrieve selector key: {key!r}")
     if not stations:
         stations = {"group": DEFAULT_GROUP}
-    return stations, stage, seq_type
+    if filter_mode == "domain" and domain_bbox is None:
+        raise ValueError("domain_bbox is required when filter_mode='domain'.")
+    return stations, stage, seq_type, use_limitation, filter_mode, domain_bbox
 
 
 def _run(argv: list[str], env: dict[str, str], timeout_s: int) -> str:
@@ -299,7 +327,7 @@ def fetch_data(
     increment_minutes=60,
     seq_type="surface",
     stage="prod",
-    use_limitation: int = 40,
+    use_limitation: int | None = 40,
     timeout_s=600,
 ) -> pd.DataFrame:
     """Fetch observation data; columns: station (int), termin (YYYYMMDDhhmmss),
@@ -318,7 +346,8 @@ def fetch_data(
         "csv",
         *_stations_to_argv(stations),
     ]
-    argv += ["--use-limitation", str(use_limitation)]
+    if use_limitation is not None:
+        argv += ["--use-limitation", str(use_limitation)]
     LOG.info("jretrieve data: %s", " ".join(argv))
     return _parse_csv(_run_with_retry(argv, env=_build_env(stage), timeout_s=timeout_s))
 
diff --git a/src/evalml/cli.py b/src/evalml/cli.py
index f5c8a2cd..ab3c530b 100644
--- a/src/evalml/cli.py
+++ b/src/evalml/cli.py
@@ -26,6 +26,7 @@ def _base_snakemake_command(
     command += config.profile.parsable()
     command += ["--configfile", str(configfile)]
     command += ["--cores", str(cores)]
+    command += ["--rerun-incomplete"]
     return command
 
 
diff --git a/src/evalml/config.py b/src/evalml/config.py
index 3980046a..364dda61 100644
--- a/src/evalml/config.py
+++ b/src/evalml/config.py
@@ -455,6 +455,40 @@ class Dashboard(BaseModel):
     )
 
 
+class StationHoldoutConfig(BaseModel):
+    """Station holdout settings for station-group stratified verification."""
+
+    holdout_fraction: Optional[float] = Field(
+        default=None,
+        gt=0,
+        lt=1,
+        description=(
+            "Fraction of truth stations to hold out for evaluation (exclusive 0–1). "
+            "Mutually exclusive with exclude_stations."
+        ),
+    )
+    holdout_seed: int = Field(
+        default=42,
+        description="Random seed for reproducible holdout station selection.",
+    )
+    exclude_stations: Optional[List[str]] = Field(
+        default=None,
+        description=(
+            "Explicit list of station nat_abbr to hold out. "
+            "Mutually exclusive with holdout_fraction."
+        ),
+    )
+
+    @model_validator(mode="after")
+    def _holdout_fraction_and_exclude_stations_are_exclusive(self):
+        if self.holdout_fraction is not None and self.exclude_stations is not None:
+            raise ValueError(
+                "station_holdout: holdout_fraction and exclude_stations are "
+                "mutually exclusive; set at most one of them."
+            )
+        return self
+
+
 class ExperimentConfig(BaseModel):
     """Configuration for the experiment workflow outputs."""
 
@@ -477,6 +511,15 @@ class ExperimentConfig(BaseModel):
         ...,
         description="Settings for the experiment dashboard.",
     )
+    station_holdout: Optional[Union[List[str], StationHoldoutConfig]] = Field(
+        default=None,
+        description=(
+            "When set, adds a 'station_group' dimension (all/holdout/holdin) to verification "
+            "metrics for all models and baselines. Either a bare list of station nat_abbr "
+            "(shorthand for exclude_stations), or a full object specifying holdout_fraction "
+            "(+ optional holdout_seed) instead."
+        ),
+    )
     scorecards: Optional[ExperimentScorecardConfig] = Field(
         default=None,
         description="Scorecard generation configuration. Omit or set enabled: false to disable.",
@@ -486,6 +529,14 @@ class ExperimentConfig(BaseModel):
         description="Score map plot configuration. Omit or set enabled: false to disable.",
     )
 
+    @field_validator("station_holdout", mode="before")
+    @classmethod
+    def normalize_station_holdout(cls, v):
+        """Accept a bare list of station nat_abbr as shorthand for exclude_stations."""
+        if isinstance(v, list):
+            return {"exclude_stations": v}
+        return v
+
     @field_validator("thresholds")
     @classmethod
     def validate_threshold_operators(
diff --git a/src/verification/__init__.py b/src/verification/__init__.py
index ae818d3a..0b5126d0 100644
--- a/src/verification/__init__.py
+++ b/src/verification/__init__.py
@@ -300,6 +300,34 @@ def _merge_metrics(ds: xr.Dataset, num_workers: int = 4) -> xr.Dataset:
     return out
 
 
+def _create_station_group_masks(
+    values_coord: xr.DataArray, holdout_stations: list[str] | None
+) -> xr.DataArray:
+    """Boolean masks over the values dimension for station groups.
+
+    Always includes "all" (every verify() output carries a station_group dimension,
+    so downstream consumers never need to check whether it's present). When
+    holdout_stations is empty/None, "all" is the only group — the holdout/holdin
+    split is only meaningful, and only added, when stations are actually held out.
+    """
+    all_nat_abbr = values_coord.values
+    if not holdout_stations:
+        return xr.DataArray(
+            np.ones((1, len(all_nat_abbr)), dtype=bool),
+            coords={"station_group": ["all"], "values": all_nat_abbr},
+            dims=["station_group", "values"],
+        )
+    is_holdout = np.isin(all_nat_abbr, holdout_stations)
+    return xr.DataArray(
+        np.stack([np.ones_like(is_holdout), is_holdout, ~is_holdout]),
+        coords={
+            "station_group": ["all", "holdout", "holdin"],
+            "values": all_nat_abbr,
+        },
+        dims=["station_group", "values"],
+    )
+
+
 def _bbox_polygon(lon_min, lon_max, lat_min, lat_max) -> Polygon:
     return Polygon(
         [
@@ -321,6 +349,7 @@ def verify(
     dim: list[str] | None = None,
     threshold_dict: dict[str, dict[str, list[float]]] | None = None,
     num_workers: int | None = None,
+    holdout_stations: list[str] | None = None,
     max_missing_fraction: float = 0.0,
 ) -> xr.Dataset:
     """
@@ -354,6 +383,11 @@ def verify(
         If None, no thresholds used.
     num_workers : int or None, optional
         Number of parallel workers for computation. If None, uses available CPU cores minus 2.
+    holdout_stations : list[str] or None, optional
+        Station nat_abbr forming the "holdout" station group. The output always
+        carries a station_group dimension — ["all"] when this is None/empty, or
+        ["all", "holdout", "holdin"] when stations are given — so callers never need to check
+        for its presence.
     max_missing_fraction : float, optional
         Maximum allowed fraction of missing forecast values among obs-valid in-region points
         before a metric is set to NaN. Computed per region and time step over the reduction
@@ -396,6 +430,14 @@ def verify(
         lon=obs_aligned["longitude"], lat=obs_aligned["latitude"]
     )
 
+    station_masks = _create_station_group_masks(obs_aligned["values"], holdout_stations)
+    if holdout_stations:
+        LOG.info(
+            "Station group masks created: %d holdout, %d holdin stations",
+            int(station_masks.sel(station_group="holdout").sum()),
+            int(station_masks.sel(station_group="holdin").sum()),
+        )
+
     scores = []
     for param in fcst_aligned.data_vars:
         if param not in obs_aligned.data_vars:
@@ -417,6 +459,10 @@ def verify(
         fcst_param = fcst_aligned[param].where(masks)
         obs_param = obs_aligned[param].where(masks)
 
+        # Apply station group masks — adds "station_group" dim, always present
+        # (["all"], or ["all", "holdout", "holdin"] when holdout_stations is set)
+        fcst_param = fcst_param.where(station_masks)
+        obs_param = obs_param.where(station_masks)
         # Missing fraction: among obs-valid in-region points, fraction where fcst is missing.
         # Normalising by obs availability avoids penalising parameters with fewer stations.
         missing_fraction = (
diff --git a/tests/unit/test_jretrieve.py b/tests/unit/test_jretrieve.py
index 7539dd0c..04094963 100644
--- a/tests/unit/test_jretrieve.py
+++ b/tests/unit/test_jretrieve.py
@@ -37,11 +37,17 @@ def test_parse_selection_default_group():
         {"group": "SwissMetNet"},
         "prod",
         "surface",
+        40,
+        None,
+        None,
     )
     assert jr.parse_selection("jretrievedwh:SwissMetNet") == (
         {"group": "SwissMetNet"},
         "prod",
         "surface",
+        40,
+        None,
+        None,
     )
 
 
@@ -50,6 +56,40 @@ def test_parse_selection_keyvalue_and_stage():
         {"locations": "ARO,KLO"},
         "devt",
         "surface",
+        40,
+        None,
+        None,
+    )
+
+
+def test_parse_selection_use_limitation_and_filter_mode_switzerland():
+    assert jr.parse_selection(
+        "jretrievedwh:bbox=40.5,53.0,0.0,17.5;use_limitation=40;filter_mode=switzerland"
+    ) == (
+        {"bbox": "40.5,53.0,0.0,17.5"},
+        "prod",
+        "surface",
+        40,
+        "switzerland",
+        None,
+    )
+
+
+def test_parse_selection_filter_mode_domain_requires_domain_bbox():
+    with pytest.raises(ValueError, match="domain_bbox"):
+        jr.parse_selection("jretrievedwh:bbox=40.5,53.0,0.0,17.5;filter_mode=domain")
+
+
+def test_parse_selection_filter_mode_domain_with_bbox():
+    assert jr.parse_selection(
+        "jretrievedwh:bbox=40.5,53.0,0.0,17.5;filter_mode=domain;domain_bbox=45.7,48.0,5.8,10.8"
+    ) == (
+        {"bbox": "40.5,53.0,0.0,17.5"},
+        "prod",
+        "surface",
+        40,
+        "domain",
+        [45.7, 48.0, 5.8, 10.8],
     )
 
 
diff --git a/uv.lock b/uv.lock
index eea30160..73af7eb6 100644
--- a/uv.lock
+++ b/uv.lock
@@ -10,17 +10,16 @@ resolution-markers = [
 
 [[package]]
 name = "adjusttext"
-version = "1.3.0"
+version = "1.4.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "matplotlib" },
     { name = "numpy" },
     { name = "scipy" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/4c/d4/6585f3b6fdb75648bca294664af4becc8aa2fb3fb08f4e4e9fd27e10d773/adjusttext-1.3.0.tar.gz", hash = "sha256:4ab75cd4453af4828876ac3e964f2c49be642ea834f0c1f7449558d5f12cbca1", size = 15724, upload-time = "2024-10-31T16:45:36.101Z" }
+sdist = { url = "https://files.pythonhosted.org/packages/b5/5c/496e506ad3313664df24a79f801719cadcd62af5999fcb299c9b08ff0d4b/adjusttext-1.4.0.tar.gz", hash = "sha256:1f73860ced8cccce3f85ee6989ca133c2579b67a7453f63dbeb38f39bf123154", size = 15852, upload-time = "2026-06-08T16:48:32.726Z" }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/53/1c/8feedd607cc14c5df9aef74fe3af9a99bf660743b842a9b5b1865326b4aa/adjustText-1.3.0-py3-none-any.whl", hash = "sha256:da23d7b24b6db5ffa039bb136bfa556207365e32f48ac74b07ad26dd485bc691", size = 13154, upload-time = "2024-10-31T16:45:35.227Z" },
-    { url = "https://files.pythonhosted.org/packages/2d/80/7ad35ee5321a86b842f9e8516c8ae4c86f58db7b40e82ce9759f94517a50/adjusttext-1.3.0-py3-none-any.whl", hash = "sha256:bc6c118cd9d7caf6ae37f9355e51d840a2d7f64b4fb2956b8401de27c5af803b", size = 13264, upload-time = "2026-06-08T16:40:05.041Z" },
+    { url = "https://files.pythonhosted.org/packages/b8/2c/897bdd17b05724c894a5b831c6b2e9853adcc2d07a70d6246c0cd5cd3912/adjusttext-1.4.0-py3-none-any.whl", hash = "sha256:6febd6484c0d45c39a22f44b2c1f4a8cd01ef58fada565cab4b629c771df79b5", size = 13262, upload-time = "2026-06-08T16:48:31.765Z" },
 ]
 
 [[package]]
@@ -3230,15 +3229,15 @@ wheels = [
 
 [[package]]
 name = "plotly"
-version = "6.7.0"
+version = "6.8.0"
 source = { registry = "https://pypi.org/simple" }
 dependencies = [
     { name = "narwhals" },
     { name = "packaging" },
 ]
-sdist = { url = "https://files.pythonhosted.org/packages/3a/7f/0f100df1172aadf88a929a9dbb902656b0880ba4b960fe5224867159d8f4/plotly-6.7.0.tar.gz", hash = "sha256:45eea0ff27e2a23ccd62776f77eb43aa1ca03df4192b76036e380bb479b892c6", size = 6911286, upload-time = "2026-04-09T20:36:45.738Z" }
+sdist = { url = "https://files.pythonhosted.org/packages/94/fd/d72c292d78aadb93d1a9bcd76bf3c678271040c7cf10abe5788b33040a39/plotly-6.8.0.tar.gz", hash = "sha256:e088e7ddc68d4f70e3d66659224727a45296d71d2b8284181862d3d8f1f0d88f", size = 6915161, upload-time = "2026-06-03T18:33:40.226Z" }
 wheels = [
-    { url = "https://files.pythonhosted.org/packages/90/ad/cba91b3bcf04073e4d1655a5c1710ef3f457f56f7d1b79dcc3d72f4dd912/plotly-6.7.0-py3-none-any.whl", hash = "sha256:ac8aca1c25c663a59b5b9140a549264a5badde2e057d79b8c772ae2920e32ff0", size = 9898444, upload-time = "2026-04-09T20:36:39.812Z" },
+    { url = "https://files.pythonhosted.org/packages/f9/14/abe5ce876ab5b66ee3c691bf537fcd43d037aea55d447aacf74630a8f31e/plotly-6.8.0-py3-none-any.whl", hash = "sha256:13c5c4a0f70b74cab1913eda0de49b826df5931708eb6f9c3010040614700ec8", size = 9902055, upload-time = "2026-06-03T18:33:34.26Z" },
 ]
 
 [[package]]
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
index acf76746..8edf0087 100644
--- a/workflow/rules/common.smk
+++ b/workflow/rules/common.smk
@@ -398,14 +398,16 @@ def verif_hash(full_config: dict) -> str:
     """Hash of all settings that affect verification outputs.
 
     Combines the truth source with verification-method settings so that
-    changing either (e.g. switching lapse_rate_correction on/off) produces
-    new output paths and unconditionally triggers a rerun.
+    changing any of them (e.g. switching lapse_rate_correction on/off, or
+    changing the station_holdout selection) produces new output
+    paths and unconditionally triggers a rerun.
     """
     truth_cfg = {
         k: v for k, v in full_config["truth"].items() if k not in TRUTH_HASH_EXCLUDE
     }
     experiment_verif_cfg = {
         "lapse_rate_correction": full_config.get("lapse_rate_correction", True),
+        "station_holdout": full_config.get("experiment", {}).get("station_holdout"),
     }
     return generate_json_hash({"truth": truth_cfg, "verif": experiment_verif_cfg})
 
@@ -423,11 +425,15 @@ def truth_file_dep(_):
 if "jretrieve" in str(config["truth"]["root"]):
     from data_input.jretrieve import check_prerequisites, parse_selection
 
-    _, _jretrieve_stage, _ = parse_selection(config["truth"]["root"])
+    _, _jretrieve_stage, _, _, _, _ = parse_selection(config["truth"]["root"])
     check_prerequisites(_jretrieve_stage)
 
 
 TRUTH_HASH = truth_hash(config["truth"])
+_station_holdout_raw = config.get("experiment", {}).get("station_holdout") or {}
+STATION_HOLDOUT_CFG = (
+    _station_holdout_raw if isinstance(_station_holdout_raw, dict) else {}
+)
 REGIONS = parse_regions()
 VERIF_HASH = verif_hash(config)
 _showcase = config.get("showcase", {})
diff --git a/workflow/rules/inference.smk b/workflow/rules/inference.smk
index 13d7ccc0..e0dd94a9 100644
--- a/workflow/rules/inference.smk
+++ b/workflow/rules/inference.smk
@@ -238,6 +238,12 @@ rule inference_prepare_temporal_downscaler:
             if RUN_CONFIGS[wc.run_id].get("forecaster") is None
             else _get_forecaster_run_id(wc.run_id)
         ),
+        # Single source of truth for the holdout station set: experiment.station_holdout
+        # in the top-level experiment config. Injected into any nudge_toward_observation
+        # filter found in the inference config (a no-op if the config has none, or if
+        # station_holdout isn't configured for this experiment) — see
+        # inference_prepare.py::_inject_nudging_station_holdout.
+        station_holdout_cfg=STATION_HOLDOUT_CFG,
     script:
         "../scripts/inference_prepare.py"
 
@@ -298,7 +304,7 @@ rule inference_execute:
             cd {params.workdir}
 
             _run_inference() {{
-                local VENV=$1
+                local VENV=/user-environment
                 source "$VENV/bin/activate"
 
                 if [ "{params.disable_local_definitions}" = "False" ]; then
@@ -312,20 +318,22 @@ rule inference_execute:
                     CMD_ARGS+=(runner.parallel.cluster=slurm)
                 fi
 
-                srun \
-                    --unbuffered \
-                    {params.account_flag} \
-                    --partition={resources.slurm_partition} \
-                    --cpus-per-task={resources.cpus_per_task} \
-                    --mem-per-cpu={resources.mem_mb_per_cpu} \
-                    --time={resources.runtime} \
-                    --gres={resources.gres} \
-                    --ntasks={resources.ntasks} \
-                    anemoi-inference run config.yaml "${{CMD_ARGS[@]}}"
+                anemoi-inference run config.yaml "${{CMD_ARGS[@]}}"
             }}
             export -f _run_inference
 
-            squashfs-mount {params.env_path}:/user-environment -- bash -c '_run_inference /user-environment'
+            srun \
+                --job-name=anemoi-inference \
+                --uenv={params.env_path}:/user-environment \
+                --unbuffered \
+                {params.account_flag} \
+                --partition={resources.slurm_partition} \
+                --cpus-per-task={resources.cpus_per_task} \
+                --mem-per-cpu={resources.mem_mb_per_cpu} \
+                --time={resources.runtime} \
+                --gres={resources.gres} \
+                --ntasks={resources.ntasks} \
+                bash -c '_run_inference'
         ) >{log} 2>&1
         touch {output.okfile}
         """
diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk
index 0ce78881..76c31b76 100644
--- a/workflow/rules/report.smk
+++ b/workflow/rules/report.smk
@@ -39,6 +39,7 @@ rule report_experiment_dashboard:
         sources=",".join(list(EXPERIMENT_PARTICIPANTS.keys())),
         header_text=make_header_text(),
         stratification=" ".join(config["experiment"]["dashboard"]["stratification"]),
+        station_holdout_flag="--station_holdout" if STATION_HOLDOUT_CFG else "",
         label_map=",".join(
             "{}:{}".format(
                 sid,
@@ -60,6 +61,7 @@ rule report_experiment_dashboard:
             --header_text "{params.header_text}" \
             --configfile "{input.configfile}" \
             --stratification {params.stratification} \
+            {params.station_holdout_flag} \
             --labels "{params.label_map}" \
             --output {output} >{log} 2>&1
         """
diff --git a/workflow/rules/verification.smk b/workflow/rules/verification.smk
index 425476ba..b2037614 100644
--- a/workflow/rules/verification.smk
+++ b/workflow/rules/verification.smk
@@ -1,6 +1,7 @@
 # ----------------------------------------------------- #
 # VERIFICATION WORKFLOW                                 #
 # ----------------------------------------------------- #
+import json
 from datetime import datetime
 
 import pandas as pd
@@ -32,6 +33,7 @@ rule verification_metrics_baseline:
         regions=REGIONS,
         experiment_params=",".join(EXPERIMENT_PARAMS),
         threshold_dict=config["experiment"]["thresholds"],
+        station_holdout_cfg=json.dumps(STATION_HOLDOUT_CFG),
         lapse_rate_flag=(
             "--lapse_rate_correction"
             if config.get("lapse_rate_correction", True)
@@ -42,7 +44,7 @@ rule verification_metrics_baseline:
         export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions)
         uv run {input.script} \
             --forecast {input.forecast} \
-            --truth {params.truth} \
+            --truth "{params.truth}" \
             --reftime {wildcards.init_time} \
             --steps "{params.baseline_steps}" \
             --source_id "{wildcards.baseline_id}" \
@@ -50,6 +52,7 @@ rule verification_metrics_baseline:
             --regions '{params.regions}' \
             --params "{params.experiment_params}" \
             --threshold_dict "{params.threshold_dict}" \
+            --station_holdout_cfg '{params.station_holdout_cfg}' \
             --member "{params.member}" \
             {params.lapse_rate_flag} \
             --output {output} >{log} 2>&1
@@ -91,6 +94,7 @@ rule verification_metrics:
         ).resolve(),
         experiment_params=",".join(EXPERIMENT_PARAMS),
         threshold_dict=config["experiment"]["thresholds"],
+        station_holdout_cfg=json.dumps(STATION_HOLDOUT_CFG),
         lapse_rate_flag=(
             "--lapse_rate_correction"
             if config.get("lapse_rate_correction", True)
@@ -101,7 +105,7 @@ rule verification_metrics:
         export ECCODES_DEFINITION_PATH=$(realpath .venv/share/eccodes-cosmo-resources/definitions)
         uv run {input.script} \
             --forecast {params.grib_out_dir} \
-            --truth {params.truth} \
+            --truth "{params.truth}" \
             --reftime {wildcards.init_time} \
             --steps "{params.fcst_steps}" \
             --source_id "{wildcards.run_id}" \
@@ -109,6 +113,7 @@ rule verification_metrics:
             --regions '{params.regions}' \
             --params "{params.experiment_params}" \
             --threshold_dict "{params.threshold_dict}" \
+            --station_holdout_cfg '{params.station_holdout_cfg}' \
             {params.lapse_rate_flag} \
             --output {output} >{log} 2>&1
         """
@@ -129,6 +134,18 @@ rule verification_metrics_aggregation:
             init_time=_restrict_reftimes_to_hours(REFTIMES),
             allow_missing=True,
         ),
+        inference_okfiles=lambda wc: expand(
+            rules.inference_execute.output.okfile,
+            run_id=(
+                [wc.run_id]
+                + (
+                    [RUN_CONFIGS[wc.run_id]["forecaster"]["run_id"]]
+                    if RUN_CONFIGS[wc.run_id].get("forecaster") is not None
+                    else []
+                )
+            ),
+            init_time=_restrict_reftimes_to_hours(REFTIMES),
+        ),
     output:
         OUT_ROOT / f"data/runs/{{run_id}}/verif_aggregated_{VERIF_HASH}.nc",
     log:
@@ -228,7 +245,7 @@ rule verification_scoremaps:
         uv run {input.script} \
             --run_root {params.run_root} \
             --reftimes {params.reftimes} \
-            --truth {input.truth} \
+            --truth "{input.truth}" \
             --step {wildcards.leadtime} \
             --steps "{params.fcst_steps}" \
             --param {wildcards.param} \
@@ -263,7 +280,7 @@ rule verification_scoremaps_baseline:
         uv run {input.script} \
             --baseline_root {input.forecast} \
             --reftimes {params.reftimes} \
-            --truth {input.truth} \
+            --truth "{input.truth}" \
             --step {wildcards.leadtime} \
             --steps "{params.baseline_steps}" \
             --param {wildcards.param} \
diff --git a/workflow/scripts/inference_prepare.py b/workflow/scripts/inference_prepare.py
index 59ba193d..a7038e19 100644
--- a/workflow/scripts/inference_prepare.py
+++ b/workflow/scripts/inference_prepare.py
@@ -8,7 +8,12 @@
 from evalml.helpers import setup_logger
 
 
-def prepare_config(default_config_path: str, output_config_path: str, params: dict):
+def prepare_config(
+    default_config_path: str,
+    output_config_path: str,
+    params: dict,
+    station_holdout_cfg: dict | None = None,
+):
     """Prepare the configuration file for the inference run.
 
     Overrides default configuration parameters with those provided in params
@@ -22,12 +27,22 @@ def prepare_config(default_config_path: str, output_config_path: str, params: di
         Path where the updated configuration file will be written.
     params : dict
         Dictionary of parameters to override in the default configuration.
+    station_holdout_cfg : dict, optional
+        The experiment's ``experiment.station_holdout`` settings. Only used if
+        the config contains a ``forward_transform_filter:
+        nudge_toward_observation`` block, in which case its
+        ``exclude_stations``/``holdout_fraction``/``holdout_seed`` are
+        injected into that block — see ``_inject_nudging_station_holdout``.
+        Configs without a nudging filter are left untouched.
     """
 
     with open(default_config_path, "r") as f:
         config = yaml.safe_load(f)
 
     config = _override_recursive(config, params)
+    nudging_filter = _find_nudging_filter(config)
+    if nudging_filter is not None:
+        _inject_nudging_station_holdout(nudging_filter, station_holdout_cfg or {})
 
     with open(output_config_path, "w") as f:
         yaml.safe_dump(config, f, sort_keys=False)
@@ -93,7 +108,14 @@ def prepare_temporal_downscaler(smk):
 
     # prepare config
     overrides = _overrides_from_params(smk)
-    prepare_config(smk.input.config, smk.output.config, overrides)
+    station_holdout_cfg = getattr(smk.params, "station_holdout_cfg", None)
+    prepare_config(
+        smk.input.config,
+        smk.output.config,
+        overrides,
+        station_holdout_cfg=station_holdout_cfg,
+    )
+
     LOG.info("Wrote config file at %s", smk.output.config)
     with open(smk.output.config, "r") as f:
         config_content = f.read()
@@ -170,6 +192,68 @@ def _override_recursive(original: dict, updates: dict) -> dict:
     return original
 
 
+def _find_nudging_filter(config: dict) -> dict | None:
+    """Return the ``forward_transform_filter: nudge_toward_observation`` block
+    of config, or None if there is none. By design an inference config
+    declares at most one nudging filter; more than one raises ValueError.
+    Unlike _override_recursive (dict-into-dict only), this walks into lists
+    too, since the filter typically sits inside a pre_processors list.
+    """
+    found = []
+
+    def _walk(node):
+        if isinstance(node, dict):
+            transform_filter = node.get("forward_transform_filter")
+            if isinstance(transform_filter, dict) and isinstance(
+                transform_filter.get("nudge_toward_observation"), dict
+            ):
+                found.append(transform_filter["nudge_toward_observation"])
+            for value in node.values():
+                _walk(value)
+        elif isinstance(node, list):
+            for item in node:
+                _walk(item)
+
+    _walk(config)
+    if len(found) > 1:
+        raise ValueError(
+            f"Found {len(found)} nudge_toward_observation filters in the "
+            "inference config; at most one is supported."
+        )
+    return found[0] if found else None
+
+
+def _inject_nudging_station_holdout(
+    nudging_filter: dict, station_holdout_cfg: dict
+) -> None:
+    """Set exclude_stations/holdout_fraction/holdout_seed on the
+    nudge_toward_observation block nudging_filter, in place, from the
+    experiment's station_holdout settings, so the holdout station set is
+    defined once, in the experiment config.
+
+    A no-op if station_holdout_cfg has neither exclude_stations nor
+    holdout_fraction set (e.g. station holdout isn't configured for this
+    experiment) — any hand-written value already in the config is then left
+    untouched.
+    """
+    exclude_stations = station_holdout_cfg.get("exclude_stations")
+    holdout_fraction = station_holdout_cfg.get("holdout_fraction")
+    if exclude_stations is None and holdout_fraction is None:
+        return
+    holdout_seed = station_holdout_cfg.get("holdout_seed", 42)
+
+    # Mutually exclusive in NudgeTowardObservation itself — clear both
+    # before setting the one station_holdout_cfg actually specifies,
+    # so a stale hand-written value of the other can never linger.
+    nudging_filter.pop("exclude_stations", None)
+    nudging_filter.pop("holdout_fraction", None)
+    if exclude_stations is not None:
+        nudging_filter["exclude_stations"] = list(exclude_stations)
+    else:
+        nudging_filter["holdout_fraction"] = holdout_fraction
+        nudging_filter["holdout_seed"] = holdout_seed
+
+
 def main(smk):
     """Main function to run the Snakemake workflow."""
     if smk.rule == "inference_prepare_forecaster":
diff --git a/workflow/scripts/report_experiment_dashboard.py b/workflow/scripts/report_experiment_dashboard.py
index b8d3abf8..c44b9e4e 100644
--- a/workflow/scripts/report_experiment_dashboard.py
+++ b/workflow/scripts/report_experiment_dashboard.py
@@ -90,7 +90,6 @@ def main(args):
     if "init_hour" not in stratification:
         df = df[df["init_hour"] == "all"]
 
-    # create a new column for line styles and shapes in dashboard
     df.dropna(inplace=True)
     LOG.info("Loaded verification data frame: \n%s", df)
 
@@ -102,6 +101,10 @@ def main(args):
     seasons = df["season"].unique() if "season" in stratification else []
     init_hours = df["init_hour"].unique() if "init_hour" in stratification else []
 
+    # station_group selector: auto-detect from data (shown when station_holdout is enabled)
+    station_groups = sorted(df["station_group"].unique())
+    station_holdout = args.station_holdout and station_groups != ["all"]
+
     # Columnar JSON: store columns + data array (no repeated keys per row).
     # region_season_init is a derived column — computed in JS at parse time.
     # Round float values to 6 significant digits to avoid unnecessary precision.
@@ -122,6 +125,7 @@ def _round_sig(x, sig=6):
         "region",
         "season",
         "init_hour",
+        "station_group",
     ]
     df_export = df[export_cols].copy()
     df_export["value"] = df_export["value"].apply(
@@ -164,6 +168,8 @@ def _sanitize(v):
         seasons=seasons,
         init_hours=init_hours,
         stratification=stratification,
+        station_holdout=station_holdout,
+        station_groups=station_groups,
         header_text=args.header_text,
         configfile_content=open(args.configfile, "r").read()
         if args.configfile.is_file()
@@ -211,6 +217,12 @@ def _sanitize(v):
         default=["region", "season", "init_hour"],
         help="Stratification dimensions to include in the dashboard (any of region, season, init_hour).",
     )
+    parser.add_argument(
+        "--station_holdout",
+        action="store_true",
+        default=False,
+        help="When set, shows a station-group selector (All / Holdout / Hold-in) in the dashboard.",
+    )
     parser.add_argument(
         "--configfile",
         type=Path,
diff --git a/workflow/scripts/report_scorecard.py b/workflow/scripts/report_scorecard.py
index 93a4a8d6..924e409d 100644
--- a/workflow/scripts/report_scorecard.py
+++ b/workflow/scripts/report_scorecard.py
@@ -228,6 +228,11 @@ def _load_relative_diff(cfg: dict) -> xr.Dataset:
     model_ds = xr.open_dataset(cfg["model"]["path"])
     baseline_ds = xr.open_dataset(cfg["baseline"]["path"])
 
+    # Scorecards always compare the "all" station group (verify() always carries
+    # this dimension, holdout or not) — select it and drop the dim.
+    model_ds = model_ds.sel(station_group="all", drop=True)
+    baseline_ds = baseline_ds.sel(station_group="all", drop=True)
+
     if strat_dim != "region":
         sel_coords["region"] = model_ds["region"].values[0]
 
diff --git a/workflow/scripts/verification_metrics.py b/workflow/scripts/verification_metrics.py
index 85a79879..277185d5 100644
--- a/workflow/scripts/verification_metrics.py
+++ b/workflow/scripts/verification_metrics.py
@@ -5,6 +5,7 @@
 from datetime import datetime
 from pathlib import Path
 
+import numpy as np
 
 from verification import verify, apply_lapse_rate_correction_inplace  # noqa: E402
 from verification.spatial import map_forecast_to_truth  # noqa: E402
@@ -31,6 +32,30 @@ class ScriptConfig(Namespace):
     steps: list[int] = parse_steps("0/120/6")
 
 
+def compute_holdout_stations(
+    all_stations: list, station_holdout_cfg: dict
+) -> list[str]:
+    """Return nat_abbr list of holdout stations derived from the truth dataset's station list.
+
+    Either the explicit exclude_stations present in the truth dataset, or a random
+    holdout_fraction of the truth stations drawn with holdout_seed.
+
+    holdout_fraction's (0, 1) range is validated once, at config-load time, by
+    StationHoldoutConfig (evalml.config) — not re-checked here.
+    """
+    exclude_stations = station_holdout_cfg.get("exclude_stations")
+    holdout_fraction = station_holdout_cfg.get("holdout_fraction")
+    holdout_seed = station_holdout_cfg.get("holdout_seed", 42)
+
+    if exclude_stations is not None:
+        return [s for s in exclude_stations if s in all_stations]
+    if holdout_fraction is not None:
+        n_holdout = round(len(all_stations) * float(holdout_fraction))
+        rng = np.random.default_rng(holdout_seed)
+        return list(rng.choice(all_stations, size=n_holdout, replace=False))
+    return []
+
+
 def program_summary_log(args):
     """Log a welcome message with the script information."""
     LOG.info("=" * 80)
@@ -88,6 +113,29 @@ def main(args: ScriptConfig):
     if args.lapse_rate_correction:
         apply_lapse_rate_correction_inplace(fcst, truth, args.params)
 
+    # determine holdout stations for station-holdout stratification
+    # holdout stations are derived from the truth dataset's station list so the
+    # same partition is used consistently across all models and baselines.
+    holdout_stations = None
+    station_holdout_cfg = args.station_holdout_cfg
+    if station_holdout_cfg and "values" in truth.dims:
+        all_stations = list(truth["values"].values)
+        holdout_stations = compute_holdout_stations(all_stations, station_holdout_cfg)
+        if holdout_stations:
+            LOG.info(
+                "Station holdout: %d / %d stations in holdout group",
+                len(holdout_stations),
+                len(all_stations),
+            )
+        else:
+            LOG.warning(
+                "station_holdout_cfg set but no holdout stations selected (check holdout_fraction / exclude_stations)."
+            )
+    elif station_holdout_cfg:
+        LOG.warning(
+            "station_holdout_cfg set but truth dataset has no 'values' dimension; station stratification skipped."
+        )
+
     # compute metrics and statistics
     now = datetime.now()
     results = verify(
@@ -97,6 +145,7 @@ def main(args: ScriptConfig):
         args.truth_source_id,
         regions=args.regions,
         threshold_dict=args.threshold_dict,
+        holdout_stations=holdout_stations,
     )
     LOG.info(
         "Computed verification metrics in %s seconds",
@@ -187,6 +236,16 @@ def main(args: ScriptConfig):
         help="Dictionary of thresholds for each parameter in the format '{param: [threshold1, threshold2, ...]}' (default: None).",
         default=None,
     )
+    parser.add_argument(
+        "--station_holdout_cfg",
+        type=lambda s: json.loads(s) if s else None,
+        default=None,
+        help=(
+            "Station holdout config as a JSON dict with keys: holdout_fraction, holdout_seed, "
+            "exclude_stations. When set, adds station_group stratification (all/holdout/holdin) "
+            "to all models and baselines using the truth dataset's station list."
+        ),
+    )
     parser.add_argument(
         "--member",
         type=str,
diff --git a/workflow/scripts/verification_plot_metrics.py b/workflow/scripts/verification_plot_metrics.py
index 64d316e1..8b9e19fe 100644
--- a/workflow/scripts/verification_plot_metrics.py
+++ b/workflow/scripts/verification_plot_metrics.py
@@ -93,6 +93,8 @@ def main(args: Namespace) -> None:
 
     # remove duplicated but not identical values from analyses (rounding errors)
     dfs = [xr.open_dataset(f) for f in args.verif_files]
+    # verify() always carries a station_group dim; select "all" for standard plots.
+    dfs = [d.sel(station_group="all", drop=True) for d in dfs]
     # 1) Ensure each dataset has unique lead_time values
     dfs = [_ensure_unique_lead_time(d) for d in dfs]
     # 2) For sources present in multiple datasets, keep the one with most lead_times
diff --git a/workflow/tools/config.schema.json b/workflow/tools/config.schema.json
index ea2b84ba..2cd8e7e0 100644
--- a/workflow/tools/config.schema.json
+++ b/workflow/tools/config.schema.json
@@ -296,6 +296,25 @@
           "$ref": "#/$defs/Dashboard",
           "description": "Settings for the experiment dashboard."
         },
+        "station_holdout": {
+          "anyOf": [
+            {
+              "items": {
+                "type": "string"
+              },
+              "type": "array"
+            },
+            {
+              "$ref": "#/$defs/StationHoldoutConfig"
+            },
+            {
+              "type": "null"
+            }
+          ],
+          "default": null,
+          "description": "When set, adds a 'station_group' dimension (all/holdout/holdin) to verification metrics for all models and baselines. Either a bare list of station nat_abbr (shorthand for exclude_stations), or a full object specifying holdout_fraction (+ optional holdout_seed) instead.",
+          "title": "Station Holdout"
+        },
         "scorecards": {
           "anyOf": [
             {
@@ -865,6 +884,50 @@
       "title": "ShowcaseConfig",
       "type": "object"
     },
+    "StationHoldoutConfig": {
+      "description": "Station holdout settings for station-group stratified verification.",
+      "properties": {
+        "holdout_fraction": {
+          "anyOf": [
+            {
+              "exclusiveMaximum": 1,
+              "exclusiveMinimum": 0,
+              "type": "number"
+            },
+            {
+              "type": "null"
+            }
+          ],
+          "default": null,
+          "description": "Fraction of truth stations to hold out for evaluation (exclusive 0\u20131). Mutually exclusive with exclude_stations.",
+          "title": "Holdout Fraction"
+        },
+        "holdout_seed": {
+          "default": 42,
+          "description": "Random seed for reproducible holdout station selection.",
+          "title": "Holdout Seed",
+          "type": "integer"
+        },
+        "exclude_stations": {
+          "anyOf": [
+            {
+              "items": {
+                "type": "string"
+              },
+              "type": "array"
+            },
+            {
+              "type": "null"
+            }
+          ],
+          "default": null,
+          "description": "Explicit list of station nat_abbr to hold out. Mutually exclusive with holdout_fraction.",
+          "title": "Exclude Stations"
+        }
+      },
+      "title": "StationHoldoutConfig",
+      "type": "object"
+    },
     "Stratification": {
       "description": "Stratification settings for the analysis.",
       "properties": {