diff --git a/CHANGELOG.md b/CHANGELOG.md
index ee7b0f5..bad9034 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,23 @@
 # Changelog
 
+## [Unreleased]
+
+- **Fix:** `Tub` kept its counter in `_step`, the attribute where dissmodel's
+  `Model` stores the interval between ticks, so the interval grew at every
+  execution (a run from t=1 to 30 executed 8 times). The counter is now
+  `_tick`; a test checks that every model executes once per tick.
+- Examples, notebooks, demo and docs use `Environment(start_time=1, end_time=N)`,
+  which runs N steps at t = 1..N like a TerraME `Timer` with `finalTime = N`.
+  The default `start_time=0` runs one step more and shifts the clock by one hour
+  in `RoomTemperature`. A note in the docs explains the convention.
+- `tests/test_terrame_parity.py` compares the full trajectory of each model
+  with a step-by-step transcription of the TerraME `sysdyn` package
+  (`terrame/sysdyn`, commit 9e671cb). Eleven models match within 1e-9.
+  Two known differences are kept as strict `xfail`: `MonoLake` (TerraME's
+  sampled Catmull-Rom spline vs scipy `CubicSpline`) and `Daisyworld`
+  (TerraME recomputes the planet albedo after updating the white area; same
+  equilibrium, different trajectory). `RandomWalk` is not compared.
+
 ## [0.1.1] — 2026-10-08
 
 First citable release. The `v0.1.0` tag was created before the release
diff --git a/demo/app.py b/demo/app.py
index 17b6e81..98fdf51 100644
--- a/demo/app.py
+++ b/demo/app.py
@@ -64,7 +64,7 @@
 # ---------------------------------------------------------------------------
 
 # 1. Environment
-env = Environment(start_time=0, end_time=steps)
+env = Environment(start_time=1, end_time=steps)
 
 # 2. Model
 model = ModelClass()
diff --git a/docs/architecture.md b/docs/architecture.md
index 84094ff..3ae0317 100644
--- a/docs/architecture.md
+++ b/docs/architecture.md
@@ -76,7 +76,7 @@ class Decay(Model):
         self.amount -= decay
 
 # Running the simulation
-env = Environment(end_time=50)
+env = Environment(start_time=1, end_time=50)
 model = Decay()
 env.run()
 ```
diff --git a/docs/models.md b/docs/models.md
index b58e75a..b8fe103 100644
--- a/docs/models.md
+++ b/docs/models.md
@@ -17,9 +17,16 @@ from dissmodel_sysdyn.models import SIR
 import matplotlib.pyplot as plt
 
 # Create the simulation environment
-env = Environment(end_time=60)
+env = Environment(start_time=1, end_time=60)
 ```
 
+!!! note "Counting steps as in TerraME"
+    `end_time` is inclusive. With `start_time=1, end_time=N` the model runs
+    N steps at t = 1, ..., N, like a TerraME `Timer` with `finalTime = N`.
+    The `Environment` default, `start_time=0`, runs one step more (t = 0, ..., N)
+    and shifts the clock by one, which matters for models that read `env.now()`,
+    such as `RoomTemperature`.
+
 ### Step 2: Model Configuration
 We instantiate the SIR model. We start with a nearly susceptible population and a small number of "Patient Zeros".
 
@@ -76,7 +83,7 @@ from dissmodel.core import Environment
 from dissmodel_sysdyn.models import PredatorPrey
 import matplotlib.pyplot as plt
 
-env = Environment(end_time=200)
+env = Environment(start_time=1, end_time=200)
 ```
 
 ### Step 2: Define Biological Parameters
diff --git a/docs/reference.md b/docs/reference.md
index 50c308c..333e4ad 100644
--- a/docs/reference.md
+++ b/docs/reference.md
@@ -29,7 +29,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import ChaoticGrowth
 
-env = Environment(end_time=100)
+env = Environment(start_time=1, end_time=100)
 model = ChaoticGrowth(pop=0.1, rate=3.9)
 env.run()
 ```
@@ -62,7 +62,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Coffee
 
-env = Environment(end_time=30)
+env = Environment(start_time=1, end_time=30)
 coffee = Coffee(temperature=90, room_temperature=22, cooling_rate=0.15)
 env.run()
 ```
@@ -105,7 +105,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Daisyworld
 
-env = Environment(end_time=200)
+env = Environment(start_time=1, end_time=200)
 model = Daisyworld(sun_luminosity=1.0)
 env.run()
 ```
@@ -138,7 +138,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Homeostasis
 
-env = Environment(end_time=50)
+env = Environment(start_time=1, end_time=50)
 model = Homeostasis(stock=10, gain=5, rate=-0.2)
 env.run()
 ```
@@ -171,7 +171,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import LimitedGrowth
 
-env = Environment(end_time=100)
+env = Environment(start_time=1, end_time=100)
 model = LimitedGrowth(pop=10, rate=0.2, capacity=1000)
 env.run()
 ```
@@ -211,7 +211,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Lorenz
 
-env = Environment(end_time=5000)
+env = Environment(start_time=1, end_time=5000)
 model = Lorenz()
 env.run()
 ```
@@ -245,7 +245,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import MonoLake
 
-env = Environment(end_time=50)
+env = Environment(start_time=1, end_time=50)
 model = MonoLake(export=150)
 env.run()
 ```
@@ -279,7 +279,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import PopulationGrowth
 
-env = Environment(end_time=20)
+env = Environment(start_time=1, end_time=20)
 model = PopulationGrowth(population=100, growth=0.05, growth_change=1.1)
 env.run()
 ```
@@ -311,7 +311,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import PredatorPrey
 
-env = Environment(end_time=200)
+env = Environment(start_time=1, end_time=200)
 model = PredatorPrey()
 env.run()
 ```
@@ -347,7 +347,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import RandomWalk
 
-env = Environment(end_time=100)
+env = Environment(start_time=1, end_time=100)
 model = RandomWalk(prob=0.55, seed=42)
 env.run()
 ```
@@ -383,7 +383,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import RoomTemperature
 
-env = Environment(end_time=24)
+env = Environment(start_time=1, end_time=24)
 model = RoomTemperature()
 env.run()
 ```
@@ -420,7 +420,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import SIR
 
-env = Environment(end_time=60)
+env = Environment(start_time=1, end_time=60)
 sir = SIR(infected=5, probability=0.3)
 env.run()
 ```
@@ -454,7 +454,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Tub
 
-env = Environment(end_time=50)
+env = Environment(start_time=1, end_time=50)
 tub = Tub(water=100, out_flow=2, in_flow=50, in_period=20)
 env.run()
 ```
@@ -486,7 +486,7 @@ def execute(self) -> None:
 from dissmodel.core import Environment
 from dissmodel_sysdyn.models import Yeast
 
-env = Environment(end_time=15)
+env = Environment(start_time=1, end_time=15)
 yeast = Yeast()
 env.run()
 ```
diff --git a/examples/cli/sysdyn_sir.py b/examples/cli/sysdyn_sir.py
index 4ff160d..3845df4 100644
--- a/examples/cli/sysdyn_sir.py
+++ b/examples/cli/sysdyn_sir.py
@@ -26,7 +26,7 @@
 # ---------------------------------------------------------------------------
 
 # 1. Environment
-env = Environment()
+env = Environment(start_time=1, end_time=30)
 
 # 2. Model
 SIR(
@@ -44,4 +44,4 @@
 # ---------------------------------------------------------------------------
 # Run
 # ---------------------------------------------------------------------------
-env.run(30)
+env.run()
diff --git a/examples/notebooks/sysdyn_chaotic_growth.ipynb b/examples/notebooks/sysdyn_chaotic_growth.ipynb
index fc0e85b..05d9cb8 100644
--- a/examples/notebooks/sysdyn_chaotic_growth.ipynb
+++ b/examples/notebooks/sysdyn_chaotic_growth.ipynb
@@ -65,7 +65,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=100)"
    ]
   },
   {
@@ -147,7 +147,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(100)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_coffee.ipynb b/examples/notebooks/sysdyn_coffee.ipynb
index 1d04a2c..2fba009 100644
--- a/examples/notebooks/sysdyn_coffee.ipynb
+++ b/examples/notebooks/sysdyn_coffee.ipynb
@@ -62,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=50)"
    ]
   },
   {
@@ -147,7 +147,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(50)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_daisyworld.ipynb b/examples/notebooks/sysdyn_daisyworld.ipynb
index e2e8ac4..5a08911 100644
--- a/examples/notebooks/sysdyn_daisyworld.ipynb
+++ b/examples/notebooks/sysdyn_daisyworld.ipynb
@@ -65,7 +65,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=100)"
    ]
   },
   {
@@ -178,7 +178,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(100)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_homeostasis.ipynb b/examples/notebooks/sysdyn_homeostasis.ipynb
index a557647..a3d26a3 100644
--- a/examples/notebooks/sysdyn_homeostasis.ipynb
+++ b/examples/notebooks/sysdyn_homeostasis.ipynb
@@ -62,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=30)"
    ]
   },
   {
@@ -147,7 +147,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(30)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_limited_growth.ipynb b/examples/notebooks/sysdyn_limited_growth.ipynb
index 588e348..3c28b73 100644
--- a/examples/notebooks/sysdyn_limited_growth.ipynb
+++ b/examples/notebooks/sysdyn_limited_growth.ipynb
@@ -62,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=100)"
    ]
   },
   {
@@ -147,7 +147,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(100)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_lorenz.ipynb b/examples/notebooks/sysdyn_lorenz.ipynb
index 4e4b173..9c87f2e 100644
--- a/examples/notebooks/sysdyn_lorenz.ipynb
+++ b/examples/notebooks/sysdyn_lorenz.ipynb
@@ -64,7 +64,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=2000)"
    ]
   },
   {
@@ -150,7 +150,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(2000)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_mono_lake.ipynb b/examples/notebooks/sysdyn_mono_lake.ipynb
index d372dbe..3d66b34 100644
--- a/examples/notebooks/sysdyn_mono_lake.ipynb
+++ b/examples/notebooks/sysdyn_mono_lake.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=50)"
    ]
   },
   {
@@ -156,7 +156,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(50)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_population_growth.ipynb b/examples/notebooks/sysdyn_population_growth.ipynb
index e2d8bcb..b156e25 100644
--- a/examples/notebooks/sysdyn_population_growth.ipynb
+++ b/examples/notebooks/sysdyn_population_growth.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=20)"
    ]
   },
   {
@@ -142,7 +142,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(20)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_predator_prey.ipynb b/examples/notebooks/sysdyn_predator_prey.ipynb
index faa61d6..9c9c487 100644
--- a/examples/notebooks/sysdyn_predator_prey.ipynb
+++ b/examples/notebooks/sysdyn_predator_prey.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=500)"
    ]
   },
   {
@@ -146,7 +146,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(500)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_random_walk.ipynb b/examples/notebooks/sysdyn_random_walk.ipynb
index c795c63..cfd230b 100644
--- a/examples/notebooks/sysdyn_random_walk.ipynb
+++ b/examples/notebooks/sysdyn_random_walk.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=100)"
    ]
   },
   {
@@ -145,7 +145,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(100)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_room_temperature.ipynb b/examples/notebooks/sysdyn_room_temperature.ipynb
index 084736c..97eef2e 100644
--- a/examples/notebooks/sysdyn_room_temperature.ipynb
+++ b/examples/notebooks/sysdyn_room_temperature.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=24)"
    ]
   },
   {
@@ -149,7 +149,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(24)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_sir.ipynb b/examples/notebooks/sysdyn_sir.ipynb
index e7c8a1c..87d96cc 100644
--- a/examples/notebooks/sysdyn_sir.ipynb
+++ b/examples/notebooks/sysdyn_sir.ipynb
@@ -68,7 +68,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=30)"
    ]
   },
   {
@@ -156,7 +156,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(30)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_tub.ipynb b/examples/notebooks/sysdyn_tub.ipynb
index cced38f..4f0c424 100644
--- a/examples/notebooks/sysdyn_tub.ipynb
+++ b/examples/notebooks/sysdyn_tub.ipynb
@@ -62,7 +62,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=100)"
    ]
   },
   {
@@ -144,7 +144,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(100)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/notebooks/sysdyn_yeast.ipynb b/examples/notebooks/sysdyn_yeast.ipynb
index 8c95efb..ba938d7 100644
--- a/examples/notebooks/sysdyn_yeast.ipynb
+++ b/examples/notebooks/sysdyn_yeast.ipynb
@@ -63,7 +63,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env = Environment()"
+    "env = Environment(start_time=1, end_time=15)"
    ]
   },
   {
@@ -150,7 +150,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "env.run(15)"
+    "env.run()"
    ]
   },
   {
diff --git a/examples/streamlit/sysdyn_all.py b/examples/streamlit/sysdyn_all.py
index 6f53708..2bd4cad 100644
--- a/examples/streamlit/sysdyn_all.py
+++ b/examples/streamlit/sysdyn_all.py
@@ -61,7 +61,7 @@
 # ---------------------------------------------------------------------------
 
 # 1. Environment
-env = Environment(start_time=0, end_time=steps)
+env = Environment(start_time=1, end_time=steps)
 
 # 2. Model — instantiated dynamically from sidebar selection
 ModelClass = model_classes[model_name]
diff --git a/examples/streamlit/sysdyn_sir.py b/examples/streamlit/sysdyn_sir.py
index b4b986f..22f95b7 100644
--- a/examples/streamlit/sysdyn_sir.py
+++ b/examples/streamlit/sysdyn_sir.py
@@ -42,7 +42,7 @@
 # ---------------------------------------------------------------------------
 
 # 1. Environment
-env = Environment(start_time=0, end_time=steps)
+env = Environment(start_time=1, end_time=steps)
 
 # 2. Model
 sir = SIR(
diff --git a/src/dissmodel_sysdyn/models/tub.py b/src/dissmodel_sysdyn/models/tub.py
index b2e4007..53c6273 100644
--- a/src/dissmodel_sysdyn/models/tub.py
+++ b/src/dissmodel_sysdyn/models/tub.py
@@ -39,7 +39,7 @@ class Tub(Model):
     Examples
     --------
     >>> from dissmodel.core import Environment
-    >>> env = Environment(end_time=8)
+    >>> env = Environment(start_time=1, end_time=30)
     >>> tub = Tub(water=40, out_flow=5, in_flow=20, in_period=10)
     >>> env.run()
     """
@@ -81,7 +81,10 @@ def setup(
         self.out_flow  = out_flow
         self.in_flow   = in_flow
         self.in_period = in_period
-        self._step     = 0
+        # Executions counter. Must not be named ``_step``: dissmodel's Model
+        # keeps the scheduling interval in ``_step``, so incrementing it here
+        # made the interval grow (ticks 1, 2, 4, 7, 11, ...).
+        self._tick     = 0
 
     def execute(self) -> None:
         """
@@ -90,13 +93,13 @@ def execute(self) -> None:
         Applies the outflow every step and the inflow every
         ``in_period`` steps. Water is clamped to zero from below.
         """
-        self._step += 1
+        self._tick += 1
 
         self.water -= self.out_flow
         if self.water < 0.0:
             self.water = 0.0
 
-        if self._step % self.in_period == 0:
+        if self._tick % self.in_period == 0:
             self.water += self.in_flow
 
 
diff --git a/tests/test_models_smoke.py b/tests/test_models_smoke.py
index c06657b..09ace54 100644
--- a/tests/test_models_smoke.py
+++ b/tests/test_models_smoke.py
@@ -44,3 +44,26 @@ def test_model_runs(name, cls):
     assert numeric, f"{name} exposes no numeric state"
     for key, value in numeric.items():
         assert math.isfinite(value), f"{name}.{key} = {value}"
+
+
+@pytest.mark.parametrize("name, cls", MODELS, ids=[n for n, _ in MODELS])
+def test_model_executes_once_per_tick(name, cls, monkeypatch):
+    """A model must not change its own scheduling.
+
+    ``dissmodel.core.Model`` keeps the interval between ticks in ``_step``;
+    a model that reuses that name (as ``Tub`` did up to 0.1.1) silently
+    skips ticks.
+    """
+    ticks: list[float] = []
+    original = cls.execute
+
+    def counting(self):
+        ticks.append(self.env.now())
+        original(self)
+
+    monkeypatch.setattr(cls, "execute", counting)
+    np.random.seed(0)
+    env = Environment(start_time=1, end_time=STEPS)
+    cls()
+    env.run()
+    assert ticks == list(range(1, STEPS + 1))
diff --git a/tests/test_terrame_parity.py b/tests/test_terrame_parity.py
new file mode 100644
index 0000000..eeff55d
--- /dev/null
+++ b/tests/test_terrame_parity.py
@@ -0,0 +1,287 @@
+"""
+Parity with the TerraME ``sysdyn`` package.
+
+Each ``ref_*`` function below transcribes, step by step, the update rule of a
+model in https://github.com/TerraME/sysdyn (``lua/<Model>.lua``, commit
+9e671cb, 2020-08-20), with the TerraME default parameters and ``finalTime``.
+A TerraME ``Timer`` runs t = 1, ..., finalTime, so the DisSModel model runs in
+``Environment(start_time=1, end_time=finalTime)``.
+
+The test compares the full trajectory, not only the final value. The reference
+is a transcription, not output of TerraME itself: it checks that the Python
+models follow the Lua equations and the TerraME clock.
+
+Not covered: ``RandomWalk`` (different random generators); the SIR
+``maximum`` event and the PredatorPrey ``period`` option (not ported, and
+inactive with the default values).
+
+Known differences, kept as strict ``xfail`` so that a fix is noticed:
+
+* ``MonoLake``: TerraME's ``sci.Spline`` is a Catmull-Rom curve sampled at
+  ``steps`` points per segment and read by linear interpolation; the Python
+  model uses ``scipy.interpolate.CubicSpline``. Up to 26 KAF of difference in
+  50 steps.
+* ``Daisyworld``: the Lua model recomputes the planet albedo for the black
+  daisies after updating the white area; the Python model uses the albedo
+  from the start of the step. The trajectories differ (up to 0.61 in area)
+  and reach the same equilibrium.
+"""
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import dissmodel_sysdyn.models as sd
+from dissmodel.core import Environment
+
+RTOL = 1e-9
+ATOL = 1e-9
+
+
+def run_dissmodel(cls, final_time, attrs, **params):
+    """Run ``cls`` for t = 1..final_time and return the state after each step."""
+    env = Environment(start_time=1, end_time=final_time)
+    model = cls(**params)  # noqa: F841  (connects to env on creation)
+    trace = []
+    original = cls.execute
+
+    def recording(self):
+        original(self)
+        trace.append([getattr(self, a) for a in attrs])
+
+    # record the state after each execute(), for this run only
+    try:
+        cls.execute = recording
+        env.run()
+    finally:
+        cls.execute = original
+    return np.array(trace, dtype=float)
+
+
+# ---------------------------------------------------------------------------
+# Transcriptions of terrame/sysdyn lua/*.lua
+# ---------------------------------------------------------------------------
+
+def ref_coffee(final_time=20):
+    temperature, room, out = 80.0, 20.0, []
+    for _ in range(final_time):
+        temperature = temperature - (temperature - room) * 0.1
+        out.append([temperature])
+    return out
+
+
+def ref_chaotic_growth(final_time=300):
+    pop, rate, out = 0.10, 4.0, []
+    for _ in range(final_time):
+        pop = rate * pop * (1 - pop)
+        out.append([pop])
+    return out
+
+
+def ref_homeostasis(final_time=30):
+    stock, gain, rate, out = 0.0, 2.0, -0.4, []
+    for _ in range(final_time):
+        stock = stock + gain + rate * stock
+        out.append([stock])
+    return out
+
+
+def ref_limited_growth(final_time=100):
+    pop, rate, capacity, out = 300.0, 0.105, 20000.0, []
+    for _ in range(final_time):
+        pop = pop + pop * rate * (1 - pop / capacity)
+        out.append([pop])
+    return out
+
+
+def ref_lorenz(final_time=10000):
+    x = y = z = 1.0
+    delta, rho, sigma, beta = 0.01, 28.0, 10.0, 8.0 / 3.0
+    out = []
+    for _ in range(final_time):
+        nx = sigma * (y - x)
+        ny = x * (rho - z) - y
+        nz = x * y - beta * z
+        x, y, z = x + delta * nx, y + delta * ny, z + delta * nz
+        out.append([x, y, z])
+    return out
+
+
+def ref_population_growth(final_time=7):
+    population, growth, growth_change, out = 60.0, 0.5, 1.0, []
+    for _ in range(final_time):
+        population = population * (1 + growth)
+        growth = growth * growth_change
+        out.append([population])
+    return out
+
+
+def ref_predator_prey(final_time=500):
+    predator, prey, out = 40.0, 1000.0, []
+    for _ in range(final_time):
+        prey = prey + 0.08 * prey - 0.001 * prey * predator
+        predator = predator - 0.02 * predator + 0.00002 * prey * predator  # uses the new prey
+        out.append([prey, predator])
+    return out
+
+
+def ref_room_temperature(final_time=24):
+    temp_set, inside, out = 20.0, 15.0, []
+    for t in range(1, final_time + 1):
+        outside = 15 - 0.1 * (t - 12) ** 2
+        inflow = 0.33 * (temp_set - inside)
+        outflow = 0.30 * (inside - outside)
+        inside = inside + inflow - outflow
+        out.append([inside, outside])
+    return out
+
+
+def ref_sir(final_time=30):
+    s, i, r = 9998.0, 2.0, 0.0
+    duration, contacts, probability = 2, 6, 0.25
+    out = []
+    for _ in range(final_time):
+        proportion = s / (s + i + r)
+        new_infected = i * contacts * probability * proportion
+        new_recovered = i / duration
+        s = s - new_infected
+        r = r + new_recovered
+        i = i + new_infected - new_recovered
+        out.append([s, i, r])
+    return out
+
+
+def ref_tub(final_time=8, in_flow=0.0):
+    water, out_flow, out = 40.0, 5.0, []
+    for t in range(1, final_time + 1):
+        water = water - out_flow
+        if water < 0:
+            water = 0.0
+        if t >= 10 and t % 10 == 0:        # Event{start = 10, period = 10}
+            water = water + in_flow
+        out.append([water])
+    return out
+
+
+def ref_yeast(final_time=9):
+    cells, capacity, rate, out = 9.6, 665.0, 1.1, []
+    for _ in range(final_time):
+        cells = cells + cells * rate * (1 - cells / capacity)
+        if cells > capacity:
+            cells = capacity
+        out.append([cells])
+    return out
+
+
+_SIGMA, _SOLAR_FLUX, _HEAT = 5.67e-08, 3668.0, 2.06e09
+
+
+def _daisy_growth(temp_k):
+    temp = temp_k - 273
+    return 1 - 0.003265 * (22.5 - temp) ** 2 if 5.0 < temp < 40.0 else 0.0
+
+
+def ref_daisyworld(final_time=100):
+    lum, white, black, empty = 0.7, 0.40, 0.273, 0.327
+
+    def albedo():
+        return white * 0.75 + black * 0.25 + empty * 0.50
+
+    out = []
+    for _ in range(final_time):
+        ave = (lum * _SOLAR_FLUX * (1 - albedo()) / (4 * _SIGMA)) ** 0.25
+        t_white = (_HEAT * (albedo() - 0.75) + ave ** 4) ** 0.25
+        white = white + white * (_daisy_growth(t_white) * empty - 0.3)
+        t_black = (_HEAT * (albedo() - 0.25) + ave ** 4) ** 0.25   # albedo with the new white area
+        black = black + black * (_daisy_growth(t_black) * empty - 0.3)
+        empty = 1.0 - (black + white)
+        out.append([white, black, empty])
+    return out
+
+
+def _terrame_spline(xs, ys, steps=1000):
+    """TerraME sci.Spline (commit 42b0b3e): sampled Catmull-Rom + linear interpolation."""
+    pts = list(zip(xs, ys))
+    count = len(pts) - 1
+    values = []
+    for i in range(1, count + 1):
+        if i == 1:
+            p0, p1, p2, p3 = pts[0], pts[0], pts[1], pts[2]
+        elif i == count:
+            p0, p1, p2, p3 = pts[-3], pts[-2], pts[-1], pts[-1]
+        else:
+            p0, p1, p2, p3 = pts[i - 2], pts[i - 1], pts[i], pts[i + 1]
+        t = 0.0
+        while t <= 1:
+            def cr(a, b, c, d):
+                return 0.5 * (2 * b + (c - a) * t + (2 * a - 5 * b + 4 * c - d) * t * t
+                              + (3 * b - a - 3 * c + d) * t ** 3)
+            x = cr(p0[0], p1[0], p2[0], p3[0])
+            y = cr(p0[1], p1[1], p2[1], p3[1])
+            if not values or (values[-1][0] != x and values[-1][1] != y):
+                values.append((x, y))
+            t += 1 / steps
+    vx, vy = (np.array(v) for v in zip(*values))
+    return lambda v: float(np.interp(v, vx, vy))
+
+
+def ref_mono_lake(final_time=50):
+    volume = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]
+    surface = _terrame_spline(volume, [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9])
+    elevation = _terrame_spline(volume, [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477])
+    water, out = 2228.0, []
+    for _ in range(final_time):
+        area = surface(water)
+        inflow = area * 0.67 + 150 + 47.6
+        outflow = 3.75 * area + 100 + 33.6
+        water = water + inflow - outflow
+        out.append([water, elevation(water)])
+    return out
+
+
+# ---------------------------------------------------------------------------
+
+KNOWN_DIFFERENCE = {
+    "MonoLake": "TerraME uses a sampled Catmull-Rom spline; the model uses scipy CubicSpline",
+    "Daisyworld": "TerraME recomputes the planet albedo after updating the white area",
+}
+
+CASES = [
+    # name, class, finalTime, attributes, reference, parameters
+    ("Coffee", sd.Coffee, 20, ["temperature"], ref_coffee, {}),
+    ("ChaoticGrowth", sd.ChaoticGrowth, 300, ["pop"], ref_chaotic_growth, {}),
+    ("Homeostasis", sd.Homeostasis, 30, ["stock"], ref_homeostasis, {}),
+    ("LimitedGrowth", sd.LimitedGrowth, 100, ["pop"], ref_limited_growth, {}),
+    ("Lorenz", sd.Lorenz, 10000, ["x", "y", "z"], ref_lorenz, {}),
+    ("PopulationGrowth", sd.PopulationGrowth, 7, ["population"], ref_population_growth, {}),
+    ("PredatorPrey", sd.PredatorPrey, 500, ["prey", "predator"], ref_predator_prey, {}),
+    ("RoomTemperature", sd.RoomTemperature, 24, ["inside", "outside"], ref_room_temperature, {}),
+    ("SIR", sd.SIR, 30, ["susceptible", "infected", "recovered"], ref_sir, {}),
+    ("Tub", sd.Tub, 8, ["water"], ref_tub, {}),
+    ("Yeast", sd.Yeast, 9, ["cells"], ref_yeast, {}),
+    ("Daisyworld", sd.Daisyworld, 100, ["white_area", "black_area", "empty_area"],
+     ref_daisyworld, {}),
+    ("MonoLake", sd.MonoLake, 50, ["water_in_lake", "level"], ref_mono_lake, {}),
+]
+
+
+def _param(case):
+    name = case[0]
+    marks = ([pytest.mark.xfail(reason=KNOWN_DIFFERENCE[name], strict=True)]
+             if name in KNOWN_DIFFERENCE else [])
+    return pytest.param(*case, id=name, marks=marks)
+
+
+@pytest.mark.parametrize("name, cls, final_time, attrs, reference, params",
+                         [_param(c) for c in CASES])
+def test_trajectory_matches_terrame(name, cls, final_time, attrs, reference, params):
+    got = run_dissmodel(cls, final_time, attrs, **params)
+    expected = np.array(reference(final_time), dtype=float)
+    assert got.shape == expected.shape
+    np.testing.assert_allclose(got, expected, rtol=RTOL, atol=ATOL)
+
+
+def test_tub_periodic_inflow():
+    """Exercises the TerraME Event{start = 10, period = 10} (inactive by default)."""
+    got = run_dissmodel(sd.Tub, 30, ["water"], in_flow=20.0)
+    np.testing.assert_allclose(got, np.array(ref_tub(30, in_flow=20.0)), rtol=RTOL, atol=ATOL)