diff --git a/raster/r.sim/r.sim.sediment/main.c b/raster/r.sim/r.sim.sediment/main.c
index 6734d0cf381..3795af88541 100644
--- a/raster/r.sim/r.sim.sediment/main.c
+++ b/raster/r.sim/r.sim.sediment/main.c
@@ -57,6 +57,7 @@
 /* INCLUDES                     */
 
 /********************************/
+#include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -102,7 +103,6 @@ int main(int argc, char *argv[])
     struct Cell_head cellhd;
     struct options parm;
     struct flags flag;
-    long seed_value;
 
     G_gisinit(argv[0]);
 
@@ -279,21 +279,16 @@ int main(int argc, char *argv[])
     parm.diffc->description = _("Water diffusion constant");
     parm.diffc->guisection = _("Parameters");
 
-    parm.seed = G_define_option();
+    parm.seed = G_define_standard_option(G_OPT_M_SEED);
     parm.seed->key = "random_seed";
-    parm.seed->type = TYPE_INTEGER;
-    parm.seed->required = NO;
-    parm.seed->label = _("Seed for random number generator");
-    parm.seed->description =
-        _("The same seed can be used to obtain same results"
-          " or random seed can be generated by other means.");
 
     flag.generateSeed = G_define_flag();
     flag.generateSeed->key = 's';
-    flag.generateSeed->label = _("Generate random seed");
+    flag.generateSeed->label =
+        _("Generate random seed (result is non-deterministic) [deprecated]");
     flag.generateSeed->description =
-        _("Automatically generates random seed for random number"
-          " generator (use when you don't want to provide the seed option)");
+        _("This flag is deprecated and will be removed in a future release. "
+          "Seeding is automatic or use parameter random_seed.");
 
     parm.threads = G_define_option();
     parm.threads->key = "nprocs";
@@ -329,20 +324,7 @@ int main(int argc, char *argv[])
                              ? SUMMARY_JSON
                              : SUMMARY_PLAIN;
 
-    if (flag.generateSeed->answer) {
-        seed_value = G_srand48_auto();
-        G_verbose_message(_("Generated random seed (-s): %ld"), seed_value);
-    }
-    else if (parm.seed->answer) {
-        seed_value = atol(parm.seed->answer);
-        G_srand48(seed_value);
-        G_verbose_message(_("Read random seed from %s option: %ld"),
-                          parm.seed->key, seed_value);
-    }
-    else {
-        /* default as it used to be */
-        G_srand48(12345);
-    }
+    int64_t seed = simwe_seed(parm.seed, flag.generateSeed);
 
     G_get_set_window(&cellhd);
 
@@ -353,6 +335,7 @@ int main(int argc, char *argv[])
     ObservationPoints points = {0};
     settings.hhmax = settings.halpha = settings.hbeta = 0;
     settings.ts = false;
+    settings.seed = seed;
     Inputs inputs = {0};
     Outputs outputs = {0};
     Grids grids = {0};
diff --git a/raster/r.sim/r.sim.sediment/r.sim.sediment.html b/raster/r.sim/r.sim.sediment/r.sim.sediment.html
index a9ab98b89b6..acb90dcedd3 100644
--- a/raster/r.sim/r.sim.sediment/r.sim.sediment.html
+++ b/raster/r.sim/r.sim.sediment/r.sim.sediment.html
@@ -44,8 +44,9 @@ <h2>DESCRIPTION</h2>
 <br>
 
 <p>
-Increasing the number of threads with <b>nprocs</b> does not really
-speed up the simulation.
+Increasing the number of threads with <b>nprocs</b> speeds up the
+simulation. The random numbers, the seed and the reproducibility of the
+results with more than one thread are as in <em>r.sim.water</em>.
 
 <h2>NOTES</h2>
 
@@ -93,6 +94,7 @@ <h3>Run summary</h3>
     "walkers_requested": 60,
     "walkers_generated": 78,
     "walkers_remaining": 43,
+    "seed": 1,
     "duration": 60,
     "simulated_time": 54.891343113611583,
     "time_step": 5.4891343113611581,
diff --git a/raster/r.sim/r.sim.sediment/r.sim.sediment.md b/raster/r.sim/r.sim.sediment/r.sim.sediment.md
index 5d665094f64..87e4afe0a6a 100644
--- a/raster/r.sim/r.sim.sediment/r.sim.sediment.md
+++ b/raster/r.sim/r.sim.sediment/r.sim.sediment.md
@@ -37,8 +37,9 @@ controlled by *duration* \[minutes\] parameter. If the resulting
 erosion/deposition map is noisy, higher number of walkers, given by
 *nwalkers* should be used.  
 
-Increasing the number of threads with **nprocs** does not really speed
-up the simulation.
+Increasing the number of threads with **nprocs** speeds up the
+simulation. The random numbers, the seed and the reproducibility of the
+results with more than one thread are as in *r.sim.water*.
 
 ## NOTES
 
@@ -83,6 +84,7 @@ r.sim.sediment elevation=elevation water_depth=water_depth detachment_coeff=deta
     "walkers_requested": 60,
     "walkers_generated": 78,
     "walkers_remaining": 43,
+    "seed": 1,
     "duration": 60,
     "simulated_time": 54.891343113611583,
     "time_step": 5.4891343113611581,
diff --git a/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_random_test.py b/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_random_test.py
new file mode 100644
index 00000000000..21c0f2c1379
--- /dev/null
+++ b/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_random_test.py
@@ -0,0 +1,117 @@
+# SPDX-FileCopyrightText: 2026 GRASS Development Team
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+"""Tests of the random numbers of r.sim.sediment.
+
+Each walker draws from a random number state of its own, set by the seed and
+the walker's number, so its random numbers do not depend on the number of
+threads. The sediment flux still does at nprocs > 1, because walkers in the
+same cell add to the concentration grid without synchronization, so the test
+which compares numbers of threads compares the walker positions instead. A
+walker's path depends only on the slope and its random numbers.
+"""
+
+import os
+from io import StringIO
+
+import numpy as np
+import pytest
+
+import grass.script as gs
+from grass.exceptions import CalledModuleError
+from grass.tools import Tools
+
+
+@pytest.fixture(scope="module")
+def session(tmp_path_factory):
+    """Session in an XY project with a bowl as elevation and constant inputs."""
+    project = tmp_path_factory.mktemp("r_sim_sediment_random") / "xy_test"
+    gs.create_project(project)
+    with (
+        gs.setup.init(project, env=os.environ.copy()) as session,
+        Tools(session=session) as tools,
+    ):
+        tools.g_region(s=0, n=40, w=0, e=40, res=1)
+        tools.r_mapcalc(
+            expression="elevation = 0.01 * ((col() - 20.5)^2 + (row() - 20.5)^2)"
+        )
+        tools.r_mapcalc(expression="water_depth = 0.1")
+        tools.r_mapcalc(expression="detachment = 0.001")
+        tools.r_mapcalc(expression="transport = 0.001")
+        tools.r_mapcalc(expression="shear_stress = 0")
+        yield session
+
+
+def simulate(session, walkers=None, env=None, **kwargs):
+    """Run the simulation, return the flux and the walker coordinates."""
+    tools = Tools(env=env) if env else Tools(session=session)
+    if walkers:
+        kwargs["walkers_output"] = walkers
+    flux = tools.r_sim_sediment(
+        elevation="elevation",
+        water_depth="water_depth",
+        detachment_coeff="detachment",
+        transport_coeff="transport",
+        shear_stress="shear_stress",
+        sediment_flux=np.array,
+        duration=2,
+        **kwargs,
+    )
+    if not walkers:
+        return flux, None
+    text = tools.v_out_ascii(input=walkers, format="point", precision=15).text
+    return flux, np.loadtxt(StringIO(text), delimiter="|")
+
+
+def test_same_seed_gives_same_flux(session):
+    """Two runs with the same seed give the same sediment flux."""
+    first, _ = simulate(session, random_seed=1)
+    second, _ = simulate(session, random_seed=1)
+    assert np.array_equal(first, second)
+
+
+def test_different_seeds_give_different_results(session):
+    """Runs with different seeds give different fluxes and walker positions."""
+    first_flux, first_walkers = simulate(session, "walkers_seed_1", random_seed=1)
+    second_flux, second_walkers = simulate(session, "walkers_seed_2", random_seed=2)
+    assert not np.array_equal(first_flux, second_flux)
+    assert not np.array_equal(first_walkers, second_walkers)
+
+
+def test_generated_seed_is_the_seed_option(session):
+    """Without a seed, GRASS_RANDOM_SEED gives what the seed option gives."""
+    env = session.env.copy()
+    env["GRASS_RANDOM_SEED"] = "3"
+    generated, _ = simulate(session, env=env)
+    given, _ = simulate(session, random_seed=3)
+    assert np.array_equal(generated, given)
+
+
+@pytest.mark.parametrize("nprocs", [2, 4])
+def test_results_do_not_depend_on_nprocs(session, nprocs):
+    """The flux and the walkers are identical with any number of threads.
+
+    The water depth is an input here, so nothing a walker reacts to is
+    changed by the other walkers of the same step.
+    """
+    serial_flux, serial = simulate(session, f"walkers_{nprocs}_serial", random_seed=5)
+    parallel_flux, parallel = simulate(
+        session, f"walkers_{nprocs}_parallel", random_seed=5, nprocs=nprocs
+    )
+    assert serial.size
+    assert np.array_equal(parallel, serial)
+    assert np.array_equal(parallel_flux, serial_flux)
+
+
+@pytest.mark.parametrize("seed", [-(2**31) - 1, 2**32])
+def test_seed_outside_range_is_an_error(session, seed):
+    """A seed the generator cannot use is refused, not silently wrapped."""
+    with pytest.raises(CalledModuleError, match="outside the range"):
+        simulate(session, random_seed=seed)
+
+
+@pytest.mark.parametrize("seed", ["12abc", "1.5", "-", "99999999999999999999"])
+def test_seed_which_is_not_an_integer_is_an_error(session, seed):
+    """A seed which the parser lets through but which is not an integer is refused."""
+    with pytest.raises(CalledModuleError, match="Invalid random seed"):
+        simulate(session, random_seed=seed)
diff --git a/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_summary_test.py b/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_summary_test.py
index 65b1c9b98bb..4fc1ee25067 100644
--- a/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_summary_test.py
+++ b/raster/r.sim/r.sim.sediment/tests/r_sim_sediment_summary_test.py
@@ -15,6 +15,7 @@
     "walkers_requested",
     "walkers_generated",
     "walkers_remaining",
+    "seed",
     "duration",
     "simulated_time",
     "time_step",
diff --git a/raster/r.sim/r.sim.sediment/testsuite/data/reference_erdep.pack b/raster/r.sim/r.sim.sediment/testsuite/data/reference_erdep.pack
index dde1747fe5d..9971824af2b 100644
Binary files a/raster/r.sim/r.sim.sediment/testsuite/data/reference_erdep.pack and b/raster/r.sim/r.sim.sediment/testsuite/data/reference_erdep.pack differ
diff --git a/raster/r.sim/r.sim.sediment/testsuite/data/reference_sedflux.pack b/raster/r.sim/r.sim.sediment/testsuite/data/reference_sedflux.pack
index 027ed194698..87c9ccba74c 100644
Binary files a/raster/r.sim/r.sim.sediment/testsuite/data/reference_sedflux.pack and b/raster/r.sim/r.sim.sediment/testsuite/data/reference_sedflux.pack differ
diff --git a/raster/r.sim/r.sim.water/main.c b/raster/r.sim/r.sim.water/main.c
index 6a7e5bef98c..e1055c16cf7 100644
--- a/raster/r.sim/r.sim.water/main.c
+++ b/raster/r.sim/r.sim.water/main.c
@@ -64,6 +64,7 @@
 /* INCLUDES                     */
 
 /********************************/
+#include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -102,7 +103,6 @@ int main(int argc, char *argv[])
     struct Cell_head cellhd;
     struct options parm;
     struct flags flag;
-    long seed_value;
 
     G_gisinit(argv[0]);
 
@@ -314,21 +314,16 @@ int main(int argc, char *argv[])
     flag.tserie->description = _("Time-series output");
     flag.tserie->guisection = _("Output");
 
-    parm.seed = G_define_option();
+    parm.seed = G_define_standard_option(G_OPT_M_SEED);
     parm.seed->key = "random_seed";
-    parm.seed->type = TYPE_INTEGER;
-    parm.seed->required = NO;
-    parm.seed->label = _("Seed for random number generator");
-    parm.seed->description =
-        _("The same seed can be used to obtain same results"
-          " or random seed can be generated by other means.");
 
     flag.generateSeed = G_define_flag();
     flag.generateSeed->key = 's';
-    flag.generateSeed->label = _("Generate random seed");
+    flag.generateSeed->label =
+        _("Generate random seed (result is non-deterministic) [deprecated]");
     flag.generateSeed->description =
-        _("Automatically generates random seed for random number"
-          " generator (use when you don't want to provide the seed option)");
+        _("This flag is deprecated and will be removed in a future release. "
+          "Seeding is automatic or use parameter random_seed.");
 
     parm.threads = G_define_option();
     parm.threads->key = "nprocs";
@@ -365,20 +360,7 @@ int main(int argc, char *argv[])
                              ? SUMMARY_JSON
                              : SUMMARY_PLAIN;
 
-    if (flag.generateSeed->answer) {
-        seed_value = G_srand48_auto();
-        G_verbose_message(_("Generated random seed (-s): %ld"), seed_value);
-    }
-    else if (parm.seed->answer) {
-        seed_value = atol(parm.seed->answer);
-        G_srand48(seed_value);
-        G_verbose_message(_("Read random seed from %s option: %ld"),
-                          parm.seed->key, seed_value);
-    }
-    else {
-        /* default as it used to be */
-        G_srand48(12345);
-    }
+    int64_t seed = simwe_seed(parm.seed, flag.generateSeed);
 
     G_get_set_window(&cellhd);
 
@@ -417,6 +399,7 @@ int main(int argc, char *argv[])
     G_debug(3, "xmax: %f, ymax: %f", geometry.xmax, geometry.ymax);
 
     settings.ts = flag.tserie->answer;
+    settings.seed = seed;
 
     inputs.elevin = parm.elevin->answer;
     inputs.dxin = parm.dxin->answer;
diff --git a/raster/r.sim/r.sim.water/r.sim.water.html b/raster/r.sim/r.sim.water/r.sim.water.html
index 7e125f78243..b9f2fbc0ab1 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.html
+++ b/raster/r.sim/r.sim.water/r.sim.water.html
@@ -206,6 +206,7 @@ <h3>Run summary</h3>
  <tr><td><code>walkers_requested</code></td><td>Number of walkers from <b>nwalkers</b>, by default twice the number of cells</td><td>count</td></tr>
  <tr><td><code>walkers_generated</code></td><td>Walkers created, at least one per cell and more where the source rate is higher</td><td>count</td></tr>
  <tr><td><code>walkers_remaining</code></td><td>Walkers still in the domain at the end of the run</td><td>count</td></tr>
+ <tr><td><code>seed</code></td><td>Seed of the random numbers, given or generated</td><td></td></tr>
  <tr><td><code>duration</code></td><td>Requested simulation length (<b>duration</b>)</td><td>s</td></tr>
  <tr><td><code>simulated_time</code></td><td>Simulated time reached at the end of the run</td><td>s</td></tr>
  <tr><td><code>time_step</code></td><td>Simulated time per iteration</td><td>s</td></tr>
@@ -286,6 +287,7 @@ <h3>Run summary</h3>
     "walkers_requested": 100000,
     "walkers_generated": 120000,
     "walkers_remaining": 112724,
+    "seed": 3,
     "duration": 1200,
     "simulated_time": 1199.2085202681737,
     "time_step": 1.0631281208051186,
@@ -320,6 +322,16 @@ <h3>Run summary</h3>
 }
 </pre></div>
 
+<h3>Random numbers and parallel processing</h3>
+
+The walkers are placed and moved using pseudo-random numbers. The seed
+is given by <b>random_seed</b>; without it, a seed is generated and
+recorded in the history of the output maps and in the run summary as
+<code>seed</code>, so that the run can be repeated. With more than one thread, the
+results differ slightly between thread counts and between repeated runs,
+since the order in which the walkers reach a cell depends on the
+threads. Use <b>nprocs=1</b> when results must be reproducible.
+
 <h2>EXAMPLE</h2>
 
 This example uses the
diff --git a/raster/r.sim/r.sim.water/r.sim.water.md b/raster/r.sim/r.sim.water/r.sim.water.md
index ec2fded33af..755767d98be 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.md
+++ b/raster/r.sim/r.sim.water/r.sim.water.md
@@ -196,6 +196,7 @@ stored in the history of the output raster maps under the same keys (see
 | `walkers_requested` | Number of walkers from **nwalkers**, by default twice the number of cells | count |
 | `walkers_generated` | Walkers created, at least one per cell and more where the source rate is higher | count |
 | `walkers_remaining` | Walkers still in the domain at the end of the run | count |
+| `seed` | Seed of the random numbers, given or generated | |
 | `duration` | Requested simulation length (**duration**) | s |
 | `simulated_time` | Simulated time reached at the end of the run | s |
 | `time_step` | Simulated time per iteration | s |
@@ -277,6 +278,7 @@ The printed summary:
     "walkers_requested": 100000,
     "walkers_generated": 120000,
     "walkers_remaining": 112724,
+    "seed": 3,
     "duration": 1200,
     "simulated_time": 1199.2085202681737,
     "time_step": 1.0631281208051186,
@@ -311,6 +313,16 @@ The printed summary:
 }
 ```
 
+### Random numbers and parallel processing
+
+The walkers are placed and moved using pseudo-random numbers. The seed
+is given by **random_seed**; without it, a seed is generated and
+recorded in the history of the output maps and in the run summary as
+`seed`, so that the run can be repeated. With more than one thread, the
+results differ slightly between thread counts and between repeated runs,
+since the order in which the walkers reach a cell depends on the
+threads. Use **nprocs=1** when results must be reproducible.
+
 ## EXAMPLE
 
 This example uses the
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_landforms_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_landforms_test.py
index 57f97cea1ca..52c10ba66e8 100644
--- a/raster/r.sim/r.sim.water/tests/r_sim_water_landforms_test.py
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_landforms_test.py
@@ -173,11 +173,11 @@ def test_pit_collects_water(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [0.50, 0.39, 0.45, 0.64, 0.39],
-            [0.45, 2.90, 14.30, 2.43, 0.80],
-            [0.45, 14.34, 87.17, 14.28, 0.45],
-            [0.50, 2.56, 14.76, 2.85, 0.68],
-            [0.64, 0.50, 0.26, 0.80, 0.50],
+            [0.55, 0.60, 0.45, 0.33, 0.55],
+            [0.45, 2.84, 14.42, 2.67, 0.64],
+            [0.39, 14.55, 87.43, 13.70, 0.45],
+            [0.39, 2.45, 14.30, 2.36, 0.76],
+            [0.60, 0.39, 0.55, 0.60, 0.50],
         ],
     )
 
@@ -225,13 +225,13 @@ def test_peak_sheds_water(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [0.89, 0.68, 0.81, 0.89, 0.82, 0.70, 0.90],
-            [0.69, 0.52, 0.48, 0.57, 0.49, 0.52, 0.74],
-            [0.79, 0.49, 0.38, 0.42, 0.39, 0.52, 0.79],
-            [0.92, 0.58, 0.40, 0.69, 0.42, 0.58, 0.90],
-            [0.81, 0.50, 0.41, 0.42, 0.37, 0.49, 0.80],
-            [0.71, 0.50, 0.48, 0.60, 0.51, 0.51, 0.69],
-            [0.86, 0.67, 0.80, 0.92, 0.80, 0.74, 0.87],
+            [0.86, 0.69, 0.79, 0.92, 0.81, 0.73, 0.86],
+            [0.72, 0.52, 0.50, 0.58, 0.48, 0.52, 0.70],
+            [0.84, 0.50, 0.38, 0.43, 0.40, 0.51, 0.81],
+            [0.89, 0.59, 0.41, 0.73, 0.40, 0.58, 0.90],
+            [0.79, 0.47, 0.36, 0.41, 0.36, 0.49, 0.83],
+            [0.69, 0.51, 0.49, 0.58, 0.49, 0.50, 0.71],
+            [0.87, 0.70, 0.84, 0.89, 0.82, 0.70, 0.86],
         ],
     )
 
@@ -267,11 +267,11 @@ def test_ridge_sheds_water_to_flanks(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [1.00, 0.66, 0.40, 0.57, 0.42, 0.68, 1.02],
-            [1.09, 0.76, 0.58, 1.37, 0.56, 0.77, 1.11],
-            [1.11, 0.73, 0.52, 1.12, 0.50, 0.72, 1.09],
-            [1.13, 0.78, 0.58, 1.45, 0.58, 0.75, 1.11],
-            [1.00, 0.64, 0.40, 0.59, 0.40, 0.66, 1.00],
+            [0.96, 0.66, 0.39, 0.59, 0.40, 0.65, 0.98],
+            [1.15, 0.77, 0.58, 1.37, 0.57, 0.77, 1.14],
+            [1.09, 0.71, 0.50, 1.07, 0.51, 0.72, 1.09],
+            [1.14, 0.77, 0.58, 1.39, 0.58, 0.78, 1.10],
+            [0.98, 0.65, 0.41, 0.58, 0.39, 0.63, 0.98],
         ],
     )
 
@@ -308,16 +308,16 @@ def test_valley_concentrates_flow_downstream(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [0.34, 0.49, 0.77, 2.17, 0.74, 0.47, 0.36],
-            [0.44, 0.60, 0.77, 1.55, 0.80, 0.59, 0.39],
-            [0.43, 0.63, 0.96, 2.42, 0.97, 0.59, 0.43],
-            [0.41, 0.63, 1.03, 3.19, 1.03, 0.64, 0.41],
-            [0.41, 0.65, 1.11, 3.88, 1.08, 0.62, 0.42],
-            [0.43, 0.61, 1.16, 4.46, 1.16, 0.64, 0.41],
-            [0.45, 0.63, 1.21, 5.02, 1.23, 0.64, 0.41],
-            [0.45, 0.65, 1.32, 5.50, 1.27, 0.61, 0.42],
-            [0.44, 0.64, 1.37, 6.04, 1.32, 0.63, 0.47],
-            [0.44, 0.73, 1.54, 6.61, 1.55, 0.77, 0.50],
+            [0.37, 0.48, 0.73, 2.11, 0.69, 0.47, 0.35],
+            [0.43, 0.60, 0.80, 1.54, 0.79, 0.56, 0.40],
+            [0.43, 0.64, 0.97, 2.38, 0.95, 0.63, 0.43],
+            [0.41, 0.62, 1.02, 3.17, 1.04, 0.64, 0.44],
+            [0.41, 0.60, 1.11, 3.84, 1.11, 0.60, 0.41],
+            [0.42, 0.64, 1.16, 4.44, 1.18, 0.64, 0.42],
+            [0.43, 0.64, 1.21, 5.01, 1.23, 0.61, 0.41],
+            [0.43, 0.63, 1.27, 5.53, 1.25, 0.60, 0.44],
+            [0.44, 0.63, 1.33, 6.03, 1.30, 0.68, 0.42],
+            [0.50, 0.75, 1.51, 6.58, 1.51, 0.74, 0.48],
         ],
     )
 
@@ -360,13 +360,13 @@ def test_saddle_bifurcates_flow(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [0.26, 0.28, 0.15, 0.20, 0.24, 0.15, 0.26],
-            [0.41, 0.30, 0.28, 0.39, 0.22, 0.51, 0.49],
-            [0.72, 0.51, 0.46, 0.81, 0.33, 0.56, 0.85],
-            [1.30, 0.83, 0.66, 3.08, 1.07, 0.95, 1.39],
-            [0.78, 0.54, 0.30, 0.78, 0.47, 0.61, 0.84],
-            [0.43, 0.35, 0.30, 0.41, 0.28, 0.43, 0.45],
-            [0.30, 0.20, 0.26, 0.20, 0.24, 0.28, 0.24],
+            [0.38, 0.18, 0.35, 0.15, 0.24, 0.18, 0.33],
+            [0.45, 0.43, 0.42, 0.24, 0.38, 0.30, 0.46],
+            [0.85, 0.64, 0.38, 0.64, 0.51, 0.59, 0.68],
+            [1.41, 0.83, 0.77, 3.21, 0.77, 0.90, 1.42],
+            [0.69, 0.70, 0.41, 0.91, 0.45, 0.66, 0.74],
+            [0.43, 0.33, 0.41, 0.33, 0.42, 0.20, 0.46],
+            [0.36, 0.30, 0.22, 0.20, 0.26, 0.18, 0.45],
         ],
     )
 
@@ -414,32 +414,32 @@ def test_hollow_concentrates_more_than_spur(tmp_path):
     assert_depth_pinned(
         depth_hollow,
         [
-            [0.34, 0.54, 0.79, 0.67, 0.67, 0.52, 0.36],
-            [0.41, 0.57, 0.96, 1.22, 0.90, 0.64, 0.46],
-            [0.47, 0.71, 1.25, 2.08, 1.20, 0.73, 0.37],
-            [0.41, 0.60, 1.46, 2.84, 1.33, 0.57, 0.42],
-            [0.48, 0.65, 1.46, 3.42, 1.43, 0.71, 0.46],
-            [0.38, 0.73, 1.56, 4.01, 1.56, 0.76, 0.42],
-            [0.43, 0.68, 1.77, 4.53, 1.72, 0.74, 0.39],
-            [0.46, 0.65, 1.79, 4.99, 1.68, 0.65, 0.43],
-            [0.37, 0.71, 1.87, 5.54, 1.87, 0.80, 0.42],
-            [0.43, 0.79, 2.18, 6.95, 2.21, 0.77, 0.43],
+            [0.32, 0.44, 0.66, 0.94, 0.77, 0.47, 0.40],
+            [0.42, 0.63, 0.89, 1.43, 0.97, 0.64, 0.44],
+            [0.50, 0.69, 1.24, 2.02, 1.23, 0.66, 0.38],
+            [0.47, 0.66, 1.34, 2.83, 1.41, 0.60, 0.43],
+            [0.43, 0.85, 1.57, 3.41, 1.37, 0.65, 0.44],
+            [0.42, 0.59, 1.54, 4.10, 1.55, 0.81, 0.43],
+            [0.45, 0.70, 1.57, 4.55, 1.78, 0.66, 0.47],
+            [0.47, 0.67, 1.65, 5.12, 1.75, 0.70, 0.37],
+            [0.44, 0.77, 1.79, 5.59, 1.77, 0.63, 0.48],
+            [0.52, 0.77, 2.24, 6.88, 2.12, 0.84, 0.49],
         ],
     )
 
     assert_depth_pinned(
         depth_spur,
         [
-            [0.73, 0.46, 0.29, 0.45, 0.36, 0.51, 0.55],
-            [0.74, 0.63, 0.56, 0.64, 0.45, 0.56, 0.80],
-            [0.97, 0.77, 0.64, 0.78, 0.56, 0.72, 1.09],
-            [1.08, 0.77, 0.61, 0.89, 0.54, 0.68, 1.00],
-            [1.13, 0.80, 0.64, 0.92, 0.71, 0.66, 0.98],
-            [1.08, 0.83, 0.69, 0.96, 0.67, 0.80, 1.11],
-            [1.19, 0.75, 0.77, 1.07, 0.77, 0.76, 1.16],
-            [1.14, 0.85, 0.73, 1.01, 0.81, 0.75, 1.09],
-            [1.20, 0.72, 0.81, 1.05, 0.79, 0.91, 1.26],
-            [1.40, 1.09, 0.93, 1.24, 0.81, 0.88, 1.29],
+            [0.68, 0.41, 0.48, 0.51, 0.43, 0.47, 0.66],
+            [0.89, 0.67, 0.49, 0.59, 0.54, 0.70, 0.81],
+            [0.87, 0.59, 0.71, 0.78, 0.48, 0.71, 1.02],
+            [1.01, 0.84, 0.74, 0.91, 0.60, 0.68, 1.00],
+            [1.23, 0.83, 0.74, 0.91, 0.62, 0.73, 0.94],
+            [1.17, 0.79, 0.74, 0.94, 0.73, 0.77, 1.20],
+            [1.07, 0.77, 0.78, 0.91, 0.75, 0.80, 1.07],
+            [1.11, 0.91, 0.61, 0.95, 0.74, 0.84, 1.22],
+            [1.27, 0.76, 0.71, 0.94, 0.67, 0.81, 1.15],
+            [1.38, 1.04, 0.87, 1.20, 0.78, 0.94, 1.40],
         ],
     )
 
@@ -482,16 +482,16 @@ def test_terrace_accumulates_on_steps(tmp_path):
     assert_depth_pinned(
         depth,
         [
-            [0.83, 0.92, 2.27, 0.31, 0.48, 1.10, 4.06, 0.78, 1.01, 1.93, 8.25, 4.35],
-            [1.78, 2.47, 2.73, 0.97, 1.05, 1.96, 6.18, 1.57, 1.33, 3.01, 12.26, 7.19],
-            [1.98, 2.97, 3.50, 0.78, 0.67, 1.43, 4.73, 1.43, 1.33, 3.37, 13.64, 8.68],
-            [3.56, 5.51, 3.31, 0.83, 0.88, 1.78, 6.33, 1.36, 1.43, 3.96, 15.10, 9.38],
-            [3.82, 5.70, 3.64, 0.97, 0.97, 2.45, 7.33, 1.93, 1.43, 3.96, 17.13, 10.67],
-            [4.35, 6.22, 3.58, 0.72, 0.83, 2.17, 5.73, 1.22, 1.40, 4.03, 18.60, 11.42],
-            [3.76, 5.41, 3.58, 0.83, 0.67, 1.72, 5.00, 1.33, 1.22, 3.50, 16.33, 10.49],
-            [3.52, 3.64, 3.17, 0.67, 0.72, 1.78, 5.99, 1.10, 1.05, 2.97, 12.78, 9.85],
-            [3.35, 3.37, 2.77, 0.55, 0.67, 1.98, 5.48, 1.47, 1.50, 2.81, 10.44, 8.40],
-            [2.50, 1.50, 2.27, 0.55, 0.48, 1.22, 3.54, 1.18, 1.14, 2.17, 7.91, 4.13],
+            [1.72, 3.47, 3.01, 0.72, 0.72, 1.60, 2.73, 0.88, 0.48, 1.69, 6.92, 3.87],
+            [1.96, 4.40, 2.52, 0.55, 0.88, 1.66, 3.96, 1.22, 1.40, 3.31, 13.03, 8.10],
+            [2.59, 4.41, 3.31, 0.92, 0.88, 2.07, 7.22, 1.43, 1.43, 3.31, 14.28, 9.66],
+            [2.17, 3.50, 2.68, 0.67, 0.83, 2.01, 6.20, 1.50, 1.53, 3.78, 16.37, 10.49],
+            [3.03, 4.53, 2.45, 0.72, 0.78, 2.04, 7.47, 1.43, 1.29, 4.15, 17.83, 11.21],
+            [3.47, 4.72, 2.75, 0.92, 0.97, 2.04, 6.24, 1.60, 1.47, 4.10, 17.20, 11.12],
+            [3.56, 4.91, 3.58, 0.61, 0.78, 1.96, 5.18, 1.14, 0.92, 3.64, 16.76, 10.80],
+            [1.69, 3.92, 1.66, 0.48, 0.72, 2.01, 6.33, 1.47, 1.36, 3.47, 12.64, 8.50],
+            [3.50, 2.47, 2.68, 0.83, 0.92, 2.12, 5.88, 1.63, 1.36, 2.77, 9.52, 6.49],
+            [2.30, 2.15, 1.87, 0.83, 0.48, 1.40, 4.43, 1.14, 1.14, 2.01, 5.54, 3.64],
         ],
     )
 
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_random_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_random_test.py
new file mode 100644
index 00000000000..ee6064d1eec
--- /dev/null
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_random_test.py
@@ -0,0 +1,171 @@
+# SPDX-FileCopyrightText: 2026 GRASS Development Team
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+"""Tests of the random numbers of r.sim.water.
+
+Each walker draws from a random number state of its own, set by the seed and
+the walker's number, so its random numbers do not depend on the number of
+threads. The water depth still does at nprocs > 1, because walkers in the same
+cell add to the depth grid without synchronization, so the tests which compare
+numbers of threads compare the walker positions instead. With a large hmax and
+no infiltration, a walker's path depends only on the slope and its random
+numbers, not on the depth. The elevation is a bowl, so that the walkers stay
+in the region until the end of the simulation.
+"""
+
+import os
+import re
+from io import StringIO
+
+import numpy as np
+import pytest
+
+import grass.script as gs
+import grass.script.array as garray
+from grass.exceptions import CalledModuleError
+from grass.tools import Tools
+
+# Larger than any depth the simulation reaches, so the diffusion does not
+# depend on the depth.
+NO_DEPTH_FEEDBACK = {"hmax": 1e6}
+
+
+@pytest.fixture(scope="module")
+def session(tmp_path_factory):
+    """Session in an XY project with a bowl as elevation."""
+    project = tmp_path_factory.mktemp("r_sim_water_random") / "xy_test"
+    gs.create_project(project)
+    with (
+        gs.setup.init(project, env=os.environ.copy()) as session,
+        Tools(session=session) as tools,
+    ):
+        tools.g_region(s=0, n=40, w=0, e=40, res=1)
+        tools.r_mapcalc(
+            expression="elevation = 0.01 * ((col() - 20.5)^2 + (row() - 20.5)^2)"
+        )
+        yield session
+
+
+def simulate_depth(session, **kwargs):
+    """Run the simulation and return the water depth as an array."""
+    tools = Tools(session=session)
+    return tools.r_sim_water(
+        elevation="elevation", depth=np.array, duration=5, **kwargs
+    )
+
+
+def simulate_walkers(session, name, **kwargs):
+    """Run the simulation and return the final walker coordinates as an array."""
+    tools = Tools(session=session)
+    tools.r_sim_water(elevation="elevation", walkers_output=name, duration=5, **kwargs)
+    text = tools.v_out_ascii(input=name, format="point", precision=15).text
+    return np.loadtxt(StringIO(text), delimiter="|")
+
+
+def test_same_seed_gives_same_depth(session):
+    """Two runs with the same seed give the same water depth."""
+    first = simulate_depth(session, random_seed=1)
+    second = simulate_depth(session, random_seed=1)
+    assert np.array_equal(first, second)
+
+
+def test_different_seeds_give_different_depth(session):
+    """Runs with different seeds give different water depths."""
+    first = simulate_depth(session, random_seed=1)
+    second = simulate_depth(session, random_seed=2)
+    assert not np.array_equal(first, second)
+
+
+def test_generated_seed_is_the_seed_option(session):
+    """Without a seed, GRASS_RANDOM_SEED gives what the seed option gives."""
+    env = session.env.copy()
+    env["GRASS_RANDOM_SEED"] = "3"
+    generated = Tools(env=env).r_sim_water(
+        elevation="elevation", depth=np.array, duration=5
+    )
+    given = simulate_depth(session, random_seed=3)
+    assert np.array_equal(generated, given)
+
+
+def test_generated_seed_is_recorded(session):
+    """A run without a seed records the seed it used, and the seed repeats it."""
+    tools = Tools(session=session)
+    tools.r_sim_water(elevation="elevation", depth="generated", duration=5)
+    history = tools.r_info(map="generated", flags="h").text
+    seed = int(re.search(r"\bseed=(-?\d+)", history).group(1))
+    repeated = simulate_depth(session, random_seed=seed)
+    generated = garray.array("generated", env=session.env)
+    assert np.array_equal(repeated, generated)
+
+
+def test_deprecated_flag_generates_a_seed(session):
+    """The deprecated -s flag still gives what GRASS_RANDOM_SEED says."""
+    env = session.env.copy()
+    env["GRASS_RANDOM_SEED"] = "3"
+    generated = Tools(env=env).r_sim_water(
+        elevation="elevation", depth=np.array, duration=5, flags="s"
+    )
+    given = simulate_depth(session, random_seed=3)
+    assert np.array_equal(generated, given)
+
+
+def test_seed_and_flag_are_exclusive(session):
+    """The seed option and the flag to generate a seed cannot be combined."""
+    with pytest.raises(CalledModuleError, match="mutually exclusive"):
+        simulate_depth(session, random_seed=1, flags="s")
+
+
+@pytest.mark.parametrize("nprocs", [2, 4])
+def test_depth_does_not_depend_on_nprocs_below_hmax(session, nprocs):
+    """The depth is identical with any number of threads while it stays below hmax.
+
+    Above hmax a walker reacts to the depth the walkers before it in the
+    same step left, and that order depends on the threads.
+    """
+    serial = simulate_depth(session, random_seed=5, **NO_DEPTH_FEEDBACK)
+    parallel = simulate_depth(
+        session, random_seed=5, nprocs=nprocs, **NO_DEPTH_FEEDBACK
+    )
+    assert np.array_equal(parallel, serial)
+
+
+@pytest.mark.parametrize("nprocs", [2, 4])
+def test_walkers_do_not_depend_on_nprocs(session, nprocs):
+    """Walkers end at the same positions with any number of threads."""
+    serial = simulate_walkers(
+        session, f"walkers_{nprocs}_serial", random_seed=5, **NO_DEPTH_FEEDBACK
+    )
+    parallel = simulate_walkers(
+        session,
+        f"walkers_{nprocs}_parallel",
+        random_seed=5,
+        nprocs=nprocs,
+        **NO_DEPTH_FEEDBACK,
+    )
+    assert serial.size
+    assert np.array_equal(parallel, serial)
+
+
+def test_walkers_depend_on_seed(session):
+    """The walker comparison above can detect a difference."""
+    first = simulate_walkers(
+        session, "walkers_seed_5", random_seed=5, **NO_DEPTH_FEEDBACK
+    )
+    second = simulate_walkers(
+        session, "walkers_seed_6", random_seed=6, **NO_DEPTH_FEEDBACK
+    )
+    assert not np.array_equal(first, second)
+
+
+@pytest.mark.parametrize("seed", [-(2**31) - 1, 2**32])
+def test_seed_outside_range_is_an_error(session, seed):
+    """A seed the generator cannot use is refused, not silently wrapped."""
+    with pytest.raises(CalledModuleError, match="outside the range"):
+        simulate_depth(session, random_seed=seed)
+
+
+@pytest.mark.parametrize("seed", ["12abc", "1.5", "-", "99999999999999999999"])
+def test_seed_which_is_not_an_integer_is_an_error(session, seed):
+    """A seed which the parser lets through but which is not an integer is refused."""
+    with pytest.raises(CalledModuleError, match="Invalid random seed"):
+        simulate_depth(session, random_seed=seed)
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
index d62e3b5abc6..17cfb806d86 100644
--- a/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_summary_test.py
@@ -24,6 +24,7 @@
     "walkers_requested",
     "walkers_generated",
     "walkers_remaining",
+    "seed",
     "duration",
     "simulated_time",
     "time_step",
@@ -91,6 +92,7 @@ def test_json_summary(session_tools):
     assert summary["mean_source_rate"] > 0
     assert summary["mean_infiltration"] == 0
     assert summary["threads"] == 1
+    assert summary["seed"] == 1
 
     assert len(summary["outputs"]) == 1
     output = summary["outputs"][0]
@@ -195,6 +197,7 @@ def test_history(session_tools):
             "walkers_requested",
             "walkers_generated",
             "walkers_remaining",
+            "seed",
             "duration",
         ):
             assert history_value(history, key) == summary[key]
diff --git a/raster/r.sim/r.sim.water/tests/r_sim_water_test.py b/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
index 0e6f67a61f3..e0b471ae981 100644
--- a/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
+++ b/raster/r.sim/r.sim.water/tests/r_sim_water_test.py
@@ -750,14 +750,18 @@ def test_flow_control_increases_depth(east_slope_session):
     Higher trapping probability means more accumulation.
     """
     tools = Tools(session=east_slope_session)
-    sum_no_control = float(np.sum(run_sim(east_slope_session)))
+    # The ten default walkers leave the total with trapping at 0.3 within
+    # noise of the total without trapping for about a quarter of the seeds;
+    # a thousand separate no trapping, 0.3 and 0.8 for every seed tried.
+    walkers = {"nwalkers": 1000}
+    sum_no_control = float(np.sum(run_sim(east_slope_session, **walkers)))
     tools.r_mapcalc(expression="flow_ctrl_low = 0.3")
     tools.r_mapcalc(expression="flow_ctrl_high = 0.8")
     sum_low_trap = float(
-        np.sum(run_sim(east_slope_session, flow_control="flow_ctrl_low"))
+        np.sum(run_sim(east_slope_session, flow_control="flow_ctrl_low", **walkers))
     )
     sum_high_trap = float(
-        np.sum(run_sim(east_slope_session, flow_control="flow_ctrl_high"))
+        np.sum(run_sim(east_slope_session, flow_control="flow_ctrl_high", **walkers))
     )
     assert sum_low_trap > sum_no_control, (
         f"Trapping should increase depth: "
@@ -907,8 +911,10 @@ def test_walkers_output(long_slope_session):
 def test_nprocs_gives_result_within_noise(east_slope_session):
     """Multiple threads produce a result close to a single thread.
 
-    Threads share the random number generator state, so multi-threaded
-    results vary between runs even with a fixed seed. Total depth should
+    Each walker draws the same random numbers with any number of threads,
+    but threads add to the water depth without synchronization, so
+    multi-threaded results vary between runs even with a fixed seed. Total
+    depth should
     still agree with the single-threaded result within Monte Carlo noise,
     using the same tolerance as test_results_consistent_across_seeds.
     """
diff --git a/raster/r.sim/r.sim.water/testsuite/data/depth_complex.pack b/raster/r.sim/r.sim.water/testsuite/data/depth_complex.pack
index 30c7b86cead..09d95e3884e 100644
Binary files a/raster/r.sim/r.sim.water/testsuite/data/depth_complex.pack and b/raster/r.sim/r.sim.water/testsuite/data/depth_complex.pack differ
diff --git a/raster/r.sim/r.sim.water/testsuite/data/depth_default.pack b/raster/r.sim/r.sim.water/testsuite/data/depth_default.pack
index e8831bd0a23..c3f629b55dd 100644
Binary files a/raster/r.sim/r.sim.water/testsuite/data/depth_default.pack and b/raster/r.sim/r.sim.water/testsuite/data/depth_default.pack differ
diff --git a/raster/r.sim/r.sim.water/testsuite/data/discharge_complex.pack b/raster/r.sim/r.sim.water/testsuite/data/discharge_complex.pack
index 0025abf887d..eef42412f17 100644
Binary files a/raster/r.sim/r.sim.water/testsuite/data/discharge_complex.pack and b/raster/r.sim/r.sim.water/testsuite/data/discharge_complex.pack differ
diff --git a/raster/r.sim/r.sim.water/testsuite/data/discharge_default.pack b/raster/r.sim/r.sim.water/testsuite/data/discharge_default.pack
index f2c2861a141..4a713cbd478 100644
Binary files a/raster/r.sim/r.sim.water/testsuite/data/discharge_default.pack and b/raster/r.sim/r.sim.water/testsuite/data/discharge_default.pack differ
diff --git a/raster/r.sim/r.sim.water/testsuite/test_r_sim_water.py b/raster/r.sim/r.sim.water/testsuite/test_r_sim_water.py
index b765f5b6d2a..5d9f0dafca6 100644
--- a/raster/r.sim/r.sim.water/testsuite/test_r_sim_water.py
+++ b/raster/r.sim/r.sim.water/testsuite/test_r_sim_water.py
@@ -252,7 +252,7 @@ def test_default(self):
             random_seed=1,
         )
         self.assertRasterFitsUnivar(
-            self.depth, reference="sum=30423.190201", precision=1e-6
+            self.depth, reference="sum=30424.022289", precision=1e-6
         )
 
 
diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c
index 457872cf865..e40479ee07a 100644
--- a/raster/r.sim/simlib/hydro.c
+++ b/raster/r.sim/simlib/hydro.c
@@ -27,9 +27,6 @@
  *
  */
 
-struct point2D;
-struct point3D;
-
 /* **************************************************** */
 /*       create walker representation of si */
 /* ******************************************************** */
@@ -71,6 +68,31 @@ void main_loop(const Setup *setup, const Geometry *geometry,
     G_debug(2, " maxwa, nblock %d %d", sim->maxwa, nblock);
     G_debug(2, "rwalk, sisum: %f %f", sim->rwalk, setup->sisum);
 
+    // Each walker draws from a random number state of its own, so its
+    // values depend on the seed and the walker's number, not on the
+    // thread which moves it. The bound is what a walker draws: two values
+    // when placed and, in each time step, a normal pair by the polar
+    // method (8 / pi values on average, two per attempt at a pair) and one
+    // more in a trap. Eight per step allow three attempts at the pair,
+    // which the sum over many steps rarely exceeds, and the 48 more allow
+    // a run of 24 rejected attempts within one step: the probability per
+    // walker is below 1e-16, so a run of 2^31 walkers (the most a 32-bit
+    // int count holds) exceeds the bound with a probability below 1e-6,
+    // one in a million.
+    struct G_random_layout layout;
+
+    G_random_init_layout_bounded(&layout, settings->seed, sim->max_walkers,
+                                 2 + 8LL * setup->miter + 48);
+    if (G_random_layout_batches(&layout) < 1)
+        G_warning(_("%d walkers over %d time steps may draw more random "
+                    "numbers than one seed provides; the random numbers "
+                    "of some walkers then repeat those of others shifted "
+                    "by a constant"),
+                  sim->max_walkers, setup->miter);
+#pragma omp parallel for
+    for (int lw = 0; lw < sim->max_walkers; lw++)
+        G_random_state_for_unit(&sim->w[lw].state, &layout, lw);
+
     for (iblock = 1; iblock <= nblock; iblock++) {
         int lw = 0;
         double walkwe = 0.;
@@ -92,10 +114,12 @@ void main_loop(const Setup *setup, const Geometry *geometry,
 
                     for (int iw = 1; iw <= mgen + 1;
                          iw++) { /* assign walkers */
-                        sim->w[lw].x =
-                            x + geometry->stepx * (simwe_rand() - 0.5);
-                        sim->w[lw].y =
-                            y + geometry->stepy * (simwe_rand() - 0.5);
+                        struct G_random_state *state = &sim->w[lw].state;
+
+                        sim->w[lw].x = x + geometry->stepx *
+                                               (G_random_double(state) - 0.5);
+                        sim->w[lw].y = y + geometry->stepy *
+                                               (G_random_double(state) - 0.5);
                         sim->w[lw].m = wei;
 
                         walkwe += sim->w[lw].m;
@@ -235,12 +259,7 @@ void main_loop(const Setup *setup, const Geometry *geometry,
 
                             double d1 = gama * conn;
                             double gaux, gauy;
-#if defined(_OPENMP)
-                            gasdev_for_paralel(&gaux, &gauy);
-#else
-                            gaux = gasdev();
-                            gauy = gasdev();
-#endif
+                            gasdev(&sim->w[lw].state, &gaux, &gauy);
                             double hhc = pow(d1, 3. / 5.);
                             double velx, vely;
                             /* Diffusion coefficient of this walker's move */
@@ -261,7 +280,7 @@ void main_loop(const Setup *setup, const Geometry *geometry,
                             if (inputs->traps != NULL &&
                                 grids->trap[k][l] != 0.) { /* traps */
 
-                                float eff = simwe_rand(); /* random generator */
+                                float eff = G_random_double(&sim->w[lw].state);
 
                                 if (eff <= grids->trap[k][l]) {
                                     velx = -0.1 *
diff --git a/raster/r.sim/simlib/input.c b/raster/r.sim/simlib/input.c
index ddb6298fbe8..f41fe03305b 100644
--- a/raster/r.sim/simlib/input.c
+++ b/raster/r.sim/simlib/input.c
@@ -62,8 +62,9 @@ void alloc_walkers(int max_walkers, Simulation *sim, const Outputs *outputs)
 {
     G_debug(1, "beginning memory allocation for walkers");
 
-    sim->w = (struct point3D *)G_calloc(max_walkers, sizeof(struct point3D));
+    sim->w = (struct walker *)G_calloc(max_walkers, sizeof(struct walker));
     sim->vavg = (struct point2D *)G_calloc(max_walkers, sizeof(struct point2D));
+    sim->max_walkers = max_walkers;
     if (outputs->outwalk != NULL)
         sim->stack =
             (struct point3D *)G_calloc(max_walkers, sizeof(struct point3D));
diff --git a/raster/r.sim/simlib/output.c b/raster/r.sim/simlib/output.c
index 393e2b12eb7..d7f5986ef14 100644
--- a/raster/r.sim/simlib/output.c
+++ b/raster/r.sim/simlib/output.c
@@ -1,5 +1,6 @@
 /* output.c (simlib), 20.nov.2002, JH */
 
+#include <inttypes.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -107,6 +108,7 @@ static void write_history(const char *name, double tt, const Setup *setup,
         &hist,
         "walkers_generated=%d, walkers_requested=%d, walkers_remaining=%d",
         sim->nwalk, sim->maxwa, sim->nwalka);
+    Rast_append_format_history(&hist, "seed=%" PRId64, settings->seed);
     Rast_append_format_history(&hist, "duration=%d, simulated_time=%f",
                                settings->timesec, tt);
     Rast_append_format_history(&hist, "time_step=%f, mean_velocity=%f",
diff --git a/raster/r.sim/simlib/random.c b/raster/r.sim/simlib/random.c
index 8a0774695fc..53db97a0917 100644
--- a/raster/r.sim/simlib/random.c
+++ b/raster/r.sim/simlib/random.c
@@ -1,55 +1,81 @@
 /* random.c (simlib), 20.nov.2002, JH */
 
-#include <stdio.h>
-#include <stdlib.h>
+#include <errno.h>
+#include <inttypes.h>
 #include <math.h>
+#include <stdlib.h>
+
 #include <grass/gis.h>
-#include <grass/bitmap.h>
-#include <grass/linkm.h>
+#include <grass/glocale.h>
 
-double simwe_rand(void)
-{
-    return G_drand48();
-} /* ulec */
+#include <grass/simlib.h>
 
-double gasdev(void)
+/*!
+ * \brief Return the seed of the walkers' random numbers
+ *
+ * The seed is read from the option, or generated when the option is not
+ * given, as in r.mapcalc. A seed which is not an integer or is outside the
+ * range the generator accepts is a fatal error, here, before any input is
+ * read. The flag to generate a seed is deprecated, since generating is the
+ * default; the flag together with the option is a fatal error.
+ *
+ * \param seed the seed option
+ * \param generate the deprecated flag to generate a seed
+ *
+ * \return the seed
+ */
+int64_t simwe_seed(const struct Option *seed, const struct Flag *generate)
 {
-    /* Initialized data */
+    int64_t value;
+    struct G_random_state check;
 
-    static int iset = 0;
-    static double gset = .1;
+    if (seed->answer && generate->answer)
+        G_fatal_error(_("%s= and -%c are mutually exclusive"), seed->key,
+                      generate->key);
+    if (generate->answer)
+        G_verbose_message(_("Flag '%c' is deprecated and will be removed in a "
+                            "future release. Seeding is automatic or use "
+                            "parameter %s."),
+                          generate->key, seed->key);
+    if (seed->answer) {
+        long long parsed;
+        char *end;
 
-    /* System generated locals */
-    double ret_val;
-
-    /* Local variables */
-    double r = 0.0, vv1 = 0.0, vv2 = 0.0, fac = 0.0;
-
-    if (iset == 0) {
-        while (r >= 1. || r == 0.) {
-            vv1 = simwe_rand() * 2. - 1.;
-            vv2 = simwe_rand() * 2. - 1.;
-            r = vv1 * vv1 + vv2 * vv2;
-        }
-        fac = sqrt(log(r) * -2. / r);
-        gset = vv1 * fac;
-        ret_val = vv2 * fac;
-        iset = 1;
+        errno = 0;
+        parsed = strtoll(seed->answer, &end, 10);
+        if (end == seed->answer || *end != '\0' || errno == ERANGE)
+            G_fatal_error(_("Invalid random seed <%s>"), seed->answer);
+        value = parsed;
+        G_verbose_message(_("Read random seed from %s option: %" PRId64),
+                          seed->key, value);
     }
     else {
-        ret_val = gset;
-        iset = 0;
+        value = G_random_generate_seed();
+        G_verbose_message(_("Generated random seed: %" PRId64), value);
     }
-    return ret_val;
-} /* gasdev */
+    /* The layout, which refuses a seed out of range, is built only when the
+     * number of time steps is known. */
+    G_random_state_from_seed(&check, value);
+    return value;
+}
 
-void gasdev_for_paralel(double *x, double *y)
+/*!
+ * \brief Draw a pair of independent standard normal values
+ *
+ * Uses the polar method, which rejects points outside the unit circle, so
+ * the number of values drawn varies, 8 / pi on average.
+ *
+ * \param state the random number state of the walker
+ * \param[out] x the first value
+ * \param[out] y the second value
+ */
+void gasdev(struct G_random_state *state, double *x, double *y)
 {
     double r = 0.0, vv1 = 0.0, vv2 = 0.0, fac = 0.0;
 
     while (r >= 1. || r == 0.) {
-        vv1 = simwe_rand() * 2. - 1.;
-        vv2 = simwe_rand() * 2. - 1.;
+        vv1 = G_random_double(state) * 2. - 1.;
+        vv2 = G_random_double(state) * 2. - 1.;
         r = vv1 * vv1 + vv2 * vv2;
     }
     fac = sqrt(log(r) * -2. / r);
diff --git a/raster/r.sim/simlib/simlib.h b/raster/r.sim/simlib/simlib.h
index fe71bea6012..884ef17e0b6 100644
--- a/raster/r.sim/simlib/simlib.h
+++ b/raster/r.sim/simlib/simlib.h
@@ -5,8 +5,11 @@
  * \brief This is the interface for the simlib (SIMWE) library.
  */
 #include <stdbool.h>
+#include <stdint.h>
 #include <stdio.h>
 
+#include <grass/gis.h>
+
 #define EPS         1.e-7
 #define UNDEF       -9999
 
@@ -37,6 +40,7 @@ typedef struct {
     int timesec;        // Time how long the simulation runs [seconds]
     bool ts;            // Time series output
     double mintimestep; // Minimum time step for the simulation [seconds]
+    int64_t seed;       // Seed of the walkers' random numbers
 } Settings;
 
 typedef struct {
@@ -59,10 +63,10 @@ typedef struct {
     int nstack;            // Number of output walkers
     struct point3D *stack; // Output 3D walkers
     int maxwa;             // Number of total walkers
-    double rwalk;      // Number of input walkers per block as double precision
-    struct point3D *w; // Weight of walkers
+    double rwalk; // Number of input walkers per block as double precision
+    struct walker *w;
     struct point2D *vavg; // Average velocity of walkers
-
+    int max_walkers; // Number of allocated walkers, maxwa plus one per cell
 } Simulation;
 
 typedef struct {
@@ -161,6 +165,14 @@ struct point3D {
     double m;
 };
 
+// A walker with its weight m and its own random number state
+struct walker {
+    double x;
+    double y;
+    double m;
+    struct G_random_state state;
+};
+
 void alloc_grids_water(const Geometry *geometry, const Outputs *outputs,
                        Grids *grids);
 void alloc_grids_sediment(const Geometry *geometry, const Outputs *outputs,
@@ -197,9 +209,8 @@ void create_observation_points(ObservationPoints *points);
 void derivatives(const Geometry *geometry, float **elevation, double **dx,
                  double **dy);
 
-double simwe_rand(void);
-double gasdev(void);
-void gasdev_for_paralel(double *, double *);
+int64_t simwe_seed(const struct Option *seed, const struct Flag *generate);
+void gasdev(struct G_random_state *state, double *x, double *y);
 double amax1(double, double);
 double amin1(double, double);
 int min(int, int);
diff --git a/raster/r.sim/simlib/summary.c b/raster/r.sim/simlib/summary.c
index 00ab21046b3..866c4282f7f 100644
--- a/raster/r.sim/simlib/summary.c
+++ b/raster/r.sim/simlib/summary.c
@@ -8,6 +8,7 @@
  *
  *****************************************************************************/
 
+#include <inttypes.h>
 #include <math.h>
 #include <stdio.h>
 
@@ -75,6 +76,7 @@ static void print_json(const Setup *setup, const Settings *settings,
     G_json_object_set_number(root, "walkers_requested", sim->maxwa);
     G_json_object_set_number(root, "walkers_generated", sim->nwalk);
     G_json_object_set_number(root, "walkers_remaining", sim->nwalka);
+    G_json_object_set_number(root, "seed", (double)settings->seed);
     G_json_object_set_number(root, "duration", settings->timesec);
     G_json_object_set_number(
         root, "simulated_time",
@@ -164,6 +166,7 @@ static void print_plain(const Setup *setup, const Settings *settings,
     printf("walkers_requested: %d\n", sim->maxwa);
     printf("walkers_generated: %d\n", sim->nwalk);
     printf("walkers_remaining: %d\n", sim->nwalka);
+    printf("seed: %" PRId64 "\n", settings->seed);
     printf("duration: %d\n", settings->timesec);
     print_number("simulated_time",
                  simulated_seconds(setup, summary->iterations_completed));