diff --git a/.github/workflows/monthly-release.yml b/.github/workflows/monthly-release.yml
index ef40f5e..fa3876d 100644
--- a/.github/workflows/monthly-release.yml
+++ b/.github/workflows/monthly-release.yml
@@ -59,11 +59,11 @@ jobs:
         if: steps.prep.outputs.release_needed == 'true'
         uses: astral-sh/setup-uv@v5
 
-      - name: Set up Python 3.13
+      - name: Set up Python 3.12
         if: steps.prep.outputs.release_needed == 'true'
         uses: actions/setup-python@v5
         with:
-          python-version: "3.13"
+          python-version: "3.12"
 
       - name: Run release checks
         if: steps.prep.outputs.release_needed == 'true'
@@ -78,7 +78,7 @@ jobs:
         run: |
           uv build
           uvx twine check dist/*
-          uv venv /tmp/uxarray-mcp-wheel --python 3.13
+          uv venv /tmp/uxarray-mcp-wheel --python 3.12
           uv pip install --python /tmp/uxarray-mcp-wheel/bin/python dist/uxarray_mcp-*.whl
           /tmp/uxarray-mcp-wheel/bin/uxarray-mcp --help
 
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e607f6d..e0b53e8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,20 @@ uses Semantic Versioning for public releases.
 
 ## Unreleased
 
+### Breaking
+- `gradient` and `curl` now default `scale_by_radius=True`, matching UXarray's
+  public API. Pass `scale_by_radius=False` explicitly to preserve the previous
+  MCP unit-sphere behavior. This default alignment requires a minor release.
+
 ### Added
+- Vector-calculus results now include a machine-actionable
+  `scientific_status` with `status`, `physically_interpretable`, stable warning
+  codes, and warning text.
+- `get_capabilities` now returns variable units, standard names, dimensions,
+  and a `scientific_contracts.vector_calculus` block that separates structural
+  applicability from semantic suitability.
+- Persistent JSON and NetCDF artifacts are written atomically under a process
+  lock so concurrent tool calls cannot expose partial result files.
 - `gradient`, `curl`, and `divergence` (via `run_analysis` and the
   `calculate_gradient`/`calculate_curl`/`calculate_divergence` tools) now
   accept `time_index`/`level_index` to select a single time/level slice
@@ -27,6 +40,14 @@ uses Semantic Versioning for public releases.
   reading the tool's structured JSON result.
 
 ### Fixed
+- `analyze_dataset` and `run_scientific_agent` now actually skip zonal
+  statistics after failed validation; `analyze_dataset` also skips variable
+  plotting and records the validation summary in provenance.
+- `analyze_dataset` derives aggregate execution venue from completed stage
+  provenance instead of labeling a fallback-local run as HPC solely because
+  `use_remote=True` was requested.
+- The monthly release workflow now uses supported Python 3.12 instead of Python
+  3.13, which conflicts with the package's `requires-python` constraint.
 - `calculate_area` (local and remote) silently defaulted `area_units` to
   `"m^2"` whenever a grid's `face_areas` carried no `units` attribute at
   all, fabricating a label the source file never provided. It now reports
diff --git a/README.md b/README.md
index 352c2ec..55621a2 100644
--- a/README.md
+++ b/README.md
@@ -191,10 +191,16 @@ auditable and the server actively flags common scientific pitfalls:
 - **Derivative unit convention is never hidden.** `gradient`, `curl`, and
   `divergence` echo `scale_by_radius` in both the result and provenance, so a
   unit-sphere result can never be mistaken for a physical (per-metre) one.
+  Gradient and curl default to physical scaling, matching UXarray; pass
+  `scale_by_radius=False` explicitly for unit-sphere output.
 - **Vector-calculus sanity guard.** `curl`/`divergence` warn (without blocking)
   when the two inputs are the same field, or when neither carries a
   velocity/flux-like `units` attribute — the classic "vorticity from two random
-  scalars" mistake now surfaces a warning in `_provenance.warnings`.
+  scalars" mistake now surfaces a warning in `_provenance.warnings` and a
+  machine-actionable `scientific_status` with stable warning codes.
+- **Applicability is not suitability.** `get_capabilities` reports whether
+  vector operations are structurally computable separately from whether
+  metadata supports physical interpretation.
 - **Local/remote version drift is surfaced.** Remote results record the
   worker's *actual* UXarray version (`remote_uxarray_version`) and emit a
   warning when it differs from the local version, so silent numerical
diff --git a/docs/tools.md b/docs/tools.md
index 6c6846a..2fd7f5c 100644
--- a/docs/tools.md
+++ b/docs/tools.md
@@ -73,14 +73,20 @@ missing.
 
 `gradient`, `curl`, and `divergence` echo the `scale_by_radius` convention in
 their result and provenance. `gradient` and `curl` accept `scale_by_radius`
-(default `False`). When `False`, results stay on the unit sphere (the historical
-behavior). Set it to `True` to divide by `uxgrid.sphere_radius` for physical
-units; the grid must define `sphere_radius`.
+(default `True`, matching UXarray). When `True`, results are divided by
+`uxgrid.sphere_radius` for physical units; the grid must define
+`sphere_radius`. Pass `False` explicitly to keep unit-sphere results.
 
 `curl` and `divergence` also emit **vector-component warnings**: if the two
 inputs are the same field, or neither carries a velocity/flux-like `units`
 attribute, a warning is added to `_provenance.warnings`. The computation still
-runs (the math is valid), but the result is flagged as possibly non-physical.
+runs (the math is valid), but `scientific_status` marks the result as a warning,
+sets `physically_interpretable` to false, and includes stable warning codes.
+
+`get_capabilities` distinguishes **structural applicability** from **semantic
+suitability** for vector operations. Two face-centered arrays make curl and
+divergence computable, but physical interpretation remains `unverified` unless
+metadata supplies vector-like units or standard names.
 
 `zonal_anomaly` and `remap_to_rectilinear` are backed by
 `UxDataArray.zonal_anomaly` and `UxDataArray.remap.to_rectilinear`, available in
diff --git a/src/uxarray_mcp/domain/vector_calc.py b/src/uxarray_mcp/domain/vector_calc.py
index 04526fb..76cb7ce 100644
--- a/src/uxarray_mcp/domain/vector_calc.py
+++ b/src/uxarray_mcp/domain/vector_calc.py
@@ -102,7 +102,7 @@ def _vector_component_warnings(
     u: Any,
     v: Any,
     operation: str,
-) -> list[str]:
+) -> tuple[list[str], list[str]]:
     """Return soft warnings when (u, v) do not look like real vector components.
 
     curl and divergence are only physically meaningful when ``u`` and ``v`` are
@@ -118,8 +118,10 @@ def _vector_component_warnings(
     scientist — can see that the inputs were suspicious.
     """
     warnings: list[str] = []
+    warning_codes: list[str] = []
 
     if u_variable == v_variable:
+        warning_codes.append("VECTOR_COMPONENTS_IDENTICAL")
         warnings.append(
             f"{operation}: u_variable and v_variable are the same field "
             f"('{u_variable}'). {operation} is only physically meaningful for a "
@@ -138,6 +140,7 @@ def _looks_velocity(unit: str) -> bool:
         return any(hint in unit for hint in _VELOCITY_LIKE_UNIT_HINTS)
 
     if not (_looks_velocity(u_unit) or _looks_velocity(v_unit)):
+        warning_codes.append("VECTOR_UNITS_UNVERIFIED")
         seen = ", ".join(
             f"{name}='{unit or 'unset'}'"
             for name, unit in ((u_variable, u_unit), (v_variable, v_unit))
@@ -149,13 +152,13 @@ def _looks_velocity(unit: str) -> bool:
             f"{operation} physically."
         )
 
-    return warnings
+    return warnings, warning_codes
 
 
 def compute_gradient(
     uxds: Any,
     variable_name: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
 ) -> dict:
@@ -171,10 +174,10 @@ def compute_gradient(
         Loaded UXarray dataset.
     variable_name : str
         Face-centered scalar variable to differentiate.
-    scale_by_radius : bool, default False
+    scale_by_radius : bool, default True
         When ``True``, divide the unit-sphere derivatives by
         ``uxgrid.sphere_radius`` to return physical units (requires a grid with
-        ``sphere_radius``). The default ``False`` keeps the unit-sphere result.
+        ``sphere_radius``). Pass ``False`` to keep the unit-sphere result.
     time_index : int, default 0
         Time index to select if the variable carries a leading time
         dimension. Ignored if there is no time dimension.
@@ -229,7 +232,7 @@ def _stats(arr: Any) -> dict:
 
     components = {name: _stats(grad[name]) for name in comp_names}
 
-    return {
+    result = {
         "variable_name": variable_name,
         "components": comp_names,
         "component_stats": components,
@@ -238,13 +241,20 @@ def _stats(arr: Any) -> dict:
         "interpretation": "zonal (∂/∂x) and meridional (∂/∂y) components of the gradient",
         "component_warnings": uxarray_warnings,
     }
+    from uxarray_mcp.provenance import attach_scientific_status
+
+    return attach_scientific_status(
+        result,
+        warnings=uxarray_warnings,
+        warning_codes=["SPHERE_RADIUS_UNAVAILABLE"] if uxarray_warnings else [],
+    )
 
 
 def compute_curl(
     uxds: Any,
     u_variable: str,
     v_variable: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
 ) -> dict:
@@ -263,10 +273,10 @@ def compute_curl(
         Zonal (east–west) component variable name.
     v_variable : str
         Meridional (north–south) component variable name.
-    scale_by_radius : bool, default False
+    scale_by_radius : bool, default True
         When ``True``, divide the unit-sphere result by ``uxgrid.sphere_radius``
         to return physical units (``1/s`` for wind in ``m/s``; requires a grid
-        with ``sphere_radius``). The default ``False`` keeps the unit sphere.
+        with ``sphere_radius``). Pass ``False`` to keep the unit sphere.
     time_index : int, default 0
         Time index to select if the components carry a leading time
         dimension. Ignored if there is no time dimension.
@@ -298,7 +308,7 @@ def compute_curl(
 
     import numpy as np
 
-    component_warnings = _vector_component_warnings(
+    component_warnings, warning_codes = _vector_component_warnings(
         u_variable, v_variable, u, v, "curl"
     )
 
@@ -316,6 +326,8 @@ def compute_curl(
             seen.add(w)
             deduped_uxarray_warnings.append(w)
     component_warnings = component_warnings + deduped_uxarray_warnings
+    if deduped_uxarray_warnings:
+        warning_codes.append("SPHERE_RADIUS_UNAVAILABLE")
     vals = result.values
     finite = vals[np.isfinite(vals)]
 
@@ -330,7 +342,7 @@ def compute_curl(
     else:
         stats = {"min": None, "max": None, "mean": None, "std": None}
 
-    return {
+    result = {
         "u_variable": u_variable,
         "v_variable": v_variable,
         "interpretation": "relative vorticity ζ = ∂v/∂x − ∂u/∂y",
@@ -339,6 +351,11 @@ def compute_curl(
         "stats": stats,
         "component_warnings": component_warnings,
     }
+    from uxarray_mcp.provenance import attach_scientific_status
+
+    return attach_scientific_status(
+        result, warnings=component_warnings, warning_codes=warning_codes
+    )
 
 
 def compute_divergence(
@@ -395,11 +412,16 @@ def compute_divergence(
 
     import numpy as np
 
-    component_warnings = _vector_component_warnings(
+    component_warnings, warning_codes = _vector_component_warnings(
         u_variable, v_variable, u, v, "divergence"
     )
 
-    result = u.divergence(v)
+    result, uxarray_warnings = _call_capturing_warnings(lambda: u.divergence(v))
+    for warning in uxarray_warnings:
+        if warning not in component_warnings:
+            component_warnings.append(warning)
+    if uxarray_warnings:
+        warning_codes.append("SPHERE_RADIUS_UNAVAILABLE")
     vals = result.values
     finite = vals[np.isfinite(vals)]
 
@@ -414,7 +436,7 @@ def compute_divergence(
     else:
         stats = {"min": None, "max": None, "mean": None, "std": None}
 
-    return {
+    output = {
         "u_variable": u_variable,
         "v_variable": v_variable,
         "interpretation": "horizontal divergence ∂u/∂x + ∂v/∂y",
@@ -422,6 +444,11 @@ def compute_divergence(
         "stats": stats,
         "component_warnings": component_warnings,
     }
+    from uxarray_mcp.provenance import attach_scientific_status
+
+    return attach_scientific_status(
+        output, warnings=component_warnings, warning_codes=warning_codes
+    )
 
 
 def compute_azimuthal_mean(
diff --git a/src/uxarray_mcp/provenance.py b/src/uxarray_mcp/provenance.py
index 94951b9..b0d40e1 100644
--- a/src/uxarray_mcp/provenance.py
+++ b/src/uxarray_mcp/provenance.py
@@ -77,3 +77,20 @@ def attach_provenance(
         provenance["validation_summary"] = validation_summary
     result["_provenance"] = provenance
     return result
+
+
+def attach_scientific_status(
+    result: dict[str, Any],
+    *,
+    warnings: list[str] | None = None,
+    warning_codes: list[str] | None = None,
+) -> dict[str, Any]:
+    """Attach a machine-actionable scientific interpretation status."""
+    messages = warnings or []
+    result["scientific_status"] = {
+        "status": "warning" if messages else "complete",
+        "physically_interpretable": not messages,
+        "warning_codes": warning_codes or [],
+        "warnings": messages,
+    }
+    return result
diff --git a/src/uxarray_mcp/remote/agent.py b/src/uxarray_mcp/remote/agent.py
index 31e1746..83ff330 100644
--- a/src/uxarray_mcp/remote/agent.py
+++ b/src/uxarray_mcp/remote/agent.py
@@ -218,7 +218,7 @@ async def calculate_gradient_remote(
         grid_path: str,
         data_path: str,
         variable_name: str,
-        scale_by_radius: bool = False,
+        scale_by_radius: bool = True,
         time_index: int = 0,
         level_index: int = 0,
     ) -> Dict[str, Any]:
@@ -240,7 +240,7 @@ async def calculate_curl_remote(
         data_path: str,
         u_variable: str,
         v_variable: str,
-        scale_by_radius: bool = False,
+        scale_by_radius: bool = True,
         time_index: int = 0,
         level_index: int = 0,
     ) -> Dict[str, Any]:
diff --git a/src/uxarray_mcp/remote/compute_functions.py b/src/uxarray_mcp/remote/compute_functions.py
index e91275e..8b59422 100644
--- a/src/uxarray_mcp/remote/compute_functions.py
+++ b/src/uxarray_mcp/remote/compute_functions.py
@@ -1051,7 +1051,7 @@ def remote_calculate_gradient(
     grid_path: str,
     data_path: str,
     variable_name: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
 ) -> Dict[str, Any]:
@@ -1146,6 +1146,14 @@ def _stats(arr: Any) -> Dict[str, Any]:
         "scale_by_radius": applied_scale,
         "interpretation": "zonal (d/dx) and meridional (d/dy) components of the gradient",
         "component_warnings": uxarray_warnings,
+        "scientific_status": {
+            "status": "warning" if uxarray_warnings else "complete",
+            "physically_interpretable": not uxarray_warnings,
+            "warning_codes": (
+                ["SPHERE_RADIUS_UNAVAILABLE"] if uxarray_warnings else []
+            ),
+            "warnings": uxarray_warnings,
+        },
         "_worker_uxarray_version": getattr(ux, "__version__", "unknown"),
     }
 
@@ -1155,7 +1163,7 @@ def remote_calculate_curl(
     data_path: str,
     u_variable: str,
     v_variable: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
 ) -> Dict[str, Any]:
@@ -1220,7 +1228,9 @@ def remote_calculate_curl(
     # Soft semantic guardrail: curl is only physical for genuine vector
     # components. Warn (do not block) on same-field or non-velocity inputs.
     component_warnings = []
+    warning_codes = []
     if u_variable == v_variable:
+        warning_codes.append("VECTOR_COMPONENTS_IDENTICAL")
         component_warnings.append(
             f"curl: u_variable and v_variable are the same field "
             f"('{u_variable}'); result is a mathematical artifact, not a "
@@ -1230,6 +1240,7 @@ def remote_calculate_curl(
     _v_units = str((getattr(v, "attrs", {}) or {}).get("units", "")).strip().lower()
     _vel_hints = ("m/s", "m s-1", "m s^-1", "cm/s", "km/h", "pa/s", "kg/m2/s")
     if not any(h in _u_units or h in _v_units for h in _vel_hints):
+        warning_codes.append("VECTOR_UNITS_UNVERIFIED")
         component_warnings.append(
             f"curl: neither component has a velocity-like 'units' attribute "
             f"(u='{_u_units or 'unset'}', v='{_v_units or 'unset'}'); verify "
@@ -1255,6 +1266,7 @@ def remote_calculate_curl(
             if _msg not in _seen:
                 _seen.add(_msg)
                 component_warnings.append(_msg)
+                warning_codes.append("SPHERE_RADIUS_UNAVAILABLE")
     vals = result.values
     finite = vals[np.isfinite(vals)]
     stats: Dict[str, Any] = (
@@ -1275,6 +1287,12 @@ def remote_calculate_curl(
         "scale_by_radius": applied_scale,
         "stats": stats,
         "component_warnings": component_warnings,
+        "scientific_status": {
+            "status": "warning" if component_warnings else "complete",
+            "physically_interpretable": not component_warnings,
+            "warning_codes": warning_codes,
+            "warnings": component_warnings,
+        },
         "_worker_uxarray_version": getattr(ux, "__version__", "unknown"),
     }
 
@@ -1347,7 +1365,9 @@ def remote_calculate_divergence(
     # Soft semantic guardrail (mirrors curl): warn on same-field or
     # non-velocity inputs; do not block.
     component_warnings = []
+    warning_codes = []
     if u_variable == v_variable:
+        warning_codes.append("VECTOR_COMPONENTS_IDENTICAL")
         component_warnings.append(
             f"divergence: u_variable and v_variable are the same field "
             f"('{u_variable}'); result is a mathematical artifact, not a "
@@ -1357,13 +1377,23 @@ def remote_calculate_divergence(
     _v_units = str((getattr(v, "attrs", {}) or {}).get("units", "")).strip().lower()
     _vel_hints = ("m/s", "m s-1", "m s^-1", "cm/s", "km/h", "pa/s", "kg/m2/s")
     if not any(h in _u_units or h in _v_units for h in _vel_hints):
+        warning_codes.append("VECTOR_UNITS_UNVERIFIED")
         component_warnings.append(
             f"divergence: neither component has a velocity-like 'units' "
             f"attribute (u='{_u_units or 'unset'}', v='{_v_units or 'unset'}'); "
             "verify inputs are genuine vector components before interpreting."
         )
 
-    result = u.divergence(v)
+    import warnings as _warnings_module
+
+    with _warnings_module.catch_warnings(record=True) as _caught:
+        _warnings_module.simplefilter("always")
+        result = u.divergence(v)
+        for _warning in _caught:
+            _message = str(_warning.message)
+            if _message not in component_warnings:
+                component_warnings.append(_message)
+                warning_codes.append("SPHERE_RADIUS_UNAVAILABLE")
     vals = result.values
     finite = vals[np.isfinite(vals)]
     stats: Dict[str, Any] = (
@@ -1383,6 +1413,12 @@ def remote_calculate_divergence(
         "n_face": int(uxds.uxgrid.n_face),
         "stats": stats,
         "component_warnings": component_warnings,
+        "scientific_status": {
+            "status": "warning" if component_warnings else "complete",
+            "physically_interpretable": not component_warnings,
+            "warning_codes": warning_codes,
+            "warnings": component_warnings,
+        },
         "_worker_uxarray_version": getattr(ux, "__version__", "unknown"),
     }
 
@@ -1492,7 +1528,8 @@ def remote_grid_facts(
 
         variables = []
         for var_name in uxds.data_vars:
-            dims = uxds[var_name].dims
+            variable = uxds[var_name]
+            dims = variable.dims
             if any(d in dims for d in ("n_face", "nCells")):
                 location = "faces"
             elif any(d in dims for d in ("n_node", "nVertices")):
@@ -1501,7 +1538,16 @@ def remote_grid_facts(
                 location = "edges"
             else:
                 location = "other"
-            variables.append({"name": str(var_name), "location": location})
+            variables.append(
+                {
+                    "name": str(var_name),
+                    "location": location,
+                    "dims": list(dims),
+                    "units": variable.attrs.get("units"),
+                    "standard_name": variable.attrs.get("standard_name"),
+                    "long_name": variable.attrs.get("long_name"),
+                }
+            )
         facts["variables"] = variables
 
     return facts
diff --git a/src/uxarray_mcp/state.py b/src/uxarray_mcp/state.py
index bc3a18d..eb9bd2f 100644
--- a/src/uxarray_mcp/state.py
+++ b/src/uxarray_mcp/state.py
@@ -5,6 +5,8 @@
 import json
 import os
 import shutil
+import tempfile
+import threading
 import uuid
 from contextlib import suppress
 from dataclasses import dataclass, field
@@ -14,6 +16,8 @@
 
 import xarray as xr
 
+_WRITE_LOCK = threading.RLock()
+
 
 def _now_utc() -> str:
     return datetime.now(timezone.utc).isoformat()
@@ -51,7 +55,9 @@ def _json_safe(value: Any) -> Any:
 
 
 def _write_json(path: Path, payload: dict[str, Any]) -> None:
-    path.write_text(json.dumps(_json_safe(payload), indent=2, sort_keys=True))
+    serialized = json.dumps(_json_safe(payload), indent=2, sort_keys=True)
+    with _WRITE_LOCK:
+        _atomic_write(path, lambda temporary: temporary.write_text(serialized))
 
 
 def _read_json(path: Path) -> dict[str, Any]:
@@ -74,6 +80,21 @@ def _result_path(result_id: str, suffix: str) -> Path:
     return _artifacts_dir() / f"{result_id}{suffix}"
 
 
+def _atomic_write(path: Path, writer: Any) -> None:
+    """Write through a sibling temporary file and atomically replace ``path``."""
+    _ensure_dir(path.parent)
+    fd, temporary_name = tempfile.mkstemp(
+        prefix=f".{path.name}.", suffix=".tmp", dir=path.parent
+    )
+    os.close(fd)
+    temporary = Path(temporary_name)
+    try:
+        writer(temporary)
+        os.replace(temporary, path)
+    finally:
+        temporary.unlink(missing_ok=True)
+
+
 def _sanitize_netcdf_attr_value(value: Any) -> Any:
     if isinstance(value, bool):
         return int(value)
@@ -318,13 +339,17 @@ def save_result(result: dict[str, Any]) -> dict[str, Any]:
 
 def write_grid_artifact(grid: Any, result_id: str) -> str:
     path = _result_path(result_id, ".nc")
-    grid.to_xarray().to_netcdf(path)
+    with _WRITE_LOCK:
+        _atomic_write(path, lambda temporary: grid.to_xarray().to_netcdf(temporary))
     return str(path)
 
 
 def write_dataarray_artifact(data: Any, result_id: str) -> str:
     path = _result_path(result_id, ".nc")
-    _sanitize_netcdf_attrs(data).to_netcdf(path)
+    with _WRITE_LOCK:
+        _atomic_write(
+            path, lambda temporary: _sanitize_netcdf_attrs(data).to_netcdf(temporary)
+        )
     return str(path)
 
 
@@ -342,7 +367,8 @@ def write_dataset_artifact(data: Any, result_id: str) -> str:
             str(key): _sanitize_netcdf_attr_value(value)
             for key, value in sanitized[name].attrs.items()
         }
-    sanitized.to_netcdf(path)
+    with _WRITE_LOCK:
+        _atomic_write(path, lambda temporary: sanitized.to_netcdf(temporary))
     return str(path)
 
 
diff --git a/src/uxarray_mcp/tools/capabilities.py b/src/uxarray_mcp/tools/capabilities.py
index a0feed8..ece6fd2 100644
--- a/src/uxarray_mcp/tools/capabilities.py
+++ b/src/uxarray_mcp/tools/capabilities.py
@@ -9,6 +9,29 @@
 from uxarray_mcp.provenance import attach_provenance
 from uxarray_mcp.remote.config import load_config
 
+_VECTOR_UNIT_HINTS = (
+    "m/s",
+    "m s-1",
+    "m s^-1",
+    "meter/second",
+    "meters/second",
+    "cm/s",
+    "km/h",
+    "kg/m2/s",
+    "kg m-2 s-1",
+    "pa/s",
+    "n/m2",
+)
+
+
+def _looks_vector_like(variable: Dict[str, Any]) -> bool:
+    units = str(variable.get("units") or "").strip().lower()
+    standard_name = str(variable.get("standard_name") or "").strip().lower()
+    return any(hint in units for hint in _VECTOR_UNIT_HINTS) or any(
+        hint in standard_name
+        for hint in ("eastward_", "northward_", "grid_eastward_", "grid_northward_")
+    )
+
 
 def _uxarray_supports(attr_path: str) -> bool:
     """Return True if ``UxDataArray`` provides the given (possibly dotted) attr.
@@ -77,7 +100,8 @@ def _local_grid_facts(grid_path: str, data_path: Optional[str]) -> Dict[str, Any
 
         variables = []
         for var_name in uxds.data_vars:
-            dims = uxds[var_name].dims
+            variable = uxds[var_name]
+            dims = variable.dims
             if any(d in dims for d in ("n_face", "nCells")):
                 location = "faces"
             elif any(d in dims for d in ("n_node", "nVertices")):
@@ -86,7 +110,16 @@ def _local_grid_facts(grid_path: str, data_path: Optional[str]) -> Dict[str, Any
                 location = "edges"
             else:
                 location = "other"
-            variables.append({"name": str(var_name), "location": location})
+            variables.append(
+                {
+                    "name": str(var_name),
+                    "location": location,
+                    "dims": list(dims),
+                    "units": variable.attrs.get("units"),
+                    "standard_name": variable.attrs.get("standard_name"),
+                    "long_name": variable.attrs.get("long_name"),
+                }
+            )
         facts["variables"] = variables
 
     return facts
@@ -243,6 +276,7 @@ def get_capabilities(
     has_node_centered_vars = False
     has_edge_centered_vars = False
     face_centered_var_names: List[str] = []
+    face_centered_var_facts: List[Dict[str, Any]] = []
 
     if data_path is not None:
         for var_fact in facts.get("variables", []):
@@ -252,6 +286,7 @@ def get_capabilities(
             if location == "faces":
                 has_face_centered_vars = True
                 face_centered_var_names.append(var_name)
+                face_centered_var_facts.append(var_fact)
             elif location == "nodes":
                 has_node_centered_vars = True
             elif location == "edges":
@@ -319,6 +354,10 @@ def get_capabilities(
                 {
                     "name": var_name,
                     "location": location,
+                    "dims": var_fact.get("dims", []),
+                    "units": var_fact.get("units"),
+                    "standard_name": var_fact.get("standard_name"),
+                    "long_name": var_fact.get("long_name"),
                     "applicable_mcp_tools": applicable_mcp,
                     "applicable_uxarray_methods": applicable_uxarray,
                 }
@@ -495,8 +534,8 @@ def get_capabilities(
         ]
         if len(face_centered_var_names) >= 2:
             uxarray_capabilities["vector_calculus"] += [
-                "ux_da_u.curl(ux_da_v)  [two face-centered vars required]",
-                "ux_da_u.divergence(ux_da_v)  [two face-centered vars required]",
+                "ux_da_u.curl(ux_da_v)  [structurally applicable; verify component semantics]",
+                "ux_da_u.divergence(ux_da_v)  [structurally applicable; verify component semantics]",
             ]
 
     # Topological ops — for node/edge data
@@ -572,10 +611,18 @@ def get_capabilities(
         )
 
     if has_face_centered_vars and len(face_centered_var_names) >= 2:
-        recommendations.append(
-            "Multiple face-centered variables detected — you can compute vector operations "
-            "like curl and divergence between pairs (e.g. wind u/v components)."
-        )
+        if sum(_looks_vector_like(var) for var in face_centered_var_facts) >= 2:
+            recommendations.append(
+                "Multiple face-centered variables carry vector-like metadata. Curl and "
+                "divergence are structurally applicable; still verify that the selected "
+                "variables are paired components in the same coordinate system."
+            )
+        else:
+            recommendations.append(
+                "Multiple face-centered variables make curl/divergence structurally "
+                "computable, but fewer than two carry vector-like units or standard names. "
+                "Do not interpret a result physically until component identity is verified."
+            )
 
     result: Dict[str, Any] = {
         "grid_summary": grid_summary,
@@ -585,6 +632,30 @@ def get_capabilities(
             "configured_endpoints": config.endpoint_names,
         },
         "uxarray_capabilities": uxarray_capabilities,
+        "scientific_contracts": {
+            "vector_calculus": {
+                "structurally_applicable": len(face_centered_var_names) >= 2,
+                "semantic_suitability": (
+                    "supported_by_metadata"
+                    if sum(_looks_vector_like(var) for var in face_centered_var_facts)
+                    >= 2
+                    else "unverified"
+                ),
+                "evidence": {
+                    var["name"]: {
+                        "units": var.get("units"),
+                        "standard_name": var.get("standard_name"),
+                        "vector_like": _looks_vector_like(var),
+                    }
+                    for var in face_centered_var_facts
+                },
+                "required_for_physical_interpretation": [
+                    "distinct paired components",
+                    "compatible velocity/flux units",
+                    "documented component coordinate system",
+                ],
+            }
+        },
         "recommendations": recommendations,
     }
 
diff --git a/src/uxarray_mcp/tools/frontdoor.py b/src/uxarray_mcp/tools/frontdoor.py
index 6a2ccd9..7f7b804 100644
--- a/src/uxarray_mcp/tools/frontdoor.py
+++ b/src/uxarray_mcp/tools/frontdoor.py
@@ -64,7 +64,7 @@ def run_analysis(
     session_id: str | None = None,
     dataset_handle: str | None = None,
     result_name: str | None = None,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
     lat_spec: tuple | float | list[Any] | None = None,
@@ -85,8 +85,9 @@ def run_analysis(
     ``regrid_dataset``, ``remap_to_rectilinear``, ``temporal_mean``,
     ``anomaly``, ``ensemble_mean``, ``ensemble_spread``, and ``export``.
 
-    ``gradient`` and ``curl`` accept ``scale_by_radius`` (default False keeps the
-    historical unit-sphere result). ``gradient``, ``curl``, and ``divergence``
+    ``gradient`` and ``curl`` accept ``scale_by_radius`` (default True matches
+    UXarray and returns physical units when sphere-radius metadata exists).
+    ``gradient``, ``curl``, and ``divergence``
     also accept ``time_index``/``level_index`` to select a single time/level
     slice when the input variable(s) carry those extra dimensions (e.g. real
     model output shaped ``(time, lev, n_face)``); both default to 0 and are
diff --git a/src/uxarray_mcp/tools/orchestration.py b/src/uxarray_mcp/tools/orchestration.py
index 206f21e..66eaff0 100644
--- a/src/uxarray_mcp/tools/orchestration.py
+++ b/src/uxarray_mcp/tools/orchestration.py
@@ -204,7 +204,12 @@ def analyze_dataset(
 
     # ── Stage 5: zonal mean (needs data + face-centered variable) ───────────
     zonal_mean: Optional[dict[str, Any]] = None
-    if resolved_data is not None and selected_variable is not None:
+    validation_passed = validation is None or validation.get("passed") is not False
+    if (
+        resolved_data is not None
+        and selected_variable is not None
+        and validation_passed
+    ):
         zonal_mean = _safe_call(
             "calculate_zonal_mean",
             lambda: calculate_zonal_mean(
@@ -219,6 +224,10 @@ def analyze_dataset(
         )
         if zonal_mean is not None:
             stages_run.append("calculate_zonal_mean")
+    elif resolved_data is not None and selected_variable is not None:
+        warnings.append(
+            "calculate_zonal_mean: skipped because dataset validation failed."
+        )
 
     # ── Stage 6 + 7: plots (optional) ────────────────────────────────────────
     mesh_plot: Optional[dict[str, Any]] = None
@@ -239,7 +248,11 @@ def analyze_dataset(
             mesh_plot = _png_meta(plot_items)
             stages_run.append("plot_mesh")
 
-        if resolved_data is not None and selected_variable is not None:
+        if (
+            resolved_data is not None
+            and selected_variable is not None
+            and validation_passed
+        ):
             var_plot_items = _safe_call(
                 "plot_variable",
                 lambda: plot_variable(
@@ -255,6 +268,8 @@ def analyze_dataset(
             if var_plot_items is not None:
                 variable_plot = _png_meta(var_plot_items)
                 stages_run.append("plot_variable")
+        elif resolved_data is not None and selected_variable is not None:
+            warnings.append("plot_variable: skipped because dataset validation failed.")
 
     # ── Recommended next steps ──────────────────────────────────────────────
     next_steps: list[str] = []
@@ -307,7 +322,19 @@ def analyze_dataset(
     # (used by tests that import `inspect_mesh` directly from this module).
     _ = inspect_mesh  # keep import alive
 
-    venue = "hpc" if use_remote else "local"
+    stage_results = [mesh, variables, area, zonal_mean]
+    venues = {
+        stage.get("_provenance", {}).get("execution_venue")
+        for stage in stage_results
+        if isinstance(stage, dict)
+    }
+    venues.discard(None)
+    if not venues:
+        venue = "local"
+    elif len(venues) == 1:
+        venue = venues.pop()
+    else:
+        venue = "mixed:" + ",".join(sorted(venues))
     return attach_provenance(
         result,
         tool="analyze_dataset",
@@ -323,5 +350,14 @@ def analyze_dataset(
         },
         venue=venue,
         warnings=warnings,
+        validation_summary=(
+            {
+                "passed": validation.get("passed"),
+                "n_variables_checked": validation.get("n_variables_checked"),
+                "n_variables_failed": validation.get("n_variables_failed"),
+            }
+            if validation is not None
+            else None
+        ),
         selected_variable=selected_variable,
     )
diff --git a/src/uxarray_mcp/tools/scientific_agent.py b/src/uxarray_mcp/tools/scientific_agent.py
index 31e29ba..f7d786f 100644
--- a/src/uxarray_mcp/tools/scientific_agent.py
+++ b/src/uxarray_mcp/tools/scientific_agent.py
@@ -232,7 +232,7 @@ def run_scientific_agent(
         area_results = calculate_area(file_path)
 
     zonal_mean_results = None
-    if face_centered_var and data_path is not None:
+    if face_centered_var and data_path is not None and validation_passed:
         target_var = variable_name or face_centered_var
         reasoning_trace.append(
             {"stage": "execute", "action": f"calculate_zonal_mean({target_var!r})"}
diff --git a/src/uxarray_mcp/tools/vector_calc.py b/src/uxarray_mcp/tools/vector_calc.py
index b430c7c..12010ce 100644
--- a/src/uxarray_mcp/tools/vector_calc.py
+++ b/src/uxarray_mcp/tools/vector_calc.py
@@ -85,7 +85,7 @@ def calculate_gradient(
     grid_path: str,
     data_path: str,
     variable_name: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
     use_remote: bool = False,
@@ -107,8 +107,8 @@ def calculate_gradient(
         Name of the face-centered scalar variable.
     scale_by_radius : bool
         If True, divide unit-sphere derivatives by ``uxgrid.sphere_radius`` for
-        physical units (requires a grid with ``sphere_radius``). Default False
-        preserves the unit-sphere result. Local execution passes this to the
+        physical units (requires a grid with ``sphere_radius``). Default True
+        matches UXarray. Local execution passes this to the
         pinned UXarray directly. The remote worker, which may run an older
         UXarray, applies it capability-safely and falls back to the unit sphere
         if unsupported; the result reports the ``scale_by_radius`` actually
@@ -196,7 +196,7 @@ def calculate_curl(
     data_path: str,
     u_variable: str,
     v_variable: str,
-    scale_by_radius: bool = False,
+    scale_by_radius: bool = True,
     time_index: int = 0,
     level_index: int = 0,
     use_remote: bool = False,
@@ -225,8 +225,8 @@ def calculate_curl(
         Meridional (north-south) component, e.g. ``"uReconstructMeridional"``.
     scale_by_radius : bool
         If True, divide the unit-sphere result by ``uxgrid.sphere_radius`` for
-        physical units (requires a grid with ``sphere_radius``). Default False
-        preserves the unit-sphere result. Local execution passes this to the
+        physical units (requires a grid with ``sphere_radius``). Default True
+        matches UXarray. Local execution passes this to the
         pinned UXarray directly. The remote worker, which may run an older
         UXarray, applies it capability-safely and falls back to the unit sphere
         if unsupported; the result reports the ``scale_by_radius`` actually
diff --git a/tests/test_analyze_dataset.py b/tests/test_analyze_dataset.py
index 7158b2a..2f9b2c1 100644
--- a/tests/test_analyze_dataset.py
+++ b/tests/test_analyze_dataset.py
@@ -115,3 +115,84 @@ def test_analyze_dataset_recommended_next_steps_present(synthetic_mesh_with_data
         kw in joined
         for kw in ("plot_zonal_mean", "extract_cross_section", "subset_bbox")
     )
+
+
+def test_validation_failure_blocks_downstream_statistics(
+    monkeypatch, synthetic_mesh_with_data
+):
+    """Validation failure must prevent zonal statistics and variable plots."""
+    from uxarray_mcp.tools import inspection, remote_tools
+
+    grid_file, data_file = synthetic_mesh_with_data
+    monkeypatch.setattr(
+        inspection,
+        "validate_dataset",
+        lambda *args, **kwargs: {
+            "passed": False,
+            "n_variables_checked": 2,
+            "n_variables_failed": 1,
+            "_provenance": {"warnings": ["invalid values"]},
+        },
+    )
+    monkeypatch.setattr(
+        remote_tools,
+        "calculate_zonal_mean",
+        lambda *args, **kwargs: (_ for _ in ()).throw(
+            AssertionError("zonal mean must not run")
+        ),
+    )
+    monkeypatch.setattr(
+        remote_tools,
+        "plot_variable",
+        lambda *args, **kwargs: (_ for _ in ()).throw(
+            AssertionError("variable plot must not run")
+        ),
+    )
+
+    result = analyze_dataset(grid_file, data_file, include_plots=True)
+    assert result["zonal_mean"] is None
+    assert result["variable_plot"] is None
+    assert any("validation failed" in warning for warning in result["warnings"])
+    assert result["_provenance"]["validation_summary"]["passed"] is False
+
+
+def test_requested_remote_with_local_stages_reports_actual_venue(
+    monkeypatch, synthetic_mesh_with_data
+):
+    """Aggregate venue comes from stage provenance, not the use_remote request."""
+    from uxarray_mcp.tools import remote_tools
+
+    grid_file, data_file = synthetic_mesh_with_data
+
+    def force_local(result):
+        result["_provenance"]["execution_venue"] = "local"
+        return result
+
+    original_inspect = remote_tools.inspect_mesh
+    original_variables = remote_tools.inspect_variable
+    original_area = remote_tools.calculate_area
+    original_zonal = remote_tools.calculate_zonal_mean
+    monkeypatch.setattr(
+        remote_tools,
+        "inspect_mesh",
+        lambda *args, **kwargs: force_local(original_inspect(args[0])),
+    )
+    monkeypatch.setattr(
+        remote_tools,
+        "inspect_variable",
+        lambda *args, **kwargs: force_local(
+            original_variables(args[0], args[1], args[2])
+        ),
+    )
+    monkeypatch.setattr(
+        remote_tools,
+        "calculate_area",
+        lambda *args, **kwargs: force_local(original_area(args[0])),
+    )
+    monkeypatch.setattr(
+        remote_tools,
+        "calculate_zonal_mean",
+        lambda *args, **kwargs: force_local(original_zonal(args[0], args[1], args[2])),
+    )
+    result = analyze_dataset(grid_file, data_file, use_remote=True, include_plots=False)
+    assert result["_provenance"]["execution_venue"] == "local"
diff --git a/tests/test_capabilities.py b/tests/test_capabilities.py
index 42c623a..82ea34d 100644
--- a/tests/test_capabilities.py
+++ b/tests/test_capabilities.py
@@ -272,14 +272,44 @@ def test_vector_calculus_available_with_face_data(self, synthetic_mesh_with_data
         assert any("gradient" in m for m in vc)
         assert any("zonal_mean" in m for m in vc)
 
-    def test_curl_available_with_multiple_face_vars(self, synthetic_mesh_with_data):
-        """curl and divergence appear when 2+ face-centered variables exist."""
+    def test_curl_structural_capability_does_not_claim_semantic_suitability(
+        self, synthetic_mesh_with_data
+    ):
+        """Two scalar fields are computable inputs, not verified vector components."""
         grid_file, data_file = synthetic_mesh_with_data
         result = get_capabilities(grid_file, data_file)
         vc = result["uxarray_capabilities"]["vector_calculus"]
         # synthetic_mesh_with_data has temperature + pressure (2 face vars)
         assert any("curl" in m for m in vc)
         assert any("divergence" in m for m in vc)
+        contract = result["scientific_contracts"]["vector_calculus"]
+        assert contract["structurally_applicable"] is True
+        assert contract["semantic_suitability"] == "unverified"
+        assert contract["evidence"]["temperature"]["units"] == "K"
+        assert contract["evidence"]["temperature"]["vector_like"] is False
+        recommendations = " ".join(result["recommendations"])
+        assert "Do not interpret" in recommendations
+
+    def test_vector_metadata_support_is_reported(self, tmp_path, synthetic_mesh_file):
+        data_file = tmp_path / "wind.nc"
+        xr.Dataset(
+            {
+                "u": (
+                    ["nMesh2_face"],
+                    [1.0],
+                    {"units": "m/s", "standard_name": "eastward_wind"},
+                ),
+                "v": (
+                    ["nMesh2_face"],
+                    [2.0],
+                    {"units": "m/s", "standard_name": "northward_wind"},
+                ),
+            }
+        ).to_netcdf(data_file)
+        result = get_capabilities(synthetic_mesh_file, str(data_file))
+        contract = result["scientific_contracts"]["vector_calculus"]
+        assert contract["semantic_suitability"] == "supported_by_metadata"
+        assert contract["evidence"]["u"]["vector_like"] is True
 
     def test_remapping_available_with_face_data(self, synthetic_mesh_with_data):
         """Remapping methods listed without [needs face-centered data] note."""
diff --git a/tests/test_stateful_tools.py b/tests/test_stateful_tools.py
index 2978cee..5aac7b9 100644
--- a/tests/test_stateful_tools.py
+++ b/tests/test_stateful_tools.py
@@ -1,7 +1,16 @@
 """Tests for session, workflow, and persisted result tools."""
 
+from concurrent.futures import ThreadPoolExecutor
 from pathlib import Path
 
+import xarray as xr
+
+from uxarray_mcp.state import (
+    get_result,
+    persist_result,
+    save_result,
+    write_dataset_artifact,
+)
 from uxarray_mcp.tools import (
     create_session,
     get_result_handle,
@@ -100,3 +109,28 @@ def test_reset_session_state_clears_results_workflows_and_operations(
     assert state["workflow_ids"] == []
     assert state["operation_ids"] == []
     assert not Path(state_dir, "results", f"{workflow['result_handle']}.json").exists()
+
+
+def test_concurrent_artifact_writes_remain_readable(state_dir):
+    """Concurrent result persistence must never expose partial NetCDF files."""
+
+    def persist(index: int) -> str:
+        record = persist_result(
+            kind="concurrent-test",
+            name=f"result-{index}",
+            summary={"index": index},
+        )
+        dataset = xr.Dataset({"value": ("x", [float(index)])})
+        path = write_dataset_artifact(dataset, record["result_handle"])
+        record["artifact_path"] = path
+        save_result(record)
+        return record["result_handle"]
+
+    with ThreadPoolExecutor(max_workers=8) as executor:
+        handles = list(executor.map(persist, range(24)))
+
+    assert len(set(handles)) == 24
+    for handle in handles:
+        record = get_result(handle)
+        with xr.open_dataset(record["artifact_path"]) as dataset:
+            assert dataset["value"].size == 1
diff --git a/tests/test_vector_calc.py b/tests/test_vector_calc.py
index be44f84..5530046 100644
--- a/tests/test_vector_calc.py
+++ b/tests/test_vector_calc.py
@@ -269,12 +269,17 @@ def test_missing_units_warns(self, healpix_wind_dataset):
         # Synthetic fields carry no 'units' attr → non-velocity warning.
         result = compute_curl(healpix_wind_dataset, "u", "v")
         assert any("velocity" in w for w in result["component_warnings"])
+        assert result["scientific_status"]["status"] == "warning"
+        assert result["scientific_status"]["physically_interpretable"] is False
+        assert "VECTOR_UNITS_UNVERIFIED" in result["scientific_status"]["warning_codes"]
 
     def test_velocity_units_suppress_warning(self):
         # Fields with velocity units set at creation time → no warnings.
         ds = _make_wind_dataset(with_units=True)
-        result = compute_curl(ds, "u", "v")
+        result = compute_curl(ds, "u", "v", scale_by_radius=False)
         assert result["component_warnings"] == []
+        assert result["scientific_status"]["status"] == "complete"
+        assert result["scientific_status"]["physically_interpretable"] is True
 
     def test_curl_warning_reaches_provenance(self, monkeypatch):
         """The tool layer surfaces component warnings into _provenance.warnings."""
@@ -522,18 +527,18 @@ def test_accepts_use_remote_endpoint_session_params(self):
 
 
 # ---------------------------------------------------------------------------
-# scale_by_radius opt-in
+# scale_by_radius defaults align with UXarray
 # ---------------------------------------------------------------------------
 
 
 class TestScaleByRadius:
-    def test_gradient_default_keeps_unit_sphere(self, healpix_wind_dataset):
+    def test_gradient_default_requests_radius_scaling(self, healpix_wind_dataset):
         result = compute_gradient(healpix_wind_dataset, "temperature")
-        assert result["scale_by_radius"] is False
+        assert result["scale_by_radius"] is True
 
-    def test_curl_default_keeps_unit_sphere(self, healpix_wind_dataset):
+    def test_curl_default_requests_radius_scaling(self, healpix_wind_dataset):
         result = compute_curl(healpix_wind_dataset, "u", "v")
-        assert result["scale_by_radius"] is False
+        assert result["scale_by_radius"] is True
 
     def test_gradient_records_scale_by_radius_flag(self, healpix_wind_dataset):
         with warnings.catch_warnings():