diff --git a/README.md b/README.md
index 9d83065..134198b 100644
--- a/README.md
+++ b/README.md
@@ -31,9 +31,17 @@ uv run ccoa --help
 
 `ccoa classify-warc` streams WARC files from S3 (or any fsspec URL),
 extracts plain text from each response record with trafilatura, and
-applies a HuggingFace-hosted fasttext classifier. Per-record output is a
-CSV `URL,prediction_score,warc_filename,warc_record_index`; a one-shot score-distribution summary is
-logged at the end and written to a `<output>.summary.csv` file.
+applies one or more HuggingFace-hosted fasttext classifiers in a single
+pass. Per-record output is a CSV with one `score_<label>` column per
+requested label, between `URL` and the `warc_filename`/`warc_record_index`
+tail:
+
+```
+URL,score_<label_1>,...,score_<label_N>,warc_filename,warc_record_index
+```
+
+A per-column score-distribution summary is logged at the end and written
+to a `<output>.summary.csv` file.
 
 ```bash
 uv run ccoa classify-warc \
@@ -86,10 +94,37 @@ HTTPS gateway URL
 (`https://data.commoncrawl.org/...`) with no credentials.
 
 The default classifier is
-[`ibm-granite/GneissWeb.Sci_classifier`](https://huggingface.co/ibm-granite/GneissWeb.Sci_classifier)
-(`__label__science`). The first run downloads the ~4 GB model into the
-HuggingFace cache. Override with `--model-repo`, `--model-file`, and
-`--target-label`.
+[`ibm-granite/GneissWeb.Sci_classifier`](https://huggingface.co/ibm-granite/GneissWeb.Sci_classifier).
+Without `--labels` it emits both of the model's labels —
+`score___label__science` and `score___label__cc`, which sum to 1.0 per
+record. The first run downloads the ~4 GB model into the HuggingFace
+cache. Override with `--model-repo`, `--model-file`, and `--labels`.
+
+`--model-repo` and `--model-file` are list-valued and zipped positionally,
+so you can score against multiple classifiers in one pass:
+
+```bash
+uv run ccoa classify-warc \
+  --warc-paths 's3://commoncrawl/.../*.warc.gz' \
+  --model-repo ibm-granite/GneissWeb.Sci_classifier ibm-granite/GneissWeb.Quality_annotator \
+  --model-file fasttext_science.bin <quality_model_filename.bin> \
+  --output data/classified.csv
+```
+
+`--labels` is also list-valued (one entry per model). Each entry is a
+comma-separated list of labels (`"__label__science,__label__cc"`) or the
+literal `*` to use all of that model's labels (the default when `--labels`
+is omitted). Output columns are emitted in the order: models in CLI order,
+labels in the order given (or model-internal order for `*`).
+
+Column naming depends on whether the run has one model or many:
+
+- **Single model**: `score_<label>` (e.g. `score___label__science`).
+- **Multiple models**: `score_m<idx>_<label>`, where `<idx>` is the 0-based
+  CLI position of the model (e.g. `score_m0___label__science`,
+  `score_m1___label__hq`). This namespacing means two models can share a
+  label name — Sci_classifier and Quality_annotator both emit
+  `__label__cc` — without colliding.
 
 `--output` accepts `-` for stdout, any local path, or any fsspec URL —
 including `s3://bucket/key.csv`. S3 outputs use the same `--anonymous-s3`
@@ -133,11 +168,14 @@ uv run ccoa classify-warc \
   --output data/classified__resume-2.csv
 ```
 
-The resume CSV must include the `warc_filename` and `warc_record_index`
-columns (the current output schema always does). Records matching that
-`(warc_filename, record_index)` pair are skipped on the new run; the
-new `--output` contains only the missing rows. Concatenate the two
-CSVs (drop the second header) to get a complete result.
+The resume CSV's header must match the new run's output schema
+**exactly** — same `score_<label>` columns in the same order, between
+the leading `URL` and the trailing `warc_filename`/`warc_record_index`.
+Any drift (reorder, missing, extra) is rejected fast with a structured
+diff so a concatenation (drop the second header) yields a well-formed
+CSV. Records matching that `(warc_filename, record_index)` pair are
+skipped on the new run; the new `--output` contains only the missing
+rows.
 
 With `--records-per-file-limit N` the limit is interpreted as the
 **target total** per file (resumed + new). Files already at the target
diff --git a/notebooks/compare_classifier_scores.ipynb b/notebooks/compare_classifier_scores.ipynb
index b3558d2..caa0a5a 100644
--- a/notebooks/compare_classifier_scores.ipynb
+++ b/notebooks/compare_classifier_scores.ipynb
@@ -35,10 +35,10 @@
    "id": "cell-01-imports",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:21.594052Z",
-     "iopub.status.busy": "2026-05-26T10:15:21.593943Z",
-     "iopub.status.idle": "2026-05-26T10:15:41.227039Z",
-     "shell.execute_reply": "2026-05-26T10:15:41.226661Z"
+     "iopub.execute_input": "2026-05-28T19:40:08.502488Z",
+     "iopub.status.busy": "2026-05-28T19:40:08.502328Z",
+     "iopub.status.idle": "2026-05-28T19:40:19.488256Z",
+     "shell.execute_reply": "2026-05-28T19:40:19.487890Z"
     }
    },
    "outputs": [
@@ -98,7 +98,7 @@
    "source": [
     "## 1. Load the score data\n",
     "\n",
-    "The output CSV schema is `URL,prediction_score,warc_filename,warc_record_index`; we only need `prediction_score` for distribution analysis.\n",
+    "The output CSV schema is `URL,score_<label_1>,...,score_<label_N>,warc_filename,warc_record_index` — one score column per model label. The cell below auto-detects available `score_*` columns and picks the first as `SCORE_COL`; change it (or extend `load_scores`) to compare other labels. Older single-column outputs (`prediction_score`) are accepted as a fallback so existing CSVs in `data/` still work.\n",
     "\n",
     "The focus crawl can be split across **multiple parts** when `--resume-from-output` was used to recover from a crash (e.g. `…_1m_scores.csv` for the first run plus `…_1m_scores__resume-2.csv` for the continuation). We auto-discover all parts matching `{focus name}_1m_scores*.csv` (filtering out the `.summary.csv` sidecars) and concatenate them into a single distribution.\n",
     "\n",
@@ -113,13 +113,21 @@
    "id": "cell-03-load",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:41.228341Z",
-     "iopub.status.busy": "2026-05-26T10:15:41.228232Z",
-     "iopub.status.idle": "2026-05-26T10:15:41.973842Z",
-     "shell.execute_reply": "2026-05-26T10:15:41.973414Z"
+     "iopub.execute_input": "2026-05-28T19:40:19.489927Z",
+     "iopub.status.busy": "2026-05-28T19:40:19.489801Z",
+     "iopub.status.idle": "2026-05-28T19:40:20.224124Z",
+     "shell.execute_reply": "2026-05-28T19:40:20.223694Z"
     }
    },
    "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Score columns in baseline: ['prediction_score (legacy)']\n",
+      "Analysing: prediction_score\n"
+     ]
+    },
     {
      "name": "stdout",
      "output_type": "stream",
@@ -136,9 +144,19 @@
     }
    ],
    "source": [
+    "# Auto-detect the score column to analyse. classify-warc now emits one\n",
+    "# `score_<label>` column per model label; older outputs had a single\n",
+    "# `prediction_score`. Edit SCORE_COL to compare other labels.\n",
+    "_baseline_columns = pd.read_csv(BASELINE_CSV, nrows=0).columns.tolist()\n",
+    "SCORE_COLUMNS = [c for c in _baseline_columns if c.startswith(\"score_\")]\n",
+    "SCORE_COL = SCORE_COLUMNS[0] if SCORE_COLUMNS else \"prediction_score\"\n",
+    "print(f\"Score columns in baseline: {SCORE_COLUMNS or [SCORE_COL + ' (legacy)']}\")\n",
+    "print(f\"Analysing: {SCORE_COL}\")\n",
+    "\n",
+    "\n",
     "def load_scores(csv_path: Path) -> pd.Series:\n",
-    "    \"\"\"Read the prediction_score column from a classify-warc output CSV.\"\"\"\n",
-    "    series = pd.read_csv(csv_path, usecols=[\"prediction_score\"])[\"prediction_score\"]\n",
+    "    \"\"\"Read the `SCORE_COL` column from a classify-warc output CSV.\"\"\"\n",
+    "    series = pd.read_csv(csv_path, usecols=[SCORE_COL])[SCORE_COL]\n",
     "    return series.astype(float)\n",
     "\n",
     "\n",
@@ -151,7 +169,9 @@
     "    is_toy = False\n",
     "    ratio = len(focus_scores) / max(len(baseline_scores), 1)\n",
     "    is_partial = ratio < 0.95\n",
-    "    per_part = \", \".join(f\"{p.name}: {len(s):,}\" for p, s in zip(FOCUS_CSVS, parts))\n",
+    "    per_part = \", \".join(\n",
+    "        f\"{p.name}: {len(s):,}\" for p, s in zip(FOCUS_CSVS, parts, strict=True)\n",
+    "    )\n",
     "    print(\n",
     "        f\"{FOCUS_NAME}: {len(focus_scores):,} scores loaded from {len(FOCUS_CSVS)} \"\n",
     "        f\"file(s) ({ratio:.0%} of baseline size). Parts → {per_part}\"\n",
@@ -164,7 +184,7 @@
     "else:\n",
     "    rng = np.random.default_rng(42)\n",
     "    noisy = baseline_scores.to_numpy() + 0.05 + rng.normal(0.0, 0.10, size=len(baseline_scores))\n",
-    "    focus_scores = pd.Series(np.clip(noisy, 0.0, 1.0), name=\"prediction_score\")\n",
+    "    focus_scores = pd.Series(np.clip(noisy, 0.0, 1.0), name=SCORE_COL)\n",
     "    is_toy = True\n",
     "    is_partial = False\n",
     "    print(\n",
@@ -195,10 +215,10 @@
    "id": "cell-05-stats",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:41.974946Z",
-     "iopub.status.busy": "2026-05-26T10:15:41.974894Z",
-     "iopub.status.idle": "2026-05-26T10:15:42.081238Z",
-     "shell.execute_reply": "2026-05-26T10:15:42.080957Z"
+     "iopub.execute_input": "2026-05-28T19:40:20.225163Z",
+     "iopub.status.busy": "2026-05-28T19:40:20.225105Z",
+     "iopub.status.idle": "2026-05-28T19:40:20.333180Z",
+     "shell.execute_reply": "2026-05-28T19:40:20.332861Z"
     }
    },
    "outputs": [
@@ -381,10 +401,10 @@
    "id": "cell-07-hist",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:42.082323Z",
-     "iopub.status.busy": "2026-05-26T10:15:42.082276Z",
-     "iopub.status.idle": "2026-05-26T10:15:42.642052Z",
-     "shell.execute_reply": "2026-05-26T10:15:42.641716Z"
+     "iopub.execute_input": "2026-05-28T19:40:20.334135Z",
+     "iopub.status.busy": "2026-05-28T19:40:20.334086Z",
+     "iopub.status.idle": "2026-05-28T19:40:20.602002Z",
+     "shell.execute_reply": "2026-05-28T19:40:20.601654Z"
     }
    },
    "outputs": [
@@ -403,7 +423,7 @@
     "fig, axes = plt.subplots(1, 2, figsize=(14, 5), sharex=True)\n",
     "bins = np.linspace(0.0, 1.0, 81)\n",
     "\n",
-    "for ax, log_scale in zip(axes, [False, True]):\n",
+    "for ax, log_scale in zip(axes, [False, True], strict=True):\n",
     "    ax.hist(baseline_scores, bins=bins, alpha=0.5, label=BASELINE_NAME, color=\"C0\", density=True)\n",
     "    ax.hist(focus_scores, bins=bins, alpha=0.5, label=FOCUS_NAME, color=\"C1\", density=True)\n",
     "    if log_scale:\n",
@@ -411,7 +431,7 @@
     "        ax.set_title(\"Histogram (log y, density-normalised)\")\n",
     "    else:\n",
     "        ax.set_title(\"Histogram (linear y, density-normalised)\")\n",
-    "    ax.set_xlabel(\"prediction_score\")\n",
+    "    ax.set_xlabel(SCORE_COL)\n",
     "    ax.set_ylabel(\"density\")\n",
     "    ax.legend()\n",
     "    ax.grid(True, alpha=0.3)\n",
@@ -436,10 +456,10 @@
    "id": "cell-09-ecdf",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:42.643283Z",
-     "iopub.status.busy": "2026-05-26T10:15:42.643203Z",
-     "iopub.status.idle": "2026-05-26T10:15:43.228237Z",
-     "shell.execute_reply": "2026-05-26T10:15:43.227857Z"
+     "iopub.execute_input": "2026-05-28T19:40:20.602921Z",
+     "iopub.status.busy": "2026-05-28T19:40:20.602874Z",
+     "iopub.status.idle": "2026-05-28T19:40:21.206119Z",
+     "shell.execute_reply": "2026-05-28T19:40:21.205709Z"
     }
    },
    "outputs": [
@@ -469,7 +489,7 @@
     "]:\n",
     "    x, y = ecdf(series)\n",
     "    plt.plot(x, y, label=name, color=color)\n",
-    "plt.xlabel(\"prediction_score\")\n",
+    "plt.xlabel(SCORE_COL)\n",
     "plt.ylabel(\"ECDF P(score ≤ x)\")\n",
     "plt.title(\"Empirical CDF of classifier scores\")\n",
     "plt.legend()\n",
@@ -494,10 +514,10 @@
    "id": "cell-11-box",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:43.229226Z",
-     "iopub.status.busy": "2026-05-26T10:15:43.229167Z",
-     "iopub.status.idle": "2026-05-26T10:15:43.929555Z",
-     "shell.execute_reply": "2026-05-26T10:15:43.929175Z"
+     "iopub.execute_input": "2026-05-28T19:40:21.207200Z",
+     "iopub.status.busy": "2026-05-28T19:40:21.207147Z",
+     "iopub.status.idle": "2026-05-28T19:40:21.891974Z",
+     "shell.execute_reply": "2026-05-28T19:40:21.891665Z"
     }
    },
    "outputs": [
@@ -519,14 +539,14 @@
     "labels = [BASELINE_NAME, FOCUS_NAME]\n",
     "\n",
     "axes[0].boxplot(data, tick_labels=labels, showfliers=False)\n",
-    "axes[0].set_ylabel(\"prediction_score\")\n",
+    "axes[0].set_ylabel(SCORE_COL)\n",
     "axes[0].set_title(\"Box plot (outliers hidden)\")\n",
     "axes[0].grid(True, alpha=0.3)\n",
     "\n",
     "axes[1].violinplot(data, showmedians=True)\n",
     "axes[1].set_xticks([1, 2])\n",
     "axes[1].set_xticklabels(labels)\n",
-    "axes[1].set_ylabel(\"prediction_score\")\n",
+    "axes[1].set_ylabel(SCORE_COL)\n",
     "axes[1].set_title(\"Violin plot\")\n",
     "axes[1].grid(True, alpha=0.3)\n",
     "\n",
@@ -552,10 +572,10 @@
    "id": "cell-13-qq",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:43.930547Z",
-     "iopub.status.busy": "2026-05-26T10:15:43.930496Z",
-     "iopub.status.idle": "2026-05-26T10:15:44.006062Z",
-     "shell.execute_reply": "2026-05-26T10:15:44.005773Z"
+     "iopub.execute_input": "2026-05-28T19:40:21.893014Z",
+     "iopub.status.busy": "2026-05-28T19:40:21.892965Z",
+     "iopub.status.idle": "2026-05-28T19:40:21.973772Z",
+     "shell.execute_reply": "2026-05-28T19:40:21.973431Z"
     }
    },
    "outputs": [
@@ -611,10 +631,10 @@
    "id": "cell-15-tests",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:44.007069Z",
-     "iopub.status.busy": "2026-05-26T10:15:44.007022Z",
-     "iopub.status.idle": "2026-05-26T10:15:44.240328Z",
-     "shell.execute_reply": "2026-05-26T10:15:44.239994Z"
+     "iopub.execute_input": "2026-05-28T19:40:21.974787Z",
+     "iopub.status.busy": "2026-05-28T19:40:21.974737Z",
+     "iopub.status.idle": "2026-05-28T19:40:22.210867Z",
+     "shell.execute_reply": "2026-05-28T19:40:22.210554Z"
     }
    },
    "outputs": [
@@ -714,10 +734,10 @@
    "id": "cell-17-effect",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:44.241419Z",
-     "iopub.status.busy": "2026-05-26T10:15:44.241371Z",
-     "iopub.status.idle": "2026-05-26T10:15:44.262566Z",
-     "shell.execute_reply": "2026-05-26T10:15:44.262191Z"
+     "iopub.execute_input": "2026-05-28T19:40:22.211857Z",
+     "iopub.status.busy": "2026-05-28T19:40:22.211810Z",
+     "iopub.status.idle": "2026-05-28T19:40:22.233262Z",
+     "shell.execute_reply": "2026-05-28T19:40:22.232953Z"
     }
    },
    "outputs": [
@@ -831,10 +851,10 @@
    "id": "cell-19-verdict",
    "metadata": {
     "execution": {
-     "iopub.execute_input": "2026-05-26T10:15:44.263522Z",
-     "iopub.status.busy": "2026-05-26T10:15:44.263466Z",
-     "iopub.status.idle": "2026-05-26T10:15:44.283675Z",
-     "shell.execute_reply": "2026-05-26T10:15:44.283399Z"
+     "iopub.execute_input": "2026-05-28T19:40:22.234289Z",
+     "iopub.status.busy": "2026-05-28T19:40:22.234236Z",
+     "iopub.status.idle": "2026-05-28T19:40:22.254388Z",
+     "shell.execute_reply": "2026-05-28T19:40:22.254084Z"
     }
    },
    "outputs": [
diff --git a/src/ccoa/classifier/fasttext.py b/src/ccoa/classifier/fasttext.py
index bc22c28..7ce3825 100644
--- a/src/ccoa/classifier/fasttext.py
+++ b/src/ccoa/classifier/fasttext.py
@@ -2,12 +2,13 @@
 
 from __future__ import annotations
 
+from collections.abc import Sequence
+
 import fasttext as _fasttext
 from huggingface_hub import hf_hub_download
 
 DEFAULT_MODEL_REPO = "ibm-granite/GneissWeb.Sci_classifier"
 DEFAULT_MODEL_FILE = "fasttext_science.bin"
-DEFAULT_TARGET_LABEL = "__label__science"
 
 FASTTEXT_MAX_INPUT_CHARS = 100_000
 
@@ -36,17 +37,19 @@ def clean_for_fasttext(text: str, max_len: int = FASTTEXT_MAX_INPUT_CHARS) -> st
     return cleaned
 
 
-def predict_target(model, text: str, target_label: str) -> float:
-    """Return the probability the model assigns to `target_label`.
+def predict_targets(model, text: str, target_labels: Sequence[str]) -> list[float]:
+    """Return probabilities for `target_labels` in the order given.
 
-    Cleans the input via `clean_for_fasttext` (strip newlines + NULs, clamp
-    length), then asks the model for probabilities over all labels and
-    returns the probability of `target_label` (0.0 if the model emits no
-    such label).
+    `model.predict(..., k=-1)` returns labels in descending probability order
+    (data-dependent), so we resolve scores by label name. Missing labels
+    yield 0.0.
     """
     cleaned = clean_for_fasttext(text)
     labels, probs = model.predict(cleaned, k=-1)
-    for lbl, prob in zip(labels, probs, strict=False):
-        if lbl == target_label:
-            return float(prob)
-    return 0.0
+    label_to_prob = dict(zip(labels, probs, strict=False))
+    return [float(label_to_prob.get(lbl, 0.0)) for lbl in target_labels]
+
+
+def get_model_labels(model) -> tuple[str, ...]:
+    """Return all labels the model can emit, in model-internal order."""
+    return tuple(model.get_labels())
diff --git a/src/ccoa/commands/classify_warc.py b/src/ccoa/commands/classify_warc.py
index 2e67d0e..067573d 100644
--- a/src/ccoa/commands/classify_warc.py
+++ b/src/ccoa/commands/classify_warc.py
@@ -1,12 +1,13 @@
-"""Classify WARC response records with a fasttext model.
+"""Classify WARC response records with one or more fasttext models.
 
-Streams one or more WARC files from S3 (or any fsspec URL), extracts
-plain text from each `response` record with trafilatura, and applies
-a HuggingFace-hosted fasttext classifier. Per-record output is a CSV
-`URL,prediction_score,warc_filename,warc_record_index` written to
-stdout, a local path, or an `s3://` URL (anything `fsspec.open`
-understands); a one-shot summary of the score distribution is logged
-at the end and written to a `<output>.summary.csv` sidecar.
+Streams WARC files from S3 (or any fsspec URL), extracts plain text from
+each `response` record with trafilatura, and applies one or more
+HuggingFace-hosted fasttext classifiers in a single pass. Per-record
+output is a CSV `URL,score_<label_1>,...,score_<label_N>,warc_filename,
+warc_record_index` written to stdout, a local path, or an `s3://` URL
+(anything `fsspec.open` understands). A one-shot per-column summary of
+each score distribution is logged at the end and written to a
+`<output>.summary.csv` sidecar.
 
 Example:
     ```bash
@@ -19,7 +20,9 @@
     ```
 
 The default model (`ibm-granite/GneissWeb.Sci_classifier`) is ~4 GB;
-the first run will download it into the HuggingFace cache.
+the first run will download it into the HuggingFace cache. Without
+`--labels`, all labels of each model are emitted (the science model
+contributes `score___label__science` and `score___label__cc`).
 """
 
 from __future__ import annotations
@@ -40,9 +43,9 @@
 from ccoa.classifier.fasttext import (
     DEFAULT_MODEL_FILE,
     DEFAULT_MODEL_REPO,
-    DEFAULT_TARGET_LABEL,
+    get_model_labels,
     load_classifier,
-    predict_target,
+    predict_targets,
 )
 from ccoa.commands import BaseCommand
 from ccoa.extraction.cache import (
@@ -70,33 +73,98 @@
 logger = logging.getLogger(__name__)
 
 
-_WORKER_MODEL: object = None
+# Per-process state. In thread mode, the parent populates `_WORKER_MODELS`
+# directly. In process mode, each worker's initializer rebuilds it from the
+# spec the parent pickled into the pool's `initargs`.
+_WORKER_MODELS: list[tuple[object, list[str]]] = []
 _WORKER_LOCK = threading.Lock()
 
 
+_LABELS_ALL = "*"
+
+
+def _resolve_model_slots(
+    model_repos: list[str],
+    model_files: list[str],
+    labels: list[str],
+) -> list[tuple[str, str, list[str]]]:
+    """Pair `--model-repo`/`--model-file`/`--labels` into `(repo, file, label_spec)` slots.
+
+    `label_spec` is either `["*"]` (resolve to all of the model's labels at
+    load time) or an explicit ordered list. Raises `ValueError` for mismatched
+    list lengths or empty CLI input without a default fallback.
+    """
+    if not model_repos and not model_files:
+        model_repos = [DEFAULT_MODEL_REPO]
+        model_files = [DEFAULT_MODEL_FILE]
+    if len(model_repos) != len(model_files):
+        raise ValueError(
+            f"--model-repo and --model-file must have the same length; "
+            f"got {len(model_repos)} repo(s) and {len(model_files)} file(s)."
+        )
+    if not labels:
+        labels = [_LABELS_ALL] * len(model_repos)
+    if len(labels) != len(model_repos):
+        raise ValueError(
+            f"--labels must have one entry per --model-repo (or be omitted "
+            f"to default to '*' for every model); got {len(labels)} entries "
+            f"for {len(model_repos)} model(s)."
+        )
+
+    slots: list[tuple[str, str, list[str]]] = []
+    for repo, file, label_spec in zip(model_repos, model_files, labels, strict=True):
+        if label_spec == _LABELS_ALL:
+            resolved = [_LABELS_ALL]
+        else:
+            resolved = [s.strip() for s in label_spec.split(",") if s.strip()]
+            if not resolved:
+                raise ValueError(
+                    f"--labels entry for model {repo}/{file} is empty; "
+                    f"use '*' (all labels) or a comma-separated list."
+                )
+        slots.append((repo, file, resolved))
+    return slots
+
+
 def load_resume_skipset(
     path: str,
     storage_options: dict[str, object],
+    expected_header: list[str],
 ) -> dict[str, frozenset[int]]:
     """Build `{warc_filename: frozenset(record_indices)}` from a prior output CSV.
 
     Used by `--resume-from-output` to skip records already classified. The
-    input must have `warc_filename` and `warc_record_index` columns (which
-    the current output schema always emits); a `ValueError` is raised
-    otherwise to flag stale outputs from the pre-resume schema.
+    file's header MUST equal `expected_header` exactly — same columns in the
+    same order — so a concatenation of the prior CSV and the new `--output`
+    yields a well-formed file. Any drift (added columns, removed columns,
+    reorder) raises `ValueError` with a structured diff.
     """
-    required = {"warc_filename", "warc_record_index"}
-    buckets: dict[str, set[int]] = {}
     with fsspec.open(path, "r", encoding="utf-8", **storage_options) as fh:
         reader = csv.DictReader(fh)
-        fields = set(reader.fieldnames or [])
-        if not required.issubset(fields):
-            missing = required - fields
+        actual = list(reader.fieldnames or [])
+        if actual != expected_header:
+            actual_set = set(actual)
+            expected_set = set(expected_header)
+            missing = expected_set - actual_set
+            extra = actual_set - expected_set
+            score_cols = [c for c in actual if c.startswith("score_")]
+            if not score_cols:
+                raise ValueError(
+                    f"{path} has no `score_*` columns; this looks like a "
+                    f"pre-multi-label CSV and is not resumable with the "
+                    f"current --model-repo/--model-file/--labels selection. "
+                    f"Found header: {actual}; expected: {expected_header}."
+                )
             raise ValueError(
-                f"{path} is missing required columns for --resume-from-output: "
-                f"need {sorted(required)}, missing {sorted(missing)} "
-                f"(found {sorted(fields)})."
+                f"{path} header does not match the planned output schema. "
+                f"Pass --model-repo/--model-file/--labels so the new run "
+                f"produces the same columns in the same order.\n"
+                f"  expected: {expected_header}\n"
+                f"  actual:   {actual}\n"
+                f"  missing:  {sorted(missing)}\n"
+                f"  extra:    {sorted(extra)}"
             )
+        buckets: dict[str, set[int]] = {}
         for row in reader:
             warc = row["warc_filename"]
             try:
@@ -111,12 +179,13 @@ def load_resume_skipset(
 class FileResult:
     """Aggregated per-file outcome returned by `process_one_file`.
 
-    Each entry in `rows` is `(url, score, warc_filename, warc_record_index)`,
-    matching the output CSV schema.
+    Each entry in `rows` is `(url, scores, warc_filename, warc_record_index)`,
+    where `scores` is a list of probabilities matching `args.score_columns`
+    in order (one entry per `score_<label>` column).
     """
 
     uri: str
-    rows: list[tuple[str, float, str, int]]
+    rows: list[tuple[str, list[float], str, int]]
     processed: int
     skipped_empty: int
     skipped_homepage: int
@@ -130,7 +199,7 @@ class FileResult:
 
 def process_one_file(
     uri: str,
-    model,
+    models_spec: list[tuple[object, list[str]]],
     model_lock: threading.Lock,
     args: argparse.Namespace,
     skip_indices: frozenset[int] = frozenset(),
@@ -139,11 +208,11 @@ def process_one_file(
 
     Owns its own fsspec input stream, its own per-WARC cache load/save
     (when `--cache-dir` is set), and its own `--records-per-file-limit`
-    counter. Calls `predict_target` under `model_lock` since fasttext
-    makes no thread-safety guarantees; the lock cost is sub-ms
-    compared to the per-record extract cost.
+    counter. Holds `model_lock` across every `predict_targets` call (one
+    per model) since fasttext makes no thread-safety guarantees; the lock
+    cost is sub-ms compared to the per-record extract cost.
     """
-    rows: list[tuple[str, float, str, int]] = []
+    rows: list[tuple[str, list[float], str, int]] = []
     processed = 0
     skipped_empty = 0
     skipped_homepage = 0
@@ -245,11 +314,13 @@ def process_one_file(
                     continue
 
                 t0 = time.perf_counter()
+                scores: list[float] = []
                 with model_lock:
-                    score = predict_target(model, cleaned, args.target_label)
+                    for model, target_labels in models_spec:
+                        scores.extend(predict_targets(model, cleaned, target_labels))
                 t_predict += time.perf_counter() - t0
 
-                rows.append((url, score, uri, record_index))
+                rows.append((url, scores, uri, record_index))
                 processed += 1
 
                 if progress_every and processed % progress_every == 0:
@@ -293,15 +364,21 @@ def process_one_file(
     )
 
 
-def _process_pool_initializer(model_repo: str, model_file: str) -> None:
-    """Per-worker setup for `ProcessPoolExecutor`: load the model once.
+def _process_pool_initializer(
+    spec: list[tuple[str, str, list[str]]],
+) -> None:
+    """Per-worker setup for `ProcessPoolExecutor`: load every model once.
+
+    `spec` is `[(repo, file, resolved_labels), ...]` in CLI order. The parent
+    resolved any `'*'` placeholders to actual label lists before pickling, so
+    workers reuse that exact label vocabulary and column ordering.
 
     Also silences trafilatura's chatty per-page WARN/ERROR logs in this
     worker. Each worker process has its own lxml + fasttext state, so a
     heap-corruption abort here cannot tear down the parent or its peers.
     """
-    global _WORKER_MODEL  # noqa: PLW0603
-    _WORKER_MODEL = load_classifier(model_repo, model_file)
+    global _WORKER_MODELS  # noqa: PLW0603
+    _WORKER_MODELS = [(load_classifier(repo, file), labels) for repo, file, labels in spec]
     logging.getLogger("trafilatura").setLevel(logging.CRITICAL)
 
 
@@ -310,7 +387,7 @@ def _process_pool_worker(
 ) -> FileResult:
     """Top-level pickleable adapter that runs `process_one_file` in a worker."""
     uri, args, skip_indices = payload
-    return process_one_file(uri, _WORKER_MODEL, _WORKER_LOCK, args, skip_indices)
+    return process_one_file(uri, _WORKER_MODELS, _WORKER_LOCK, args, skip_indices)
 
 
 class ClassifyWarcCommand(BaseCommand):
@@ -382,20 +459,39 @@ def add_arguments(self, parser: argparse.ArgumentParser) -> None:
         )
         parser.add_argument(
             "--model-repo",
-            default=DEFAULT_MODEL_REPO,
-            help=f"HuggingFace repo id of the fasttext model (default: {DEFAULT_MODEL_REPO}).",
+            nargs="*",
+            default=[],
+            metavar="REPO",
+            help=(
+                "HuggingFace repo id(s) of the fasttext model(s). Repeatable; "
+                "must have the same length as --model-file. When both --model-repo "
+                "and --model-file are omitted, falls back to "
+                f"{DEFAULT_MODEL_REPO}/{DEFAULT_MODEL_FILE}."
+            ),
         )
         parser.add_argument(
             "--model-file",
-            default=DEFAULT_MODEL_FILE,
-            help=f"Filename of the .bin model inside the repo (default: {DEFAULT_MODEL_FILE}).",
+            nargs="*",
+            default=[],
+            metavar="FILE",
+            help=(
+                "Filename(s) of the .bin model inside each repo, positionally "
+                "paired with --model-repo (same length)."
+            ),
         )
         parser.add_argument(
-            "--target-label",
-            default=DEFAULT_TARGET_LABEL,
+            "--labels",
+            nargs="*",
+            default=[],
+            metavar="LABELS",
             help=(
-                "Label whose probability is written as 'prediction_score' "
-                f"(default: {DEFAULT_TARGET_LABEL})."
+                "Per-model label filter, positionally paired with --model-repo. "
+                "Each entry is a comma-separated list of fasttext labels "
+                "(e.g. '__label__science,__label__cc'), or the literal '*' to "
+                "use all labels of that model (the default when --labels is "
+                "omitted entirely). Output columns are `score_<label>` in the "
+                "order: models in CLI order, labels in the order given (or in "
+                "model-internal order for '*')."
             ),
         )
         parser.add_argument(
@@ -481,6 +577,14 @@ def add_arguments(self, parser: argparse.ArgumentParser) -> None:
                 "throughput this is roughly one line per ~40s of CPU work."
             ),
         )
+        parser.add_argument(
+            "--overwrite",
+            action="store_true",
+            help=(
+                "Overwrite an existing --output (and its `.summary` sidecar) "
+                "instead of failing fast. Off by default to protect prior runs."
+            ),
+        )
         parser.add_argument(
             "--anonymous-s3",
             action="store_true",
@@ -507,6 +611,12 @@ def run(self, args: argparse.Namespace) -> int:
             )
             return 2
 
+        try:
+            model_slots = _resolve_model_slots(args.model_repo, args.model_file, args.labels)
+        except ValueError as exc:
+            logger.error("%s", exc)
+            return 2
+
         output_storage_options = s3_storage_options(
             args.output, args.anonymous_s3, args.s3_requester_pays
         )
@@ -522,12 +632,16 @@ def run(self, args: argparse.Namespace) -> int:
                 (summary_uri, summary_storage_options),
             ):
                 if output_exists(path, opts):
-                    logger.error(
-                        "Output already exists: %s. Refusing to overwrite; "
-                        "delete the file(s) or pick a fresh --output path.",
-                        path,
-                    )
-                    return 2
+                    if args.overwrite:
+                        logger.warning("Overwriting existing output: %s", path)
+                    else:
+                        logger.error(
+                            "Output already exists: %s. Refusing to overwrite; "
+                            "delete the file(s), pick a fresh --output path, or "
+                            "pass --overwrite to replace them.",
+                            path,
+                        )
+                        return 2
 
         logger.info("Resolving WARC paths: %s", args.warc_paths)
 
@@ -550,16 +664,70 @@ def run(self, args: argparse.Namespace) -> int:
 
         if not resolved:
             logger.warning("No WARC files matched --warc-paths; nothing to classify.")
-            log_summary([])
+            log_summary({})
             return 0
 
+        # Load each model once in the parent so we can resolve `*` label slots
+        # via `get_model_labels` and freeze the global column order BEFORE any
+        # workers start. In thread mode and --workers 1 the parent's models are
+        # reused for scoring; in process mode they're freed and each worker
+        # reloads them from the same `resolved_spec`.
+        runtime_models: list[tuple[object, list[str]]] = []
+        resolved_spec: list[tuple[str, str, list[str]]] = []
+        for repo, file, label_spec in model_slots:
+            logger.info("Loading model %s/%s", repo, file)
+            model_obj = load_classifier(repo, file)
+            if label_spec == [_LABELS_ALL]:
+                labels = list(get_model_labels(model_obj))
+            else:
+                labels = list(label_spec)
+            resolved_spec.append((repo, file, labels))
+            runtime_models.append((model_obj, labels))
+            logger.info("  resolved labels for %s/%s: %s", repo, file, labels)
+
+        # Single model: `score_<label>` keeps the CSV clean. Multiple models:
+        # `score_m<idx>_<label>` so each model's labels live in their own
+        # namespace (e.g. both Sci_classifier and Quality_annotator emit
+        # `__label__cc`, which would otherwise collide).
+        score_columns: list[str] = []
+        seen_columns: dict[str, tuple[str, str]] = {}
+        multi_model = len(resolved_spec) > 1
+        for idx, (repo, file, labels) in enumerate(resolved_spec):
+            prefix = f"score_m{idx}_" if multi_model else "score_"
+            for lbl in labels:
+                col = f"{prefix}{lbl}"
+                if col in seen_columns:
+                    prev_repo, prev_file = seen_columns[col]
+                    logger.error(
+                        "Column %s would come from both %s/%s and %s/%s. "
+                        "Same model+label listed twice?",
+                        col,
+                        prev_repo,
+                        prev_file,
+                        repo,
+                        file,
+                    )
+                    return 2
+                seen_columns[col] = (repo, file)
+                score_columns.append(col)
+        args.score_columns = score_columns
+        logger.info("Output score columns (%d): %s", len(score_columns), score_columns)
+
+        expected_header = ["URL", *score_columns, "warc_filename", "warc_record_index"]
+
         skip_by_warc: dict[str, frozenset[int]] = {}
         if args.resume_from_output:
             resume_storage_options = s3_storage_options(
                 args.resume_from_output, args.anonymous_s3, args.s3_requester_pays
             )
             logger.info("Loading resume skip-set from %s", args.resume_from_output)
-            skip_by_warc = load_resume_skipset(args.resume_from_output, resume_storage_options)
+            try:
+                skip_by_warc = load_resume_skipset(
+                    args.resume_from_output, resume_storage_options, expected_header
+                )
+            except ValueError as exc:
+                logger.error("%s", exc)
+                return 2
             total_skips = sum(len(v) for v in skip_by_warc.values())
             logger.info(
                 "Resume skip-set: %d records across %d distinct WARCs.",
@@ -568,20 +736,14 @@ def run(self, args: argparse.Namespace) -> int:
             )
 
         use_process_pool = args.workers > 1 and args.workers_mode == "process"
+        model_lock = threading.Lock()
         if use_process_pool:
             logger.info(
-                "Worker mode 'process'; deferring model load to each of "
-                "%d worker processes (model: %s/%s).",
+                "Worker mode 'process'; freeing parent-side models and deferring "
+                "load to each of %d worker processes.",
                 args.workers,
-                args.model_repo,
-                args.model_file,
             )
-            model = None
-            model_lock = threading.Lock()
-        else:
-            logger.info("Loading model %s/%s", args.model_repo, args.model_file)
-            model = load_classifier(args.model_repo, args.model_file)
-            model_lock = threading.Lock()
+            runtime_models = []
 
         # Silence trafilatura's per-page WARN/ERROR noise (empty / non-HTML /
         # bad-encoding pages) — these outcomes already land in `skipped_empty`.
@@ -592,7 +754,7 @@ def run(self, args: argparse.Namespace) -> int:
         traf_prev_propagate = traf_logger.propagate
         traf_logger.propagate = False
 
-        scores: list[float] = []
+        scores_by_column: dict[str, list[float]] = {col: [] for col in score_columns}
         processed = 0
         skipped_empty = 0
         skipped_homepage = 0
@@ -614,7 +776,7 @@ def run(self, args: argparse.Namespace) -> int:
         t_processing_start = time.perf_counter()
         with open_output_sink(args.output, output_storage_options) as sink:
             writer = csv.writer(sink)
-            writer.writerow(["URL", "prediction_score", "warc_filename", "warc_record_index"])
+            writer.writerow(expected_header)
 
             def _aggregate(result: FileResult) -> None:
                 """Write `result` rows to the CSV and fold its counters into the totals."""
@@ -622,9 +784,12 @@ def _aggregate(result: FileResult) -> None:
                 nonlocal skipped_resume_total, extract_errors_total
                 nonlocal cache_hits_total, cache_misses_total
                 nonlocal t_extract_total, t_predict_total, files_done
-                for url, score, warc_filename, record_index in result.rows:
-                    writer.writerow([url, f"{score:.6f}", warc_filename, record_index])
-                    scores.append(score)
+                for url, scores, warc_filename, record_index in result.rows:
+                    writer.writerow(
+                        [url, *(f"{s:.6f}" for s in scores), warc_filename, record_index]
+                    )
+                    for col, s in zip(score_columns, scores, strict=True):
+                        scores_by_column[col].append(s)
                 if sink is sys.stdout:
                     sink.flush()
                 processed += result.processed
@@ -663,7 +828,11 @@ def _aggregate(result: FileResult) -> None:
                 for uri in resolved:
                     _aggregate(
                         process_one_file(
-                            uri, model, model_lock, args, skip_by_warc.get(uri, empty_skip)
+                            uri,
+                            runtime_models,
+                            model_lock,
+                            args,
+                            skip_by_warc.get(uri, empty_skip),
                         )
                     )
                     if limit and processed >= limit:
@@ -679,7 +848,7 @@ def _aggregate(result: FileResult) -> None:
                         with ProcessPoolExecutor(
                             max_workers=args.workers,
                             initializer=_process_pool_initializer,
-                            initargs=(args.model_repo, args.model_file),
+                            initargs=(resolved_spec,),
                         ) as pool:
                             payloads = [
                                 (u, args, skip_by_warc.get(u, empty_skip)) for u in remaining_uris
@@ -723,7 +892,11 @@ def _aggregate(result: FileResult) -> None:
                 with ThreadPoolExecutor(max_workers=args.workers) as pool:
                     for result in pool.map(
                         lambda u: process_one_file(
-                            u, model, model_lock, args, skip_by_warc.get(u, empty_skip)
+                            u,
+                            runtime_models,
+                            model_lock,
+                            args,
+                            skip_by_warc.get(u, empty_skip),
                         ),
                         resolved,
                     ):
@@ -746,7 +919,7 @@ def _aggregate(result: FileResult) -> None:
                 cache_hits_total,
                 cache_misses_total,
             )
-        log_summary(scores)
+        log_summary(scores_by_column)
         log_timing(processed, t_processing, t_extract_total, t_predict_total)
 
         if traf_counter.warnings or traf_counter.errors:
@@ -763,7 +936,7 @@ def _aggregate(result: FileResult) -> None:
                 summary_storage_options,
                 args=args,
                 resolved_count=len(resolved),
-                scores=scores,
+                scores_by_column=scores_by_column,
                 processed=processed,
                 skipped_empty=skipped_empty,
                 skipped_homepage=skipped_homepage,
diff --git a/src/ccoa/utils/reporting/summary.py b/src/ccoa/utils/reporting/summary.py
index 8aafb00..05d1008 100644
--- a/src/ccoa/utils/reporting/summary.py
+++ b/src/ccoa/utils/reporting/summary.py
@@ -103,29 +103,39 @@ def compute_score_stats(scores: list[float]) -> dict[str, float | int | str]:
     return stats
 
 
-def log_summary(scores: list[float]) -> None:
-    """Log a one-line summary of the score distribution (count, percentiles, etc.)."""
-    if not scores:
+def log_summary(scores_by_column: dict[str, list[float]]) -> None:
+    """Log one INFO line per score column with count/percentiles/mean.
+
+    `scores_by_column` is keyed by output column name (e.g. `score___label__science`)
+    and preserves insertion order. An empty dict, or one whose every value list
+    is empty, logs a single warning instead.
+    """
+    if not scores_by_column or not any(scores_by_column.values()):
         logger.warning("No records classified.")
         return
 
-    stats = compute_score_stats(scores)
-    logger.info(
-        "score stats — count=%d min=%s p10=%s p25=%s p50=%s p75=%s p90=%s p95=%s p99=%s max=%s mean=%s median=%s stdev=%s",
-        stats["count"],
-        format_score(stats["min"]),
-        format_score(stats["p10"]),
-        format_score(stats["p25"]),
-        format_score(stats["p50"]),
-        format_score(stats["p75"]),
-        format_score(stats["p90"]),
-        format_score(stats["p95"]),
-        format_score(stats["p99"]),
-        format_score(stats["max"]),
-        format_score(stats["mean"]),
-        format_score(stats["median"]),
-        format_score(stats["stdev"]),
-    )
+    for column, scores in scores_by_column.items():
+        if not scores:
+            logger.warning("No scores collected for %s.", column)
+            continue
+        stats = compute_score_stats(scores)
+        logger.info(
+            "score stats [%s] — count=%d min=%s p10=%s p25=%s p50=%s p75=%s p90=%s p95=%s p99=%s max=%s mean=%s median=%s stdev=%s",
+            column,
+            stats["count"],
+            format_score(stats["min"]),
+            format_score(stats["p10"]),
+            format_score(stats["p25"]),
+            format_score(stats["p50"]),
+            format_score(stats["p75"]),
+            format_score(stats["p90"]),
+            format_score(stats["p95"]),
+            format_score(stats["p99"]),
+            format_score(stats["max"]),
+            format_score(stats["mean"]),
+            format_score(stats["median"]),
+            format_score(stats["stdev"]),
+        )
 
 
 def write_run_summary(
@@ -134,7 +144,7 @@ def write_run_summary(
     *,
     args: argparse.Namespace,
     resolved_count: int,
-    scores: list[float],
+    scores_by_column: dict[str, list[float]],
     processed: int,
     skipped_empty: int,
     skipped_homepage: int,
@@ -153,9 +163,10 @@ def write_run_summary(
     """Write a sidecar two-column CSV (`key,value`) capturing inputs + results.
 
     Sections (by key prefix): `run.*` (cli + timestamps), `arg.*` (every CLI flag),
-    `input.*` (resolved files), `count.*` (record counters), `score.*` (stats from
-    `compute_score_stats`), `time.*` (wall-clock / extract / predict). Designed to
-    be parsed back with `csv.reader` or `pandas.read_csv`.
+    `input.*` (resolved files), `count.*` (record counters),
+    `score.<column>.*` (per-column stats from `compute_score_stats`),
+    `time.*` (wall-clock / extract / predict). Designed to be parsed back with
+    `csv.reader` or `pandas.read_csv`.
     """
     rows: list[tuple[str, str]] = []
 
@@ -178,25 +189,31 @@ def write_run_summary(
     rows.append(("count.trafilatura_warnings", str(trafilatura_warnings)))
     rows.append(("count.trafilatura_errors", str(trafilatura_errors)))
 
-    stats = compute_score_stats(scores)
-    for key in (
-        "count",
-        "min",
-        "p10",
-        "p25",
-        "p50",
-        "p75",
-        "p90",
-        "p95",
-        "p99",
-        "max",
-        "mean",
-        "median",
-        "stdev",
-    ):
-        if key in stats:
-            value = stats[key]
-            rows.append((f"score.{key}", format_score(value) if key != "count" else str(value)))
+    for column, scores in scores_by_column.items():
+        stats = compute_score_stats(scores)
+        for key in (
+            "count",
+            "min",
+            "p10",
+            "p25",
+            "p50",
+            "p75",
+            "p90",
+            "p95",
+            "p99",
+            "max",
+            "mean",
+            "median",
+            "stdev",
+        ):
+            if key in stats:
+                value = stats[key]
+                rows.append(
+                    (
+                        f"score.{column}.{key}",
+                        format_score(value) if key != "count" else str(value),
+                    )
+                )
 
     rows.append(("time.total_seconds", f"{t_processing:.6f}"))
     rows.append(("time.extract_total_seconds", f"{t_extract_total:.6f}"))
diff --git a/tests/test_classify_warc.py b/tests/test_classify_warc.py
index 27f5226..e05acb8 100644
--- a/tests/test_classify_warc.py
+++ b/tests/test_classify_warc.py
@@ -9,7 +9,11 @@
 import pytest
 
 from ccoa.classifier.fasttext import FASTTEXT_MAX_INPUT_CHARS, clean_for_fasttext
-from ccoa.commands.classify_warc import ClassifyWarcCommand, load_resume_skipset
+from ccoa.commands.classify_warc import (
+    ClassifyWarcCommand,
+    _resolve_model_slots,
+    load_resume_skipset,
+)
 from ccoa.extraction.cache import (
     cache_path_for_warc,
     load_extraction_cache,
@@ -38,25 +42,32 @@ def _make_warcs(tmp_path, names):
 
 
 def test_summary_handles_empty(caplog):
-    """An empty score list logs a single warning and returns cleanly."""
+    """An empty mapping logs a single warning and returns cleanly."""
     caplog.set_level(logging.WARNING, logger="ccoa.utils.reporting.summary")
-    log_summary([])
+    log_summary({})
     warnings = [r for r in caplog.records if r.levelno == logging.WARNING]
     assert len(warnings) == 1
     assert "No records classified" in warnings[0].getMessage()
 
 
 def test_summary_computes_stats(caplog):
-    """A non-empty score list emits an INFO line containing count/min/max/mean."""
+    """Each column emits an INFO line containing count/min/max/mean for its scores."""
     caplog.set_level(logging.INFO, logger="ccoa.utils.reporting.summary")
-    log_summary([0.1, 0.5, 0.9])
+    log_summary(
+        {
+            "score___label__science": [0.1, 0.5, 0.9],
+            "score___label__cc": [0.9, 0.5, 0.1],
+        }
+    )
     messages = [r.getMessage() for r in caplog.records if r.levelno == logging.INFO]
-    assert messages, "expected an INFO summary line"
-    line = "\n".join(messages)
-    assert "count=3" in line
-    assert "min=0.100000" in line
-    assert "max=0.900000" in line
-    assert "mean=0.500000" in line
+    assert messages, "expected INFO summary lines"
+    text = "\n".join(messages)
+    assert "[score___label__science]" in text
+    assert "[score___label__cc]" in text
+    assert text.count("count=3") == 2
+    assert "min=0.100000" in text
+    assert "max=0.900000" in text
+    assert "mean=0.500000" in text
 
 
 def test_timing_zero_processed(caplog):
@@ -211,9 +222,9 @@ def _default_args(**overrides) -> argparse.Namespace:
         "files_limit": 0,
         "shuffle_files": False,
         "seed": 42,
-        "model_repo": "ibm-granite/GneissWeb.Sci_classifier",
-        "model_file": "fasttext_science.bin",
-        "target_label": "__label__science",
+        "model_repo": [],
+        "model_file": [],
+        "labels": [],
         "output": "-",
         "cache_dir": None,
         "anonymous_s3": False,
@@ -222,6 +233,7 @@ def _default_args(**overrides) -> argparse.Namespace:
         "workers_mode": "thread",
         "max_pool_restarts": 10,
         "resume_from_output": None,
+        "overwrite": False,
     }
     base.update(overrides)
     return argparse.Namespace(**base)
@@ -275,6 +287,36 @@ def test_run_aborts_when_summary_exists(tmp_path, caplog):
     assert any("already exists" in r.getMessage() for r in caplog.records)
 
 
+def test_run_overwrite_replaces_existing_output(monkeypatch, tmp_path, caplog):
+    """`--overwrite` skips the exists-guard and warns; the run proceeds and rewrites the file."""
+    warc_path = tmp_path / "empty.warc.gz"
+    _write_synthetic_warc(warc_path, "http://example.com/", b"<html></html>")
+    output_path = tmp_path / "out.csv"
+    output_path.write_text("STALE,DATA\n")  # pre-existing
+    summary_path = tmp_path / "out.summary.csv"
+    summary_path.write_text("stale,summary\n")
+
+    fake = _FakeModel(("__label__a",))
+    _patch_models(monkeypatch, {("r/x", "x.bin"): fake})
+
+    args = _default_args(
+        warc_paths=[str(warc_path)],
+        output=str(output_path),
+        model_repo=["r/x"],
+        model_file=["x.bin"],
+        overwrite=True,
+    )
+    caplog.set_level(logging.WARNING, logger="ccoa.commands.classify_warc")
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 0
+    warnings = [r.getMessage() for r in caplog.records if r.levelno == logging.WARNING]
+    assert any("Overwriting existing output" in m for m in warnings)
+    # The stale content must be replaced — fresh header on the first line.
+    assert output_path.read_text().splitlines()[0] == (
+        "URL,score___label__a,warc_filename,warc_record_index"
+    )
+
+
 def test_compute_score_stats_empty_and_singleton() -> None:
     """Empty -> {count: 0}; singleton -> n/a percentiles + stdev."""
     assert compute_score_stats([]) == {"count": 0}
@@ -286,7 +328,7 @@ def test_compute_score_stats_empty_and_singleton() -> None:
 
 
 def test_write_run_summary_round_trip(tmp_path):
-    """The sidecar CSV captures args, counters, score stats, and timings as key/value rows."""
+    """The sidecar CSV captures args, counters, per-column score stats, and timings."""
     summary_path = tmp_path / "out.summary.csv"
     args = _default_args(
         warc_paths=["s3://bucket/a.warc.gz", "s3://bucket/b.warc.gz"],
@@ -299,7 +341,10 @@ def test_write_run_summary_round_trip(tmp_path):
         {},
         args=args,
         resolved_count=2,
-        scores=[0.1, 0.5, 0.9],
+        scores_by_column={
+            "score___label__science": [0.1, 0.5, 0.9],
+            "score___label__cc": [0.9, 0.5, 0.1],
+        },
         processed=3,
         skipped_empty=1,
         skipped_homepage=2,
@@ -330,10 +375,12 @@ def test_write_run_summary_round_trip(tmp_path):
     assert rows["count.cache_hits"] == "4"
     assert rows["count.trafilatura_warnings"] == "17"
     assert rows["count.trafilatura_errors"] == "6"
-    assert rows["score.count"] == "3"
-    assert rows["score.min"] == "0.100000"
-    assert rows["score.max"] == "0.900000"
-    assert rows["score.mean"] == "0.500000"
+    assert rows["score.score___label__science.count"] == "3"
+    assert rows["score.score___label__science.min"] == "0.100000"
+    assert rows["score.score___label__science.max"] == "0.900000"
+    assert rows["score.score___label__science.mean"] == "0.500000"
+    assert rows["score.score___label__cc.count"] == "3"
+    assert rows["score.score___label__cc.mean"] == "0.500000"
     assert rows["time.total_seconds"] == "1.500000"
     assert rows["time.throughput_docs_per_sec"] == "2.000000"
     assert rows["run.started_at"] == "2026-05-24T00:00:00+00:00"
@@ -475,27 +522,52 @@ def test_clean_for_fasttext_clamps_length() -> None:
 
 
 def test_load_resume_skipset_groups_by_warc(tmp_path):
-    """A prior output CSV is grouped into {warc_filename: frozenset(record_indices)}."""
+    """A prior output CSV with the expected header is grouped by warc_filename."""
     csv_path = tmp_path / "prior.csv"
     csv_path.write_text(
-        "URL,prediction_score,warc_filename,warc_record_index\n"
+        "URL,score___label__science,warc_filename,warc_record_index\n"
         "http://a/,0.1,s3://bucket/a.warc.gz,0\n"
         "http://a/x,0.2,s3://bucket/a.warc.gz,3\n"
         "http://b/,0.3,s3://bucket/b.warc.gz,7\n"
     )
-    skipset = load_resume_skipset(str(csv_path), {})
+    expected = ["URL", "score___label__science", "warc_filename", "warc_record_index"]
+    skipset = load_resume_skipset(str(csv_path), {}, expected)
     assert skipset == {
         "s3://bucket/a.warc.gz": frozenset({0, 3}),
         "s3://bucket/b.warc.gz": frozenset({7}),
     }
 
 
-def test_load_resume_skipset_rejects_old_schema(tmp_path):
-    """CSV without warc_filename/warc_record_index raises a clear error."""
+def test_load_resume_skipset_rejects_pre_multi_label_csv(tmp_path):
+    """A pre-multi-label CSV (no `score_*` columns) is rejected with a clear message."""
     csv_path = tmp_path / "old.csv"
-    csv_path.write_text("URL,prediction_score\nhttp://a/,0.1\n")
-    with pytest.raises(ValueError, match="missing required columns"):
-        load_resume_skipset(str(csv_path), {})
+    csv_path.write_text("URL,prediction_score,warc_filename,warc_record_index\nhttp://a/,0.1,x,0\n")
+    expected = ["URL", "score___label__science", "warc_filename", "warc_record_index"]
+    with pytest.raises(ValueError, match="pre-multi-label"):
+        load_resume_skipset(str(csv_path), {}, expected)
+
+
+def test_load_resume_skipset_rejects_column_reorder(tmp_path):
+    """Same columns in a different order: rejected (concat-friendliness)."""
+    csv_path = tmp_path / "reordered.csv"
+    csv_path.write_text(
+        "URL,warc_filename,warc_record_index,score___label__science\nhttp://a/,x,0,0.1\n"
+    )
+    expected = ["URL", "score___label__science", "warc_filename", "warc_record_index"]
+    with pytest.raises(ValueError, match="header does not match"):
+        load_resume_skipset(str(csv_path), {}, expected)
+
+
+def test_load_resume_skipset_rejects_extra_columns(tmp_path):
+    """Prior CSV has columns the new run doesn't produce: rejected with named diff."""
+    csv_path = tmp_path / "extra.csv"
+    csv_path.write_text(
+        "URL,score___label__science,score___label__cc,warc_filename,warc_record_index\n"
+        "http://a/,0.1,0.9,x,0\n"
+    )
+    expected = ["URL", "score___label__science", "warc_filename", "warc_record_index"]
+    with pytest.raises(ValueError, match="score___label__cc"):
+        load_resume_skipset(str(csv_path), {}, expected)
 
 
 def _write_synthetic_warc(warc_path, url: str, body: bytes) -> None:
@@ -549,16 +621,21 @@ def test_classify_warc_end_to_end(tmp_path):
     args = _default_args(
         warc_paths=[str(warc_path)],
         output=str(output_path),
-        model_repo="facebook/fasttext-language-identification",
-        model_file="model.bin",
-        target_label="__label__eng_Latn",
+        model_repo=["facebook/fasttext-language-identification"],
+        model_file=["model.bin"],
+        labels=["__label__eng_Latn"],
     )
     rc = ClassifyWarcCommand().run(args)
     assert rc == 0
 
     with output_path.open(newline="", encoding="utf-8") as fh:
         rows = list(csv.reader(fh))
-    assert rows[0] == ["URL", "prediction_score", "warc_filename", "warc_record_index"]
+    assert rows[0] == [
+        "URL",
+        "score___label__eng_Latn",
+        "warc_filename",
+        "warc_record_index",
+    ]
     assert len(rows) == 2, f"expected 1 data row, got {len(rows) - 1}"
 
     url, score, warc_filename, record_index = rows[1]
@@ -573,3 +650,239 @@ def test_classify_warc_end_to_end(tmp_path):
         summary_rows = dict(list(csv.reader(fh))[1:])
     assert summary_rows["count.processed"] == "1"
     assert summary_rows["input.resolved_count"] == "1"
+    assert summary_rows["score.score___label__eng_Latn.count"] == "1"
+
+
+def test_resolve_model_slots_defaults_when_both_empty():
+    """No CLI input → falls back to the default science model with '*' labels."""
+    slots = _resolve_model_slots([], [], [])
+    assert slots == [
+        ("ibm-granite/GneissWeb.Sci_classifier", "fasttext_science.bin", ["*"]),
+    ]
+
+
+def test_resolve_model_slots_zips_parallel_lists():
+    """Parallel `--model-repo`/`--model-file`/`--labels` lists are zipped positionally."""
+    slots = _resolve_model_slots(
+        ["r1/m", "r2/m"],
+        ["f1.bin", "f2.bin"],
+        ["__label__a,__label__b", "*"],
+    )
+    assert slots == [
+        ("r1/m", "f1.bin", ["__label__a", "__label__b"]),
+        ("r2/m", "f2.bin", ["*"]),
+    ]
+
+
+def test_resolve_model_slots_broadcasts_missing_labels():
+    """Omitted `--labels` defaults to '*' for every model slot."""
+    slots = _resolve_model_slots(["r1/m", "r2/m"], ["f1.bin", "f2.bin"], [])
+    assert slots == [
+        ("r1/m", "f1.bin", ["*"]),
+        ("r2/m", "f2.bin", ["*"]),
+    ]
+
+
+def test_resolve_model_slots_rejects_length_mismatch():
+    """`--model-repo` and `--model-file` of different lengths fails."""
+    with pytest.raises(ValueError, match="same length"):
+        _resolve_model_slots(["r1", "r2"], ["f1.bin"], [])
+
+
+def test_resolve_model_slots_rejects_label_length_mismatch():
+    """A `--labels` list whose length doesn't match the number of models fails."""
+    with pytest.raises(ValueError, match="one entry per --model-repo"):
+        _resolve_model_slots(["r1", "r2"], ["f1.bin", "f2.bin"], ["only_one"])
+
+
+class _FakeModel:
+    """Stand-in for a loaded fasttext model.
+
+    `get_labels` returns the labels we want exposed to '*' expansion; `predict`
+    deterministically returns a 1.0 for the first label and 0.0 for the rest,
+    so the test can assert on output values without depending on the real model.
+    """
+
+    def __init__(self, labels: tuple[str, ...]):
+        self._labels = labels
+
+    def get_labels(self) -> tuple[str, ...]:
+        return self._labels
+
+    def predict(self, _text: str, k: int = -1):
+        import numpy as np
+
+        probs = np.zeros(len(self._labels), dtype=float)
+        if probs.size:
+            probs[0] = 1.0
+        return (tuple(self._labels), probs)
+
+
+def _patch_models(monkeypatch, model_for: dict[tuple[str, str], _FakeModel]) -> None:
+    """Patch `load_classifier`/`get_model_labels` to return `_FakeModel`s by (repo, file)."""
+
+    def fake_load(repo: str, file: str):
+        try:
+            return model_for[(repo, file)]
+        except KeyError as exc:
+            raise AssertionError(f"unexpected load_classifier({repo!r}, {file!r})") from exc
+
+    def fake_get_labels(model):
+        return tuple(model.get_labels())
+
+    monkeypatch.setattr("ccoa.commands.classify_warc.load_classifier", fake_load)
+    monkeypatch.setattr("ccoa.commands.classify_warc.get_model_labels", fake_get_labels)
+
+
+def test_run_expands_star_labels_from_model(monkeypatch, tmp_path):
+    """'*' expands to the model's full label set via `get_model_labels`."""
+    warc_path = tmp_path / "empty.warc.gz"
+    _write_synthetic_warc(warc_path, "http://example.com/", b"<html></html>")
+    fake = _FakeModel(("__label__alpha", "__label__beta"))
+    _patch_models(monkeypatch, {("r/x", "x.bin"): fake})
+
+    output_path = tmp_path / "out.csv"
+    args = _default_args(
+        warc_paths=[str(warc_path)],
+        output=str(output_path),
+        model_repo=["r/x"],
+        model_file=["x.bin"],
+    )
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 0
+
+    with output_path.open(newline="", encoding="utf-8") as fh:
+        rows = list(csv.reader(fh))
+    assert rows[0] == [
+        "URL",
+        "score___label__alpha",
+        "score___label__beta",
+        "warc_filename",
+        "warc_record_index",
+    ]
+
+
+def test_run_writes_multi_model_header(monkeypatch, tmp_path):
+    """Two models: columns get a `m<idx>_` prefix so per-model labels never collide."""
+    warc_path = tmp_path / "empty.warc.gz"
+    _write_synthetic_warc(warc_path, "http://example.com/", b"<html></html>")
+    fake_a = _FakeModel(("__label__sci", "__label__cc"))
+    fake_b = _FakeModel(("__label__hq", "__label__lq"))
+    _patch_models(
+        monkeypatch,
+        {
+            ("rA/sci", "sci.bin"): fake_a,
+            ("rB/q", "q.bin"): fake_b,
+        },
+    )
+
+    output_path = tmp_path / "out.csv"
+    args = _default_args(
+        warc_paths=[str(warc_path)],
+        output=str(output_path),
+        model_repo=["rA/sci", "rB/q"],
+        model_file=["sci.bin", "q.bin"],
+    )
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 0
+
+    with output_path.open(newline="", encoding="utf-8") as fh:
+        rows = list(csv.reader(fh))
+    assert rows[0] == [
+        "URL",
+        "score_m0___label__sci",
+        "score_m0___label__cc",
+        "score_m1___label__hq",
+        "score_m1___label__lq",
+        "warc_filename",
+        "warc_record_index",
+    ]
+
+
+def test_run_disambiguates_shared_label_across_models(monkeypatch, tmp_path):
+    """Two models sharing a label name (e.g. `__label__cc`) coexist via the `m<idx>_` prefix."""
+    warc_path = tmp_path / "empty.warc.gz"
+    _write_synthetic_warc(warc_path, "http://example.com/", b"<html></html>")
+    fake_a = _FakeModel(("__label__shared",))
+    fake_b = _FakeModel(("__label__shared",))
+    _patch_models(
+        monkeypatch,
+        {
+            ("rA/m", "a.bin"): fake_a,
+            ("rB/m", "b.bin"): fake_b,
+        },
+    )
+
+    output_path = tmp_path / "out.csv"
+    args = _default_args(
+        warc_paths=[str(warc_path)],
+        output=str(output_path),
+        model_repo=["rA/m", "rB/m"],
+        model_file=["a.bin", "b.bin"],
+    )
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 0
+
+    with output_path.open(newline="", encoding="utf-8") as fh:
+        rows = list(csv.reader(fh))
+    assert rows[0] == [
+        "URL",
+        "score_m0___label__shared",
+        "score_m1___label__shared",
+        "warc_filename",
+        "warc_record_index",
+    ]
+
+
+def test_run_rejects_mismatched_model_lengths(caplog):
+    """--model-repo and --model-file of different lengths: rc=2, no work done."""
+    caplog.set_level(logging.ERROR, logger="ccoa.commands.classify_warc")
+    args = _default_args(
+        model_repo=["r1/m", "r2/m"],
+        model_file=["only_one.bin"],
+        output="-",
+    )
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 2
+
+
+def test_run_scores_all_resolved_labels(monkeypatch, tmp_path):
+    """Each record is scored against every requested label; outputs land in column order."""
+    warc_path = tmp_path / "doc.warc.gz"
+    html = (
+        b"<!DOCTYPE html><html><head><title>T</title></head><body>"
+        b"<p>Long enough body to survive trafilatura's minimum-content gates. "
+        b"We just need a single response record to flow through the pipeline.</p>"
+        b"</body></html>"
+    )
+    _write_synthetic_warc(warc_path, "http://example.com/doc", html)
+    fake = _FakeModel(("__label__alpha", "__label__beta"))
+    _patch_models(monkeypatch, {("r/x", "x.bin"): fake})
+
+    output_path = tmp_path / "out.csv"
+    args = _default_args(
+        warc_paths=[str(warc_path)],
+        output=str(output_path),
+        model_repo=["r/x"],
+        model_file=["x.bin"],
+    )
+    rc = ClassifyWarcCommand().run(args)
+    assert rc == 0
+
+    with output_path.open(newline="", encoding="utf-8") as fh:
+        rows = list(csv.reader(fh))
+    assert rows[0] == [
+        "URL",
+        "score___label__alpha",
+        "score___label__beta",
+        "warc_filename",
+        "warc_record_index",
+    ]
+    assert len(rows) == 2
+    url, alpha, beta, warc_filename, record_index = rows[1]
+    assert url == "http://example.com/doc"
+    # `_FakeModel.predict` is rigged: 1.0 for the first label, 0.0 for the rest.
+    assert float(alpha) == 1.0
+    assert float(beta) == 0.0
+    assert warc_filename == str(warc_path)
+    assert record_index == "0"