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 1234fa1f092..082cf8484fd 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.html
+++ b/raster/r.sim/r.sim.water/r.sim.water.html
@@ -73,8 +73,10 @@ <h2>DESCRIPTION</h2>
 
 <p>
 Output includes a water depth raster map <b>depth</b> in [m], and a water discharge
-raster map <b>discharge</b> in [m3/s]. Error of the numerical solution can be analyzed using
-the <b>error</b> raster map (the resulting water depth is an average, and err is its RMSE).
+raster map <b>discharge</b> in [m3/s]. The <b>error</b> raster map is a Monte Carlo
+standard-deviation estimator across replicas of the particle simulation; the simulation
+currently runs a single replica, so this map is zero everywhere and is provided for
+forward compatibility with planned multiple-replica execution.
 The output vector points map <b>output_walkers</b> can be used to analyze and visualize
 spatial distribution of walkers at different simulation times (note that
 the resulting water depth is based on the density of these walkers).
@@ -83,11 +85,6 @@ <h2>DESCRIPTION</h2>
 <!-- from
 http://www.ing.unitn.it/~grass/conferences/GRASS2002/proceedings/proceedings/pdfs/Mitasova_Helena_2.pdf
 -->
-The spatial distribution of numerical error associated with path sampling solution can be
-analysed using the output error raster file [m]. This error is a function of the number
-of particles used in the simulation and can be reduced by increasing the number of walkers
-given by parameter <b>nwalkers</b>.
-<!--(<font color="#ff0000"> toto treba upresnit/zmenit, lebo nwalk ide prec</font>). -->
 Duration of simulation is controlled by the <b>duration</b> parameter. The default value
 is 10 minutes, reaching the steady-state may require much longer time,
 depending on the time step, complexity of terrain, land cover and size of the area.
@@ -240,16 +237,6 @@ <h2>EXAMPLE</h2>
 </i>
 </div>
 
-<h2>ERROR MESSAGES</h2>
-
-If the module fails with
-
-<div class="code"><pre>
-ERROR: nwalk (7000001) &gt; maxw (7000000)!
-</pre></div>
-
-then a lower <b>nwalkers</b> parameter value has to be selected.
-
 <h2>REFERENCES</h2>
 
 <ul>
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 9316a56bd77..527252899ef 100644
--- a/raster/r.sim/r.sim.water/r.sim.water.md
+++ b/raster/r.sim/r.sim.water/r.sim.water.md
@@ -64,17 +64,15 @@ defines the probability of particles to pass through the structure (the
 values will be 0-1).
 
 Output includes a water depth raster map **depth** in \[m\], and a water
-discharge raster map **discharge** in \[m3/s\]. Error of the numerical
-solution can be analyzed using the **error** raster map (the resulting
-water depth is an average, and err is its RMSE). The output vector
-points map **output_walkers** can be used to analyze and visualize
-spatial distribution of walkers at different simulation times (note that
-the resulting water depth is based on the density of these walkers). The
-spatial distribution of numerical error associated with path sampling
-solution can be analysed using the output error raster file \[m\]. This
-error is a function of the number of particles used in the simulation
-and can be reduced by increasing the number of walkers given by
-parameter **nwalkers**. Duration of simulation is controlled by the
+discharge raster map **discharge** in \[m3/s\]. The **error** raster map
+is a Monte Carlo standard-deviation estimator across replicas of the
+particle simulation; the simulation currently runs a single replica, so
+this map is zero everywhere and is provided for forward compatibility
+with planned multiple-replica execution. The output vector points map
+**output_walkers** can be used to analyze and visualize spatial
+distribution of walkers at different simulation times (note that the
+resulting water depth is based on the density of these walkers).
+Duration of simulation is controlled by the
 **duration** parameter. The default value is 10 minutes, reaching the
 steady-state may require much longer time, depending on the time step,
 complexity of terrain, land cover and size of the area. Output walker,
@@ -218,16 +216,6 @@ d.mon stop=cairo
 *Figure: Simulated water depth map in the rural area of the North
 Carolina sample dataset.*
 
-## ERROR MESSAGES
-
-If the module fails with
-
-```sh
-ERROR: nwalk (7000001) > maxw (7000000)!
-```
-
-then a lower **nwalkers** parameter value has to be selected.
-
 ## REFERENCES
 
 - Mitasova, H., Thaxton, C., Hofierka, J., McLaughlin, R., Moore, A.,
diff --git a/raster/r.sim/simlib/hydro.c b/raster/r.sim/simlib/hydro.c
index 5631c2e6d06..af46e5dc9f1 100644
--- a/raster/r.sim/simlib/hydro.c
+++ b/raster/r.sim/simlib/hydro.c
@@ -45,25 +45,32 @@ void main_loop(const Setup *setup, const Geometry *geometry,
 {
     int i, l, k;
     int iblock;
-    double conn;
+    double conn = 1.0;
     double addac;
 
+    // nblock is reserved for Monte Carlo replicas. A future
+    // change will allow nblock > 1, give each replica an
+    // rwalk-sized walker subset and potentially run replicas concurrently and
+    // reduce into the shared grids after the iblock loop. The factor and
+    // conn formulas below already encode that design (factor's denominator
+    // (rwalk * nblock) equals total walkers; conn = nblock/iblock scales a
+    // sequential cumulative partial sum to an estimator of the eventual
+    // total). With nblock = 1 today, the loop is a no-op wrapper and conn
+    // collapses to 1.0. The historical auto-split based on a static MAXW
+    // cap was removed: it broke per-block walker accounting (rwalk was not
+    // actually divided), produced biased depth/discharge for users.
     int nblock = 1;
 
     double stxm = geometry->stepx * (double)(geometry->mx + 1) - geometry->xmin;
     double stym = geometry->stepy * (double)(geometry->my + 1) - geometry->ymin;
     double deldif = sqrt(setup->deltap) * settings->frac; /* diffuse factor */
 
-    if (sim->maxwa > (MAXW - geometry->mx * geometry->my)) {
-        nblock = 1 + sim->maxwa / (MAXW - geometry->mx * geometry->my);
-        sim->maxwa = sim->maxwa / nblock;
-    }
     double factor =
         setup->deltap * setup->sisum / (sim->rwalk * (double)nblock);
 
     G_debug(2, " deldif, factor %f %e", deldif, factor);
     G_debug(2, " maxwa, nblock %d %d", sim->maxwa, nblock);
-    G_debug(2, "rwalk,sisum: %f %f", sim->rwalk, setup->sisum);
+    G_debug(2, "rwalk, sisum: %f %f", sim->rwalk, setup->sisum);
 
     for (iblock = 1; iblock <= nblock; iblock++) {
         int lw = 0;
@@ -101,12 +108,16 @@ void main_loop(const Setup *setup, const Geometry *geometry,
             }
         }
         sim->nwalk = lw;
-        G_debug(2, " nwalk, maxw %d %d", sim->nwalk, MAXW);
+        G_debug(2, " nwalk %d", sim->nwalk);
         G_debug(2, " walkwe (walk weight),frac %f %f", walkwe, settings->frac);
 
         sim->nwalka = 0;
         int nwalka = 0;
 
+        // conn scales the cumulative partial sum in gama into an estimator
+        // of the eventual total when blocks run sequentially.
+        conn = (double)nblock / (double)iblock;
+
         /* ********************************************************** */
         /*       main loop over the projection time */
         /* *********************************************************** */
@@ -131,7 +142,6 @@ void main_loop(const Setup *setup, const Geometry *geometry,
             /* ************************************************************ */
 
             addac = factor;
-            conn = (double)nblock / (double)iblock;
             if (i == 1) {
                 addac = factor * .5;
             }
@@ -330,7 +340,6 @@ void main_loop(const Setup *setup, const Geometry *geometry,
                     erod(grids->gama, setup, geometry,
                          grids); /* divergence of gama field */
 
-                conn = (double)nblock / (double)iblock;
                 int itime = (int)(i * setup->deltap * setup->timec);
                 int ii = output_data(itime, conn, setup, geometry, settings,
                                      sim, inputs, outputs, grids);
@@ -385,12 +394,18 @@ void main_loop(const Setup *setup, const Geometry *geometry,
            }
            } */
 
+        // Per-block sample for the Monte Carlo standard-deviation estimator
+        // over nblock replicas, matching the original Fortran: accumulate
+        // (gama * conn)^2 here, then finalize as sqrt(|gammas/nblock - gama^2|)
+        // after the iblock loop closes. With nblock = 1 there is only one
+        // sample and the finalized value is zero everywhere; the map becomes
+        // meaningful once nblocks > 1 will be allowed.
         if (outputs->err != NULL) {
             for (k = 0; k < geometry->my; k++) {
                 for (l = 0; l < geometry->mx; l++) {
                     if (grids->zz[k][l] != UNDEF) {
                         double d1 = grids->gama[k][l] * (double)conn;
-                        grids->gammas[k][l] += pow(d1, 3. / 5.);
+                        grids->gammas[k][l] += d1 * d1;
                     } /* DEFined area */
                 }
             }
@@ -400,9 +415,26 @@ void main_loop(const Setup *setup, const Geometry *geometry,
     }
     /*                       ........ end of iblock loop */
 
+    // Finalize the err map as the sample standard deviation of the per-block
+    // estimators of the final field: sqrt(|E[X^2] - E[X]^2|), where each X
+    // is gama * (nblock/iblock) recorded at the end of block iblock.
+    if (outputs->err != NULL) {
+        for (k = 0; k < geometry->my; k++) {
+            for (l = 0; l < geometry->mx; l++) {
+                if (grids->zz[k][l] != UNDEF) {
+                    double mean_sq = grids->gama[k][l] * grids->gama[k][l];
+                    double mean_of_sq = grids->gammas[k][l] / (double)nblock;
+                    grids->gammas[k][l] = sqrt(fabs(mean_of_sq - mean_sq));
+                }
+            }
+        }
+    }
+
     /* Write final maps here because we know the last time stamp here */
     if (!settings->ts) {
-        conn = (double)nblock / (double)iblock;
+        // All blocks have completed; gama is the eventual cumulative total,
+        // so no extrapolation is needed.
+        conn = 1.0;
         int itime = (int)(i * setup->deltap * setup->timec);
         int ii = output_data(itime, conn, setup, geometry, settings, sim,
                              inputs, outputs, grids);
diff --git a/raster/r.sim/simlib/output.c b/raster/r.sim/simlib/output.c
index 7f8c07485eb..fb858c12357 100644
--- a/raster/r.sim/simlib/output.c
+++ b/raster/r.sim/simlib/output.c
@@ -94,9 +94,12 @@ void output_walker_as_vector(int tt_minutes, int ndigit,
     return;
 }
 
-/* Soeren 8. Mar 2011 TODO:
- * This function needs to be refractured and splittet into smaller parts */
-int output_data(int tt, double ft G_UNUSED, const Setup *setup,
+/* conn is the sequential-block extrapolation factor nblock/iblock: scales
+ * the cumulative partial sum in gama into an estimator of the eventual total
+ * for snapshots taken before all blocks have run. With nblock = 1 (or after
+ * the final block) conn = 1.0 and the output formulas are unchanged. err is
+ * written from gammas, which already has conn baked into its accumulator. */
+int output_data(int tt, double conn, const Setup *setup,
                 const Geometry *geometry, const Settings *settings,
                 const Simulation *sim, const Inputs *inputs,
                 const Outputs *outputs, const Grids *grids)
@@ -234,7 +237,7 @@ int output_data(int tt, double ft G_UNUSED, const Setup *setup,
                 if (grids->zz[i][j] == UNDEF || grids->gama[i][j] == UNDEF)
                     Rast_set_f_null_value(depth_cell + j, 1);
                 else {
-                    a1 = pow(grids->gama[i][j], 3. / 5.);
+                    a1 = pow(grids->gama[i][j] * conn, 3. / 5.);
                     depth_cell[j] = (FCELL)a1; /* add conv? */
                     gmax = amax1(gmax, a1);
                 }
@@ -248,7 +251,7 @@ int output_data(int tt, double ft G_UNUSED, const Setup *setup,
                     grids->cchez[i][j] == UNDEF)
                     Rast_set_f_null_value(disch_cell + j, 1);
                 else {
-                    a2 = geometry->step * grids->gama[i][j] *
+                    a2 = geometry->step * grids->gama[i][j] * conn *
                          grids->cchez[i][j];   /* cchez incl. sqrt(sinsl) */
                     disch_cell[j] = (FCELL)a2; /* add conv? */
                     dismax = amax1(dismax, a2);
@@ -274,7 +277,7 @@ int output_data(int tt, double ft G_UNUSED, const Setup *setup,
                 if (grids->zz[i][j] == UNDEF || grids->gama[i][j] == UNDEF)
                     Rast_set_f_null_value(conc_cell + j, 1);
                 else {
-                    conc_cell[j] = (FCELL)grids->gama[i][j];
+                    conc_cell[j] = (FCELL)(grids->gama[i][j] * conn);
                     /*      gsmax = amax1(gsmax, gama[i][j]); */
                 }
             }
@@ -287,7 +290,7 @@ int output_data(int tt, double ft G_UNUSED, const Setup *setup,
                     grids->slope[i][j] == UNDEF)
                     Rast_set_f_null_value(flux_cell + j, 1);
                 else {
-                    a2 = grids->gama[i][j] * grids->slope[i][j];
+                    a2 = grids->gama[i][j] * conn * grids->slope[i][j];
                     flux_cell[j] = (FCELL)a2;
                     dismax = amax1(dismax, a2);
                 }
diff --git a/raster/r.sim/simlib/simlib.h b/raster/r.sim/simlib/simlib.h
index d516fd1b70d..edd2e7302e8 100644
--- a/raster/r.sim/simlib/simlib.h
+++ b/raster/r.sim/simlib/simlib.h
@@ -8,7 +8,6 @@
 #include <stdio.h>
 
 #define EPS         1.e-7
-#define MAXW        7000000
 #define UNDEF       -9999
 
 #define NUM_THREADS "1"
@@ -57,7 +56,7 @@ typedef struct {
     int nwalka;            // Remaining walkers in an iteration
     int nstack;            // Number of output walkers
     struct point3D *stack; // Output 3D walkers
-    int maxwa;             // Number of input walkers per block
+    int maxwa;             // Number of total walkers
     double rwalk;      // Number of input walkers per block as double precision
     struct point3D *w; // Weight of walkers
     struct point2D *vavg; // Average velocity of walkers
@@ -149,7 +148,7 @@ void main_loop(const Setup *setup, const Geometry *geometry,
                const Settings *settings, Simulation *sim,
                ObservationPoints *points, const Inputs *inputs,
                const Outputs *outputs, Grids *grids);
-int output_data(int, double, const Setup *setup, const Geometry *geometry,
+int output_data(int, double conn, const Setup *setup, const Geometry *geometry,
                 const Settings *settings, const Simulation *sim,
                 const Inputs *inputs, const Outputs *outputs,
                 const Grids *grids);