diff --git a/ashvale/features.py b/ashvale/features.py
index f1e1cd9..968a4e5 100644
--- a/ashvale/features.py
+++ b/ashvale/features.py
@@ -82,7 +82,8 @@ def _rolling(a: np.ndarray, win: int, fn) -> np.ndarray:
 
 def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
                    press_slp: np.ndarray, lux: np.ndarray,
-                   grid_s: int, latitude: float, longitude: float
+                   grid_s: int, latitude: float, longitude: float,
+                   min_days_annual: float = 120.0
                    ) -> Tuple[np.ndarray, np.ndarray]:
     """Return (X of shape (n, N_FEATURES), valid mask of shape (n,))."""
     n = grid_ts.size
@@ -127,6 +128,20 @@ def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
     hour = (grid_ts % 86400.0) / 86400.0
     doy = (grid_ts % 31557600.0) / 31557600.0
 
+    # Annual harmonics are held at zero until the record spans enough of a year
+    # to excite them, exactly as the climatology fit already gates its annual
+    # terms. Left on from day one they are near-constant, near-collinear with
+    # each other and with the bias, and RLS answers that rank-deficient system
+    # with enormous cancelling weights. Measured on a real station after 1.5
+    # days: cos_doy +1191, sin_doy +1174, ||theta|| 1680 against a median |theta|
+    # of 1.67, cond(P) 3.1e9, and a six hour forecast of 53 C in a 24 C room.
+    # Zero is the honest value: with a day and a half of data the station knows
+    # nothing whatsoever about the season.
+    span_days = float(grid_ts[-1] - grid_ts[0]) / 86400.0 if n > 1 else 0.0
+    annual_on = 1.0 if span_days >= min_days_annual else 0.0
+    sin_doy = np.sin(2 * np.pi * doy) * annual_on
+    cos_doy = np.cos(2 * np.pi * doy) * annual_on
+
     X = np.column_stack([
         np.ones(n),
         temp, temp_rate_1h, temp_rate_3h, temp_std_3h, temp_dev_24h,
@@ -136,7 +151,7 @@ def build_features(grid_ts: np.ndarray, temp: np.ndarray, hum: np.ndarray,
         log_lux, cloud, elev, np.clip(elev, 0.0, None), (elev > 0.0).astype(float),
         np.sin(2 * np.pi * hour), np.cos(2 * np.pi * hour),
         np.sin(4 * np.pi * hour), np.cos(4 * np.pi * hour),
-        np.sin(2 * np.pi * doy), np.cos(2 * np.pi * doy),
+        sin_doy, cos_doy,
         press_anom * (hum - 70.0) / 100.0,
         press_tend_3h * dep,
     ])
@@ -172,7 +187,12 @@ class Standardiser:
         if self.n < 2:
             return np.atleast_2d(X)
         std = np.sqrt(self.m2 / max(self.n - 1, 1))
-        std = np.where(std < 1e-8, 1.0, std)
+        # 1e-8 was a token guard: it only catches a bit-exactly constant column.
+        # A feature that merely barely moves sails through and gets divided by
+        # its own noise, which manufactures a large z-score out of nothing. A
+        # feature with this little spread carries no information, so scale it by
+        # one and let it stay near zero rather than amplifying it.
+        std = np.where(std < 1e-3, 1.0, std)
         out = (np.atleast_2d(X) - self.mean) / std
         out[:, 0] = 1.0            # keep the bias column intact
         return out
diff --git a/ashvale/models/nowcast.py b/ashvale/models/nowcast.py
index 9593f33..3b00161 100644
--- a/ashvale/models/nowcast.py
+++ b/ashvale/models/nowcast.py
@@ -179,6 +179,18 @@ class NowcastEnsemble:
         """Batch pass over history. Called on startup and every retrain tick."""
         if X.shape[0] < 10:
             return {"rows": 0}
+        # A refit starts from the prior. Without this, every retrain tick replays
+        # the same history into a live filter, and RLS with forgetting reads that
+        # as new evidence each time: measured on a real station after 1.5 days,
+        # 453 grid rows had produced 64,676 updates, cond(P) of 3.1e9 and a
+        # weight vector of norm 1680 whose two largest entries were the annual
+        # harmonics the record cannot yet resolve. The result was a six hour
+        # forecast of 53 C in a 24 C room, with a plus or minus of 0.43.
+        #
+        # The conformal calibrators and the Hedge weights are deliberately left
+        # alone: those are earned from scored forecasts, not from this regression.
+        for head in self.heads.values():
+            head.model.reset()
         self.scaler.partial_fit(X[valid][:: max(1, X.shape[0] // 2000)])
         Xs = self.scaler.transform(X)
 
diff --git a/ashvale/models/rls.py b/ashvale/models/rls.py
index 0d89719..44680d1 100644
--- a/ashvale/models/rls.py
+++ b/ashvale/models/rls.py
@@ -47,8 +47,9 @@ class RecursiveLeastSquares:
         self.d = int(n_features)
         self.lam = float(forgetting)
         self.p_max = float(p_max)
+        self.delta = float(delta)          # kept so a refit can return to the prior
         self.theta = np.zeros(self.d)
-        self.P = np.eye(self.d) * float(delta)
+        self.P = np.eye(self.d) * self.delta
         self.n_updates = 0
         self.ewma_sq_error = 0.0
 
@@ -102,14 +103,32 @@ class RecursiveLeastSquares:
     def noise_var(self) -> float:
         return float(max(self.ewma_sq_error, 1e-9))
 
+    def reset(self) -> None:
+        """Return to the prior, keeping the configuration.
+
+        A batch refit has to start from here rather than continuing, because
+        replaying the same history into a live filter is not the same as seeing
+        new data. RLS with forgetting treats every update as fresh evidence, so
+        feeding it the same rows on each retrain tick makes it believe it has
+        many times the data it has: P collapses, and the weights in directions
+        the data never excites drift without anything to pull them back.
+        """
+        self.theta = np.zeros(self.d)
+        self.P = np.eye(self.d) * self.delta
+        self.n_updates = 0
+        self.ewma_sq_error = 0.0
+
     def to_dict(self) -> Dict:
         return {"d": self.d, "lam": self.lam, "p_max": self.p_max,
+                "delta": self.delta,
                 "theta": self.theta.tolist(), "P": self.P.tolist(),
                 "n": self.n_updates, "ewma": self.ewma_sq_error}
 
     @classmethod
     def from_dict(cls, s: Dict) -> "RecursiveLeastSquares":
-        m = cls(s["d"], s["lam"], 1.0, s.get("p_max", 1e6))
+        # delta must survive the round trip or a refit after a restart would
+        # return to the wrong prior.
+        m = cls(s["d"], s["lam"], s.get("delta", 100.0), s.get("p_max", 1e6))
         m.theta = np.array(s["theta"], dtype=float)
         m.P = np.array(s["P"], dtype=float)
         m.n_updates = s.get("n", 0)
diff --git a/ashvale/station.py b/ashvale/station.py
index 923fc1f..50c46cd 100644
--- a/ashvale/station.py
+++ b/ashvale/station.py
@@ -505,6 +505,7 @@ class Station:
             grid_ts, cols["temperature"], cols["humidity"], cols["pressure"],
             cols["lux"], self.cfg.model.grid_s,
             self.cfg.site.latitude, self.cfg.site.longitude,
+            self.cfg.model.climatology_min_days_annual,
         )
         return grid_ts, cols, X, valid
 
diff --git a/scripts/evaluate.py b/scripts/evaluate.py
index bbdfafc..73cb0a0 100644
--- a/scripts/evaluate.py
+++ b/scripts/evaluate.py
@@ -86,7 +86,8 @@ def main() -> None:
     )
     X, valid = build_features(grid_ts, cols["temperature"], cols["humidity"],
                               cols["pressure"], cols["lux"], cfg.model.grid_s,
-                              cfg.site.latitude, cfg.site.longitude)
+                              cfg.site.latitude, cfg.site.longitude,
+                              cfg.model.climatology_min_days_annual)
 
     n = grid_ts.size
     split = int(n * args.train_frac)
diff --git a/tests/test_models.py b/tests/test_models.py
index 7242550..fb34344 100644
--- a/tests/test_models.py
+++ b/tests/test_models.py
@@ -149,3 +149,94 @@ def test_zambretti_rain_prior_rises_with_z():
     settled = zambretti(1035.0, 1.5, 6)
     stormy = zambretti(960.0, -2.5, 6)
     assert stormy["prior_rain_prob"] > settled["prior_rain_prob"]
+
+
+# ---------------------------------------------------------------- refit safety
+
+def test_rls_reset_returns_to_the_prior():
+    m = RecursiveLeastSquares(n_features=5, forgetting=0.999, delta=100.0)
+    rng = np.random.default_rng(9)
+    for _ in range(500):
+        m.update(rng.normal(size=5), float(rng.normal()))
+    assert m.n_updates == 500
+    m.reset()
+    assert m.n_updates == 0
+    assert np.allclose(m.theta, 0.0)
+    assert np.allclose(m.P, np.eye(5) * 100.0)
+
+
+def test_rls_delta_survives_serialisation():
+    """A refit after a restart must return to the same prior it started from."""
+    m = RecursiveLeastSquares(n_features=4, forgetting=0.99, delta=100.0)
+    m.update(np.ones(4), 1.0)
+    back = RecursiveLeastSquares.from_dict(m.to_dict())
+    back.reset()
+    assert np.allclose(back.P, np.eye(4) * 100.0), "reload lost the prior"
+
+
+def test_repeated_refits_do_not_accumulate():
+    """Refitting the same history must be idempotent, not cumulative.
+
+    This is the bug that put a 53 C six-hour forecast on a real station in a
+    24 C room. fit() replayed history into a live filter on every retrain tick
+    and never reset, so 453 grid rows had produced 64,676 updates in a day and a
+    half. RLS with forgetting reads each update as fresh evidence, so P
+    collapsed and the weights drifted without bound in the directions the data
+    never excited.
+    """
+    from ashvale.config import CONFIG
+    from ashvale.models.nowcast import NowcastEnsemble
+
+    rng = np.random.default_rng(3)
+    n = 400
+    g = CONFIG.model.grid_s
+    ts = np.arange(n) * g + 1.7554e9
+    cols = {
+        "temperature": 22 + 2 * np.sin(np.arange(n) / 40.0) + 0.2 * rng.normal(size=n),
+        "humidity": 50 + 5 * np.cos(np.arange(n) / 33.0),
+        "pressure": 1013 + np.sin(np.arange(n) / 77.0),
+        "lux": np.clip(300 * np.sin(np.arange(n) / 120.0), 0, None),
+    }
+    X = rng.normal(size=(n, 33))
+    X[:, 0] = 1.0
+    valid = np.ones(n, dtype=bool)
+
+    ens = NowcastEnsemble(CONFIG.model.targets, CONFIG.model.horizons_s, CONFIG.model)
+    head = ens.heads[("temperature", 21600)]
+
+    ens.fit(X, valid, cols, None, ts)
+    first_norm = float(np.linalg.norm(head.model.theta))
+    first_updates = head.model.n_updates
+
+    for _ in range(15):
+        ens.fit(X, valid, cols, None, ts)
+
+    assert head.model.n_updates == first_updates, "updates accumulated across refits"
+    assert float(np.linalg.norm(head.model.theta)) == pytest.approx(first_norm, rel=0.05)
+
+
+def test_annual_harmonics_are_zero_until_the_record_spans_a_season():
+    """Two near-constant, near-collinear columns are a rank-deficient regressor.
+
+    Left on from day one, sin_doy and cos_doy carried +1174 and +1191 on a real
+    station whose median weight was 1.67. Zero is the honest value: a day and a
+    half of data says nothing whatsoever about the season.
+    """
+    from ashvale.features import FEATURE_NAMES, build_features
+
+    n = 450
+    ts = np.arange(n) * 300.0 + 1.7554e9          # about 1.5 days
+    t = 22 + 2 * np.sin(np.arange(n) / 40.0)
+    h = 50 + 5 * np.cos(np.arange(n) / 33.0)
+    p = 1013 + np.sin(np.arange(n) / 77.0)
+    lux = np.clip(300 * np.sin(np.arange(n) / 120.0), 0, None)
+
+    si, ci = FEATURE_NAMES.index("sin_doy"), FEATURE_NAMES.index("cos_doy")
+
+    X, _ = build_features(ts, t, h, p, lux, 300, 52.2, 0.12, min_days_annual=120.0)
+    assert np.all(X[:, si] == 0.0) and np.all(X[:, ci] == 0.0)
+
+    # A record that does span the year keeps them.
+    ts_long = np.arange(n) * (200 * 86400.0 / n) + 1.7554e9
+    X2, _ = build_features(ts_long, t, h, p, lux, 300, 52.2, 0.12, min_days_annual=120.0)
+    assert X2[:, si].std() > 0.1, "annual terms should return once the record is long enough"