diff --git a/CHANGELOG.md b/CHANGELOG.md
index bad9034..37d5799 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,27 @@
 
 ## [Unreleased]
 
+- **Fix:** `Daisyworld` now follows TerraME: after updating the white area,
+  it computes the planet albedo again for the black daisies. The trajectory
+  changes (up to 0.61 in area); the equilibrium is the same. New attributes
+  `temp_near_white` and `temp_near_black`, as in TerraME.
+- **Fix:** `MonoLake` reads the lake surface and elevation with TerraME's
+  `sci.Spline` (a Catmull-Rom curve sampled at 1000 points per segment and
+  read linearly; new `dissmodel_sysdyn.models._spline.Spline`) instead of
+  scipy's `CubicSpline`, which differed by up to 26 KAF in 50 steps. As in
+  TerraME, a volume outside the curve raises an error instead of
+  extrapolating. scipy, which was used without being a declared
+  dependency, is no longer needed.
+- `PredatorPrey`, `RoomTemperature`, `SIR` and `MonoLake` do their
+  arithmetic in TerraME's order, so that the results are the same doubles.
+- `tests/test_goldens.py` (marker `golden`): every model against TerraME
+  2.0.1 runs from `terrame-goldens` v0.1.0. All 16 deterministic cases
+  (13 models, some with two parameter sets) give **the same doubles as
+  TerraME at every step**; RandomWalk is compared by distribution.
+  `tests/goldens/` holds the goldens, their SHA-256 provenance and the
+  script that copies them. The strict `xfail`s of `Daisyworld` and
+  `MonoLake` in `test_terrame_parity.py` are removed.
+
 - **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
diff --git a/README.md b/README.md
index e802c52..f66aad3 100644
--- a/README.md
+++ b/README.md
@@ -25,6 +25,16 @@ The library includes a wide range of System Dynamics models:
 *   **Environment**: Mono Lake Water Balance.
 *   **Stochasticity**: Random Walk.
 
+## Verification against TerraME
+
+Every model was run in TerraME 2.0.1 with the `sysdyn` package
+([`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0).
+For the 16 deterministic cases (13 models, some with two parameter sets),
+the Python models give **the same doubles as TerraME at every step**, bit for
+bit, for every variable TerraME records; RandomWalk is compared by
+distribution. `tests/test_goldens.py` runs these comparisons
+(`pytest -m golden`, under a second).
+
 ## Examples
 
 The project includes several examples demonstrating how to run the models via CLI, Streamlit, or Jupyter Notebooks.
diff --git a/examples/notebooks/sysdyn_mono_lake.ipynb b/examples/notebooks/sysdyn_mono_lake.ipynb
index 3d66b34..9f137fc 100644
--- a/examples/notebooks/sysdyn_mono_lake.ipynb
+++ b/examples/notebooks/sysdyn_mono_lake.ipynb
@@ -88,17 +88,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "from scipy.interpolate import CubicSpline\n",
     "from dissmodel.core import Model\n",
     "from dissmodel.visualization import track_plot\n",
+    "from dissmodel_sysdyn.models._spline import Spline  # TerraME's sci.Spline\n",
     "\n",
-    "_VOLUME_KAF  = np.array([0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000], dtype=float)\n",
-    "_SURFACE_KM2 = np.array([0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9], dtype=float)\n",
-    "_ELEV_FT     = np.array([6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477], dtype=float)\n",
+    "_VOLUME_KAF  = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]\n",
+    "_SURFACE_KM2 = [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9]\n",
+    "_ELEV_FT     = [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477]\n",
     "\n",
-    "_water_surface   = CubicSpline(_VOLUME_KAF, _SURFACE_KM2, extrapolate=True)\n",
-    "_water_elevation = CubicSpline(_VOLUME_KAF, _ELEV_FT,     extrapolate=True)\n",
+    "# a Catmull-Rom curve sampled at 1000 points per segment, read linearly\n",
+    "_water_surface   = Spline(_VOLUME_KAF, _SURFACE_KM2, steps=1000)\n",
+    "_water_elevation = Spline(_VOLUME_KAF, _ELEV_FT,     steps=1000)\n",
     "\n",
     "@track_plot(\"level\", \"blue\")\n",
     "@track_plot(\"water_in_lake\", \"cyan\")\n",
@@ -109,11 +109,11 @@
     "        self.evap_rate, self.other_out, self.export = evap_rate, other_out, export\n",
     "\n",
     "    def execute(self) -> None:\n",
-    "        surface = float(_water_surface(self.water_in_lake))\n",
+    "        surface = _water_surface.value(self.water_in_lake)\n",
     "        total_in = surface * self.prec_rate + self.runoff + self.other_in\n",
     "        total_out = surface * self.evap_rate + self.export + self.other_out\n",
-    "        self.water_in_lake += total_in - total_out\n",
-    "        self.level = float(_water_elevation(self.water_in_lake))"
+    "        self.water_in_lake = self.water_in_lake + total_in - total_out\n",
+    "        self.level = _water_elevation.value(self.water_in_lake)"
    ]
   },
   {
diff --git a/pyproject.toml b/pyproject.toml
index 808a232..15a646c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -35,3 +35,6 @@ Homepage = "https://github.com/DisSModel/dissmodel-sysdyn"
 
 [tool.pytest.ini_options]
 testpaths = ["tests"]
+markers = [
+    "golden: compares the models with TerraME runs (tests/goldens)",
+]
diff --git a/src/dissmodel_sysdyn/models/_spline.py b/src/dissmodel_sysdyn/models/_spline.py
new file mode 100644
index 0000000..f292b76
--- /dev/null
+++ b/src/dissmodel_sysdyn/models/_spline.py
@@ -0,0 +1,87 @@
+"""
+The ``Spline`` of TerraME's ``sci`` package (commit 42b0b3e, ``lua/Spline.lua``).
+
+``sci.Spline`` is not a cubic spline through the points: it samples a
+Catmull-Rom curve at ``steps`` points per segment (the first and last points
+repeated as end controls) and reads it by linear interpolation between the
+samples. :class:`Spline` reproduces both steps, with the same floating-point
+operations, so that a model that uses it matches TerraME to the last bit.
+
+One consequence of following Lua: the loop over ``t`` adds ``1 / steps``
+to ``t`` and, for some ``steps`` (1000 among them), the sum ends just below
+1, so the last sample stops short of the last point (at 7999.5 instead of
+8000 in MonoLake) and a larger ``x`` raises an error, in TerraME as here.
+"""
+from __future__ import annotations
+
+from bisect import bisect_left
+from collections.abc import Sequence
+
+
+class Spline:
+    """
+    TerraME ``sci.Spline``: a sampled Catmull-Rom curve read linearly.
+
+    Parameters
+    ----------
+    xs, ys : sequence of float
+        Control points, with increasing ``xs``.
+    steps : int, optional
+        Samples per segment, by default 10 (as in ``sci``).
+
+    Raises
+    ------
+    ValueError
+        From :meth:`value`, for an ``x`` outside the range of the points, as
+        TerraME does.
+    """
+
+    def __init__(self, xs: Sequence[float], ys: Sequence[float], steps: int = 10) -> None:
+        points = list(zip((float(x) for x in xs), (float(y) for y in ys), strict=True))
+        self.values = self._build(points, steps)
+        self.xs = [p[0] for p in self.values]
+
+    @staticmethod
+    def _build(points: list[tuple[float, float]], steps: int) -> list[tuple[float, float]]:
+        if len(points) < 3:
+            return points
+        count = len(points) - 1
+        values: list[tuple[float, float]] = []
+        for i in range(1, count + 1):  # 1-based, as in Lua
+            if i == 1:
+                p0, p1, p2, p3 = points[0], points[0], points[1], points[2]
+            elif i == count:
+                p0, p1, p2, p3 = points[-3], points[-2], points[-1], points[-1]
+            else:
+                p0, p1, p2, p3 = points[i - 2], points[i - 1], points[i], points[i + 1]
+            # Lua's ``for t = 0, 1, 1 / steps`` adds the step to t each time
+            t = 0.0
+            dt = 1 / steps
+            while t <= 1:
+                x = 0.5 * ((2 * p1[0]) + (p2[0] - p0[0]) * t
+                           + (2 * p0[0] - 5 * p1[0] + 4 * p2[0] - p3[0]) * t * t
+                           + (3 * p1[0] - p0[0] - 3 * p2[0] + p3[0]) * t * t * t)
+                y = 0.5 * ((2 * p1[1]) + (p2[1] - p0[1]) * t
+                           + (2 * p0[1] - 5 * p1[1] + 4 * p2[1] - p3[1]) * t * t
+                           + (3 * p1[1] - p0[1] - 3 * p2[1] + p3[1]) * t * t * t)
+                if not values or (values[-1][0] != x and values[-1][1] != y):
+                    values.append((x, y))
+                t += dt
+        return values
+
+    def value(self, x: float) -> float:
+        """Linear interpolation between the two samples around ``x``."""
+        xs = self.xs
+        if not xs[0] <= x <= xs[-1]:
+            raise ValueError(f"Spline: {x} is outside the range [{xs[0]}, {xs[-1]}].")
+        i = bisect_left(xs, x)
+        if xs[i] == x:
+            return self.values[i][1]
+        (x0, y0), (x1, y1) = self.values[i - 1], self.values[i]
+        w = (x - x0) / (x1 - x0)
+        return (1 - w) * y0 + w * y1
+
+    __call__ = value
+
+
+__all__ = ["Spline"]
diff --git a/src/dissmodel_sysdyn/models/daisyworld.py b/src/dissmodel_sysdyn/models/daisyworld.py
index bc56784..f80f9af 100644
--- a/src/dissmodel_sysdyn/models/daisyworld.py
+++ b/src/dissmodel_sysdyn/models/daisyworld.py
@@ -65,6 +65,10 @@ class Daisyworld(Model):
 
     Notes
     -----
+    The step follows the TerraME ``sysdyn`` model line by line: the white
+    area is updated first, and the planet albedo used for the black daisies
+    is computed again with the new white area.
+
     The ``@track_plot`` decorators register ``white_area``,
     ``black_area``, ``empty_area``, and ``ave_temp`` for automatic
     live plotting.
@@ -89,6 +93,12 @@ class Daisyworld(Model):
     #: Current average planet temperature (K).
     ave_temp: float
 
+    #: Temperature near the white daisies in the last step (K).
+    temp_near_white: float
+
+    #: Temperature near the black daisies in the last step (K).
+    temp_near_black: float
+
     #: Solar luminosity (fractional).
     sun_luminosity: float
 
@@ -138,15 +148,16 @@ def _planet_albedo(self) -> float:
         )
 
     def execute(self) -> None:
-        pa = self._planet_albedo()
-        self.ave_temp = _planet_temp(self.sun_luminosity, pa)
+        self.ave_temp = _planet_temp(self.sun_luminosity, self._planet_albedo())
 
-        temp_white     = _local_temp(self.ave_temp, pa, self.white_albedo)
-        white_growth   = _daisy_growth_rate(temp_white) * self.empty_area
+        self.temp_near_white = _local_temp(self.ave_temp, self._planet_albedo(), self.white_albedo)
+        white_growth = _daisy_growth_rate(self.temp_near_white) * self.empty_area
         self.white_area += self.white_area * (white_growth - self.decay_rate)
 
-        temp_black     = _local_temp(self.ave_temp, pa, self.black_albedo)
-        black_growth   = _daisy_growth_rate(temp_black) * self.empty_area
+        # As in TerraME, the planet albedo is computed again here, with the
+        # new white area (and the black and empty areas of the start of the step)
+        self.temp_near_black = _local_temp(self.ave_temp, self._planet_albedo(), self.black_albedo)
+        black_growth = _daisy_growth_rate(self.temp_near_black) * self.empty_area
         self.black_area += self.black_area * (black_growth - self.decay_rate)
 
         self.empty_area = self.planet_area - (self.white_area + self.black_area)
diff --git a/src/dissmodel_sysdyn/models/mono_lake.py b/src/dissmodel_sysdyn/models/mono_lake.py
index fddf159..259f86b 100644
--- a/src/dissmodel_sysdyn/models/mono_lake.py
+++ b/src/dissmodel_sysdyn/models/mono_lake.py
@@ -1,17 +1,17 @@
 from __future__ import annotations
 
-from scipy.interpolate import CubicSpline
-import numpy as np
-
 from dissmodel.core import Model
 from dissmodel.visualization import track_plot
 
-_VOLUME_KAF  = np.array([0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000], dtype=float)
-_SURFACE_KM2 = np.array([0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9], dtype=float)
-_ELEV_FT     = np.array([6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477], dtype=float)
+from ._spline import Spline
+
+_VOLUME_KAF  = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]
+_SURFACE_KM2 = [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9]
+_ELEV_FT     = [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477]
 
-_water_surface   = CubicSpline(_VOLUME_KAF, _SURFACE_KM2, extrapolate=True)
-_water_elevation = CubicSpline(_VOLUME_KAF, _ELEV_FT,     extrapolate=True)
+# As in TerraME: sci.Spline{points = ..., steps = 1000}
+_water_surface   = Spline(_VOLUME_KAF, _SURFACE_KM2, steps=1000)
+_water_elevation = Spline(_VOLUME_KAF, _ELEV_FT,     steps=1000)
 
 
 @track_plot("level", "blue")
@@ -43,6 +43,14 @@ class MonoLake(Model):
     export : float, optional
         Annual water exported in KAF/year, by default 100.0.
 
+    Notes
+    -----
+    The lake surface and elevation are read from the volume through
+    TerraME's ``sci.Spline`` (:class:`~dissmodel_sysdyn.models._spline.Spline`):
+    a Catmull-Rom curve sampled at 1000 points per segment and read by
+    linear interpolation, as in the TerraME model. As there, a volume
+    outside 0 to about 7999.5 KAF raises an error instead of extrapolating.
+
     Examples
     --------
     >>> from dissmodel.core import Environment
@@ -96,7 +104,7 @@ def setup(
         self.export        = export
 
     def _surface(self) -> float:
-        return float(_water_surface(self.water_in_lake))
+        return _water_surface.value(self.water_in_lake)
 
     def _total_input(self) -> float:
         return self._surface() * self.prec_rate + self.runoff + self.other_in
@@ -105,8 +113,9 @@ def _total_output(self) -> float:
         return self._surface() * self.evap_rate + self.export + self.other_out
 
     def execute(self) -> None:
-        self.water_in_lake += self._total_input() - self._total_output()
-        self.level = float(_water_elevation(self.water_in_lake))
+        # same order of operations as TerraME: (water + input) - output
+        self.water_in_lake = self.water_in_lake + self._total_input() - self._total_output()
+        self.level = _water_elevation.value(self.water_in_lake)
 
 
 __all__ = ["MonoLake"]
diff --git a/src/dissmodel_sysdyn/models/predatorprey.py b/src/dissmodel_sysdyn/models/predatorprey.py
index 54ce963..c2ae607 100644
--- a/src/dissmodel_sysdyn/models/predatorprey.py
+++ b/src/dissmodel_sysdyn/models/predatorprey.py
@@ -116,12 +116,14 @@ def execute(self) -> None:
         Updates prey and predator populations according to the
         discrete Lotka-Volterra equations.
         """
-        self.prey += (
-            self.prey_growth * self.prey
+        # same order of operations as TerraME, so that the results match it
+        # to the last bit; the predators use the prey already updated
+        self.prey = (
+            self.prey + self.prey_growth * self.prey
             - self.prey_death_pred * self.prey * self.predator
         )
-        self.predator += (
-            -self.pred_death * self.predator
+        self.predator = (
+            self.predator - self.pred_death * self.predator
             + self.pred_growth_kills * self.prey * self.predator
         )
 
diff --git a/src/dissmodel_sysdyn/models/room_temperature.py b/src/dissmodel_sysdyn/models/room_temperature.py
index c9ce46d..cb2928a 100644
--- a/src/dissmodel_sysdyn/models/room_temperature.py
+++ b/src/dissmodel_sysdyn/models/room_temperature.py
@@ -92,7 +92,7 @@ def execute(self) -> None:
         self.outside = climate(t)
         inflow  = self.thermal_inertia * (self.temp_set - self.inside)
         outflow = self.loss_to_outside  * (self.inside   - self.outside)
-        self.inside += inflow - outflow
+        self.inside = self.inside + inflow - outflow  # TerraME's order of operations
 
 
 __all__ = ["RoomTemperature"]
diff --git a/src/dissmodel_sysdyn/models/sir.py b/src/dissmodel_sysdyn/models/sir.py
index 94825e2..ce05f42 100644
--- a/src/dissmodel_sysdyn/models/sir.py
+++ b/src/dissmodel_sysdyn/models/sir.py
@@ -119,15 +119,16 @@ def execute(self) -> None:
         Computes new infections and recoveries, then updates the
         susceptible, infected, and recovered compartments.
         """
-        total = self.susceptible + self.infected + self.recovered
-        alpha = self.contacts * self.probability
+        # same order of operations as TerraME, so that the results match it
+        # to the last bit
+        proportion = self.susceptible / (self.susceptible + self.infected + self.recovered)
 
-        new_infected  = self.infected * alpha * (self.susceptible / total)
+        new_infected  = self.infected * self.contacts * self.probability * proportion
         new_recovered = self.infected / self.duration
 
-        self.susceptible -= new_infected
-        self.infected    += new_infected - new_recovered
-        self.recovered   += new_recovered
+        self.susceptible = self.susceptible - new_infected
+        self.recovered   = self.recovered + new_recovered
+        self.infected    = self.infected + new_infected - new_recovered
 
 
 __all__ = ["SIR"]
diff --git a/tests/goldens/README.md b/tests/goldens/README.md
new file mode 100644
index 0000000..64b6219
--- /dev/null
+++ b/tests/goldens/README.md
@@ -0,0 +1,26 @@
+# TerraME reference values for the sysdyn models
+
+The `sysdyn` suite of
+[`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0:
+what TerraME 2.0.1 computed running the models of the
+[`sysdyn`](https://github.com/TerraME/sysdyn) package (commit `9e671cb`; `sci`
+commit `42b0b3e` for MonoLake's spline).
+
+- `<case>.csv.gz`: a column `step` and one column per variable, from step 0
+  (the state after `init`) to `finalTime`, written with 17 significant digits,
+  so every double is recovered exactly. A variable that does not exist yet
+  at step 0 is empty. Lorenz keeps steps 0 to 100 and then every 100th step.
+- `random_walk.csv.gz`: `seed, step, value` for seeds 1 to 50.
+- `cases.csv`: model, arguments, final time and variables of each case.
+- `reference.json`: the source version, the package commits and the SHA-256
+  of each source file.
+
+To copy them again from a copy of `terrame-goldens`:
+
+```bash
+python tests/goldens/make_reference.py /path/to/terrame-goldens
+```
+
+The values are the output of TerraME and of the `sysdyn` and `sci` packages
+(INPE and TerraLAB/UFOP); they are offered here under CC-BY-4.0, as the
+goldens they come from.
diff --git a/tests/goldens/cases.csv b/tests/goldens/cases.csv
new file mode 100644
index 0000000..5f100f5
--- /dev/null
+++ b/tests/goldens/cases.csv
@@ -0,0 +1,18 @@
+case,model,arguments,seeds,final_time,variables
+chaotic_growth,ChaoticGrowth,,,300,pop
+coffee,Coffee,,,20,temperature
+daisyworld,Daisyworld,,,100,whiteArea;blackArea;emptyArea;aveTemp;tempNearWhite;tempNearBlack
+homeostasis,Homeostasis,,,30,stock
+homeostasis_100,Homeostasis,stock=100.0,,30,stock
+limited_growth,LimitedGrowth,,,100,pop
+lorenz,Lorenz,,,10000,x;y;z
+mono_lake,MonoLake,,,50,waterInLake;level
+mono_lake_export_50,MonoLake,export=50,,50,waterInLake;level
+population_growth,PopulationGrowth,,,7,population;growth
+predator_prey,PredatorPrey,,,500,prey;predator
+room_temperature,RoomTemperature,,,24,inside;outside
+sir,SIR,,,30,susceptible;infected;recovered;contacts
+tub,Tub,,,8,water
+yeast,Yeast,,,9,cells
+yeast_rate_2,Yeast,rate=2.0,,9,cells
+random_walk,RandomWalk,,50,100,value
diff --git a/tests/goldens/chaotic_growth.csv.gz b/tests/goldens/chaotic_growth.csv.gz
new file mode 100644
index 0000000..2f4fb3b
Binary files /dev/null and b/tests/goldens/chaotic_growth.csv.gz differ
diff --git a/tests/goldens/coffee.csv.gz b/tests/goldens/coffee.csv.gz
new file mode 100644
index 0000000..ebf3080
Binary files /dev/null and b/tests/goldens/coffee.csv.gz differ
diff --git a/tests/goldens/daisyworld.csv.gz b/tests/goldens/daisyworld.csv.gz
new file mode 100644
index 0000000..f9a7383
Binary files /dev/null and b/tests/goldens/daisyworld.csv.gz differ
diff --git a/tests/goldens/homeostasis.csv.gz b/tests/goldens/homeostasis.csv.gz
new file mode 100644
index 0000000..00f5b5d
Binary files /dev/null and b/tests/goldens/homeostasis.csv.gz differ
diff --git a/tests/goldens/homeostasis_100.csv.gz b/tests/goldens/homeostasis_100.csv.gz
new file mode 100644
index 0000000..968bd67
Binary files /dev/null and b/tests/goldens/homeostasis_100.csv.gz differ
diff --git a/tests/goldens/limited_growth.csv.gz b/tests/goldens/limited_growth.csv.gz
new file mode 100644
index 0000000..c2d721b
Binary files /dev/null and b/tests/goldens/limited_growth.csv.gz differ
diff --git a/tests/goldens/lorenz.csv.gz b/tests/goldens/lorenz.csv.gz
new file mode 100644
index 0000000..5a5d5fc
Binary files /dev/null and b/tests/goldens/lorenz.csv.gz differ
diff --git a/tests/goldens/make_reference.py b/tests/goldens/make_reference.py
new file mode 100644
index 0000000..aee2d53
--- /dev/null
+++ b/tests/goldens/make_reference.py
@@ -0,0 +1,60 @@
+"""
+Copy the reference values used by ``tests/test_goldens.py`` from the TerraME
+``sysdyn`` goldens of ``terrame-goldens``.
+
+Every golden is copied as is (``<case>.csv.gz``: a column ``step`` and one per
+variable, values with 17 significant digits), except Lorenz (10000 steps),
+which keeps steps 0 to 100 and then every 100th step. ``reference.json``
+records the source version, the package commits and the SHA-256 of each
+source file.
+
+    python tests/goldens/make_reference.py /path/to/terrame-goldens
+"""
+from __future__ import annotations
+
+import gzip
+import hashlib
+import io
+import json
+import sys
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+
+
+def write_gz(text: str, dest: Path) -> None:
+    buf = io.BytesIO()
+    with gzip.GzipFile(filename="", mode="wb", fileobj=buf, mtime=0, compresslevel=9) as gz:
+        gz.write(text.encode())
+    dest.write_bytes(buf.getvalue())
+
+
+def main(root: Path) -> None:
+    src = root / "goldens" / "sysdyn"
+    sources = {}
+    for path in sorted(src.glob("*.csv.gz")):
+        sources[path.relative_to(root).as_posix()] = hashlib.sha256(path.read_bytes()).hexdigest()
+        with gzip.open(path, "rt", newline="") as f:
+            lines = f.read().splitlines()
+        if path.name == "lorenz.csv.gz":
+            lines = [lines[0]] + [
+                ln for ln in lines[1:]
+                if (s := int(ln.split(",", 1)[0])) <= 100 or s % 100 == 0
+            ]
+        write_gz("\n".join(lines) + "\n", HERE / path.name)
+    (HERE / "cases.csv").write_bytes((src / "cases.csv").read_bytes())
+    manifest = json.loads((root / "manifest.json").read_text())
+    meta = {
+        "source": "https://github.com/LambdaGeo/terrame-goldens",
+        "version": manifest["version"],
+        "terrame": manifest["terrame"],
+        "packages": [p for p in manifest["packages"] if p["name"] in ("sysdyn", "sci")],
+        "files": sources,
+        "lorenz_steps": "0 to 100, then every 100th step",
+    }
+    (HERE / "reference.json").write_text(json.dumps(meta, indent=2) + "\n")
+    print(f"{len(sources)} goldens copied")
+
+
+if __name__ == "__main__":
+    main(Path(sys.argv[1]))
diff --git a/tests/goldens/mono_lake.csv.gz b/tests/goldens/mono_lake.csv.gz
new file mode 100644
index 0000000..ed6a7eb
Binary files /dev/null and b/tests/goldens/mono_lake.csv.gz differ
diff --git a/tests/goldens/mono_lake_export_50.csv.gz b/tests/goldens/mono_lake_export_50.csv.gz
new file mode 100644
index 0000000..0e3a23a
Binary files /dev/null and b/tests/goldens/mono_lake_export_50.csv.gz differ
diff --git a/tests/goldens/population_growth.csv.gz b/tests/goldens/population_growth.csv.gz
new file mode 100644
index 0000000..94d1e1b
Binary files /dev/null and b/tests/goldens/population_growth.csv.gz differ
diff --git a/tests/goldens/predator_prey.csv.gz b/tests/goldens/predator_prey.csv.gz
new file mode 100644
index 0000000..fe25cc1
Binary files /dev/null and b/tests/goldens/predator_prey.csv.gz differ
diff --git a/tests/goldens/random_walk.csv.gz b/tests/goldens/random_walk.csv.gz
new file mode 100644
index 0000000..10f9931
Binary files /dev/null and b/tests/goldens/random_walk.csv.gz differ
diff --git a/tests/goldens/reference.json b/tests/goldens/reference.json
new file mode 100644
index 0000000..a7cbaf6
--- /dev/null
+++ b/tests/goldens/reference.json
@@ -0,0 +1,37 @@
+{
+  "source": "https://github.com/LambdaGeo/terrame-goldens",
+  "version": "v0.1.0",
+  "terrame": "2.0.1",
+  "packages": [
+    {
+      "name": "sysdyn",
+      "repository": "https://github.com/TerraME/sysdyn",
+      "commit": "9e671cb86ee33143bd910f1243bc71f2d1c5e5a9"
+    },
+    {
+      "name": "sci",
+      "repository": "https://github.com/TerraME/sci",
+      "commit": "42b0b3e319b6b6fe7bc2d3566e23876db12a2bac"
+    }
+  ],
+  "files": {
+    "goldens/sysdyn/chaotic_growth.csv.gz": "6328b3248918b88ec4adfc1002a99ec6c4d252fa1e6a3beaa0c1477faa528c75",
+    "goldens/sysdyn/coffee.csv.gz": "cf32769b784339a21489a54ef05f1e59c403da6857fa7421dd28e05f1581a4c1",
+    "goldens/sysdyn/daisyworld.csv.gz": "da394371bf4dc84114ae5c5815797d20aaa7fd85b6e7b182446e180f7d747c6c",
+    "goldens/sysdyn/homeostasis.csv.gz": "2cb49a621d7153e98632f7f0e11bfdc47ea04a105a2f659d6976bfa3dfc0fd4f",
+    "goldens/sysdyn/homeostasis_100.csv.gz": "6896a8faed843dc0c42f58194b284381fad971e793bed8de531ef26282ac280b",
+    "goldens/sysdyn/limited_growth.csv.gz": "9972047016a135b12c2e70e26209d438e1b911dd44ac927c47f4b8eab5f28b87",
+    "goldens/sysdyn/lorenz.csv.gz": "96c83fea2d715e3f649be194c7c5a2f436b2e11d9784bf38d86cb09e915038f6",
+    "goldens/sysdyn/mono_lake.csv.gz": "0b085b0286bb7af814f51c549c0fed539e142b2b40b19a37e74d166ca8f7db80",
+    "goldens/sysdyn/mono_lake_export_50.csv.gz": "fdc1b31b6181fc4a5d4ea6b82525d521f1ab0347e1162853d31423c7a62f6524",
+    "goldens/sysdyn/population_growth.csv.gz": "ac1e0a04495db4ab8239be3df64261420fdb6714e1f824bda7c4c4a094d87736",
+    "goldens/sysdyn/predator_prey.csv.gz": "51feea1fb7383c1eea15ee7cf8579003fe526b308644faf4423338aa6380434b",
+    "goldens/sysdyn/random_walk.csv.gz": "851545dee78e65f5f7a7d1b22fe2fbd4fffede22ceef795ce7a1536a891c7a6b",
+    "goldens/sysdyn/room_temperature.csv.gz": "353331f5f08929a2c9cb707bb83644b5ae1a114dbbffb75cfeb654bcc09c3276",
+    "goldens/sysdyn/sir.csv.gz": "c3381b461ba01945be76f257bfadbd4416ae532704663580cabd992d39c331ae",
+    "goldens/sysdyn/tub.csv.gz": "2a47c516a64fca621f1c86d9ac45c9c3c4b34a24aa9fdd2098288a022810bee8",
+    "goldens/sysdyn/yeast.csv.gz": "5db25a6ed8cf4ee8d5ff0bd171543d1e32c11f2bf26f4ae596ee91bd9dc41eba",
+    "goldens/sysdyn/yeast_rate_2.csv.gz": "7052dec647e5ff9d0f84d9e70d9549fe70349e1a04f11c07697f4e58a41e9b27"
+  },
+  "lorenz_steps": "0 to 100, then every 100th step"
+}
diff --git a/tests/goldens/room_temperature.csv.gz b/tests/goldens/room_temperature.csv.gz
new file mode 100644
index 0000000..52cf01d
Binary files /dev/null and b/tests/goldens/room_temperature.csv.gz differ
diff --git a/tests/goldens/sir.csv.gz b/tests/goldens/sir.csv.gz
new file mode 100644
index 0000000..aa25677
Binary files /dev/null and b/tests/goldens/sir.csv.gz differ
diff --git a/tests/goldens/tub.csv.gz b/tests/goldens/tub.csv.gz
new file mode 100644
index 0000000..24eb898
Binary files /dev/null and b/tests/goldens/tub.csv.gz differ
diff --git a/tests/goldens/yeast.csv.gz b/tests/goldens/yeast.csv.gz
new file mode 100644
index 0000000..a9c725a
Binary files /dev/null and b/tests/goldens/yeast.csv.gz differ
diff --git a/tests/goldens/yeast_rate_2.csv.gz b/tests/goldens/yeast_rate_2.csv.gz
new file mode 100644
index 0000000..b018e74
Binary files /dev/null and b/tests/goldens/yeast_rate_2.csv.gz differ
diff --git a/tests/test_goldens.py b/tests/test_goldens.py
new file mode 100644
index 0000000..376950c
--- /dev/null
+++ b/tests/test_goldens.py
@@ -0,0 +1,96 @@
+"""
+The models against TerraME itself.
+
+``goldens/`` holds the ``sysdyn`` suite of ``terrame-goldens`` (see
+``goldens/README.md``): what TerraME 2.0.1 computed running the models of the
+``sysdyn`` package, written with 17 significant digits so that every double
+is recovered exactly. Each model runs with the same parameters and must give
+**the same doubles** as TerraME at every step (bit for bit), for every
+variable TerraME records.
+
+RandomWalk draws random numbers, so it is compared by distribution with the
+50 TerraME runs.
+"""
+from __future__ import annotations
+
+import csv
+import gzip
+import re
+from pathlib import Path
+
+import numpy as np
+import pytest
+from dissmodel.core import Environment, Model
+
+import dissmodel_sysdyn.models as sd
+
+pytestmark = pytest.mark.golden
+
+GOLDENS = Path(__file__).resolve().parent / "goldens"
+
+with open(GOLDENS / "cases.csv", newline="") as _f:
+    CASES = [c for c in csv.DictReader(_f) if not c["seeds"]]
+
+
+def snake(name: str) -> str:
+    """TerraME's camelCase attribute name in DisSModel's snake_case."""
+    return re.sub(r"(?<!^)(?=[A-Z])", "_", name).lower()
+
+
+def parse_arguments(text: str) -> dict[str, float]:
+    return {snake(k): float(v) for k, v in (a.split("=") for a in text.split(";") if a)}
+
+
+def load(case: str) -> tuple[list[str], dict[int, list[float]]]:
+    with gzip.open(GOLDENS / f"{case}.csv.gz", "rt", newline="") as f:
+        rows = list(csv.reader(f))
+    nan = float("nan")
+    return rows[0][1:], {int(r[0]): [float(v) if v else nan for v in r[1:]] for r in rows[1:]}
+
+
+class _Recorder(Model):
+    """Reads the model after each step (created after the model, so it runs after it)."""
+
+    def setup(self, model: Model, attrs: list[str]) -> None:  # type: ignore[override]
+        self.model, self.attrs, self.rows = model, attrs, {}
+        self.rows[0] = self.read()
+
+    def read(self) -> list[float]:
+        return [float(getattr(self.model, a, float("nan"))) for a in self.attrs]
+
+    def execute(self) -> None:
+        self.rows[int(self.env.now())] = self.read()
+
+
+@pytest.mark.parametrize("case", CASES, ids=[c["case"] for c in CASES])
+def test_trajectory_is_bit_identical_to_terrame(case: dict[str, str]) -> None:
+    variables, ref = load(case["case"])
+    attrs = [snake(v) for v in variables]
+    env = Environment(start_time=1, end_time=int(case["final_time"]))
+    model = getattr(sd, case["model"])(**parse_arguments(case["arguments"]))
+    rec = _Recorder(model=model, attrs=attrs)
+    env.run()
+    for t, expected in ref.items():
+        got = rec.rows[t]
+        for var, e, g in zip(variables, expected, got, strict=True):
+            if np.isnan(e):  # not defined yet in TerraME (Daisyworld at step 0)
+                continue
+            assert g == e, f"{case['case']}.{var} at step {t}: {g!r} != TerraME {e!r}"
+
+
+def test_random_walk_matches_terrame_distribution() -> None:
+    with gzip.open(GOLDENS / "random_walk.csv.gz", "rt", newline="") as f:
+        rows = list(csv.DictReader(f))
+    final_tm = np.array([int(r["value"]) for r in rows if r["step"] == "100"])
+    ours = []
+    for seed in range(1, 201):
+        env = Environment(start_time=1, end_time=100)
+        model = sd.RandomWalk(seed=seed)
+        env.run()
+        ours.append(model.value)
+    ours_a = np.array(ours)
+    # a walk of 100 steps of +-1: mean 0, standard deviation 10
+    se = np.sqrt(final_tm.var(ddof=1) / len(final_tm) + ours_a.var(ddof=1) / len(ours_a))
+    assert abs(ours_a.mean() - final_tm.mean()) < 4 * se
+    assert 7 < ours_a.std() < 13 and 7 < final_tm.std() < 13
+    assert set(np.unique(ours_a % 2)) == {0}  # even after an even number of steps
diff --git a/tests/test_spline.py b/tests/test_spline.py
new file mode 100644
index 0000000..82bda16
--- /dev/null
+++ b/tests/test_spline.py
@@ -0,0 +1,46 @@
+"""The ``sci.Spline`` of TerraME used by MonoLake."""
+from __future__ import annotations
+
+import pytest
+
+from dissmodel_sysdyn.models._spline import Spline
+
+XS = [0, 1000, 2000, 3000]
+YS = [0, 24.7, 35.3, 48.6]
+
+
+def test_passes_through_the_points():
+    s = Spline(XS, YS, steps=1000)
+    for x, y in zip(XS[:-1], YS[:-1], strict=True):
+        assert s.value(x) == pytest.approx(y, abs=1e-12)
+
+
+def test_last_sample_falls_short_of_the_last_point_as_in_lua():
+    # Lua's ``for t = 0, 1, 1 / steps`` adds 1/1000 a thousand times and ends
+    # just below 1, so the curve stops before the last point; TerraME then
+    # refuses x beyond the last sample
+    s = Spline(XS, YS, steps=1000)
+    assert 2999 < s.xs[-1] < 3000
+    with pytest.raises(ValueError):
+        s.value(3000)
+    assert Spline(XS, YS, steps=10).xs[-1] == 3000
+
+
+def test_linear_between_samples():
+    s = Spline(XS, YS, steps=10)
+    (x0, y0), (x1, y1) = s.values[3], s.values[4]
+    mid = (x0 + x1) / 2
+    assert s.value(mid) == pytest.approx((y0 + y1) / 2)
+
+
+def test_outside_the_range_raises_as_in_terrame():
+    s = Spline(XS, YS)
+    with pytest.raises(ValueError):
+        s.value(-1)
+    with pytest.raises(ValueError):
+        s.value(3000.5)
+
+
+def test_two_points_are_kept_as_given():
+    s = Spline([0, 10], [1, 3])
+    assert s.value(5) == 2
diff --git a/tests/test_terrame_parity.py b/tests/test_terrame_parity.py
index eeff55d..fd6be5e 100644
--- a/tests/test_terrame_parity.py
+++ b/tests/test_terrame_parity.py
@@ -15,16 +15,7 @@
 ``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.
+``test_goldens.py`` compares the models with the output of TerraME itself.
 """
 from __future__ import annotations
 
@@ -241,11 +232,6 @@ def ref_mono_lake(final_time=50):
 
 # ---------------------------------------------------------------------------
 
-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, {}),
@@ -266,10 +252,7 @@ def ref_mono_lake(final_time=50):
 
 
 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)
+    return pytest.param(*case, id=case[0])
 
 
 @pytest.mark.parametrize("name, cls, final_time, attrs, reference, params",