diff --git a/src/raster/r.hydro.flatten/r.hydro.flatten.html b/src/raster/r.hydro.flatten/r.hydro.flatten.html
index f43a043c1b..c0c8cc76be 100644
--- a/src/raster/r.hydro.flatten/r.hydro.flatten.html
+++ b/src/raster/r.hydro.flatten/r.hydro.flatten.html
@@ -15,6 +15,14 @@ <h2>DESCRIPTION</h2>
 output representing the standard deviation of the lidar elevation values
 along the water body's edge. Higher deviation suggests problematic areas
 that need to be further inspected.
+The optional output <b>filled_elevation</b> is a raster map of the input ground surface
+filled with the computed <b>water_elevation</b> raster map.
+
+<p>
+The <b>breaklines</b> parameter is an optional input that specifies a vector map of lines
+that represent e.g., a break between an impoundment and downstream river, allowing
+correct elevation computation.
+</p>
 
 <p>
 To keep the intermediate results for inspection, use flag <b>-k</b>.
@@ -44,7 +52,7 @@ <h2>EXAMPLE</h2>
 # convert elevation from feet to meters
 r.mapcalc "ground_m = ground * 0.304800609601219"
 # derive elevation of water bodies and standard deviation
-r.hydro.flatten input=ground_m water_elevation=water_elevation water_elevation_stddev=water_elevation_stddev percentile=10 misize=4000
+r.hydro.flatten input=ground_m water_elevation=water_elevation water_elevation_stddev=water_elevation_stddev filled_elevation=filled percentile=10 misize=4000
 </pre></div>
 
 <div align="center" style="margin: 10px">
diff --git a/src/raster/r.hydro.flatten/r.hydro.flatten.py b/src/raster/r.hydro.flatten/r.hydro.flatten.py
index 56d06096e3..8350934a25 100644
--- a/src/raster/r.hydro.flatten/r.hydro.flatten.py
+++ b/src/raster/r.hydro.flatten/r.hydro.flatten.py
@@ -28,6 +28,11 @@
 # % key: input
 # % description: Raster map of binned lidar point elevation
 # %end
+# %option G_OPT_V_INPUT
+# % key: breaklines
+# % description: Vector map of breaklines
+# % required: no
+# %end
 # %option G_OPT_R_OUTPUT
 # % key: water_elevation
 # % description: Represents single elevation value for each water body
@@ -37,6 +42,11 @@
 # % key: water_elevation_stddev
 # % description: Raster map of derived water elevation standard deviation
 # %end
+# %option G_OPT_R_OUTPUT
+# % key: filled_elevation
+# % required: no
+# % description: Raster map representing filled digital elevation model
+# %end
 # %option
 # % key: percentile
 # % type: double
@@ -92,6 +102,7 @@ def get_name(basename):
             return name
 
     ground = options["input"]
+    breaklines = options["breaklines"]
     size_threshold = options["min_size"]
     if size_threshold:
         size_threshold = int(size_threshold)
@@ -107,6 +118,17 @@ def get_name(basename):
     region_m = gs.parse_command("g.region", flags="gm")
     resolution_m = (float(region_m["nsres"]) + float(region_m["ewres"])) / 2
     tmp_rfillstats = get_name("rfillstats")
+
+    if breaklines:
+        tmp_breaklines = get_name("breaklines")
+        gs.run_command(
+            "v.to.rast",
+            input=breaklines,
+            output=tmp_breaklines,
+            use="val",
+            flags="d",
+            value=1000,
+        )
     gs.run_command(
         "r.fill.stats",
         flags="k",
@@ -125,6 +147,15 @@ def get_name(basename):
         distances=[x * resolution_m for x in range(1, buffer_last_strip)],
         units="meters",
     )
+    if breaklines:
+        tmp_buffer_with_breaklines = get_name("buffer_with_breaklines")
+        gs.run_command(
+            "r.patch",
+            input=[tmp_breaklines, tmp_buffer],
+            output=tmp_buffer_with_breaklines,
+        )
+        tmp_buffer = tmp_buffer_with_breaklines
+
     tmp_reclass_for_clump = get_name("reclass_for_clump")
     gs.write_command(
         "r.reclass",
@@ -155,30 +186,53 @@ def get_name(basename):
         method="stddev",
         output=tmp_water_stddev,
     )
-    tmp_water_elevation_dist = get_name("water_elevation_dist")
+    tmp_water_elevation_zonal = get_name("water_elevation_zonal")
     gs.run_command(
-        "r.grow.distance", input=tmp_water_elevation, value=tmp_water_elevation_dist
+        "r.stats.zonal",
+        base=tmp_clump,
+        cover=tmp_water_elevation,
+        method="average",
+        output=tmp_water_elevation_zonal,
     )
-    tmp_water_stddev_dist = get_name("water_stddev_dist")
+    tmp_water_elevation_stddev_zonal = get_name("water_elevation_stddev_zonal")
     gs.run_command(
-        "r.grow.distance", input=tmp_water_stddev, value=tmp_water_stddev_dist
-    )
-    tmp_water_elevation_dist_res = get_name("water_elevation_dist_res")
-    gs.mapcalc(
-        f"{tmp_water_elevation_dist_res} = if ({tmp_buffer} < {buffer_last_strip}, "
-        f"{tmp_water_elevation_dist}, null())"
+        "r.stats.zonal",
+        base=tmp_clump,
+        cover=tmp_water_stddev,
+        method="average",
+        output=tmp_water_elevation_stddev_zonal,
     )
-    tmp_water_stddev_dist_res = get_name("water_stddev_dist_res")
+    tmp_water_elevation_zonal_res = get_name("water_elevation_zonal_res")
+    if breaklines:
+        water_elevation_zonal_res_breaklines = get_name(
+            "water_elevation_zonal_res_breaklines"
+        )
+        gs.mapcalc(
+            f"{water_elevation_zonal_res_breaklines} = if (isnull({tmp_strip}), {tmp_water_elevation_zonal}, null())"
+        )
+        # heal the breakline holes
+        gs.run_command(
+            "r.neighbors",
+            input=water_elevation_zonal_res_breaklines,
+            selection=tmp_breaklines,
+            output=tmp_water_elevation_zonal_res,
+            size=5,
+        )
+    else:
+        gs.mapcalc(
+            f"{tmp_water_elevation_zonal_res} = if (isnull({tmp_strip}), {tmp_water_elevation_zonal}, null())"
+        )
+
+    tmp_water_elevation_stddev_zonal_res = get_name("water_elevation_stddev_zonal_res")
     gs.mapcalc(
-        f"{tmp_water_stddev_dist_res} = if ({tmp_buffer} < {buffer_last_strip}, "
-        f"{tmp_water_stddev_dist}, null())"
+        f"{tmp_water_elevation_stddev_zonal_res} = if (isnull({tmp_strip}), {tmp_water_elevation_stddev_zonal}, null())"
     )
     if size_threshold:
         size_threshold /= region["nsres"] * region["ewres"]
         tmp_reclass = get_name("reclass")
         gs.write_command(
             "r.reclass",
-            input=tmp_water_elevation_dist_res,
+            input=tmp_water_elevation_zonal_res,
             output=tmp_reclass,
             rules="-",
             stdin="* = 1",
@@ -194,18 +248,28 @@ def get_name(basename):
             output=tmp_size,
         )
         gs.mapcalc(
-            f"{options['water_elevation']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_dist_res}, null())"
+            f"{options['water_elevation']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_zonal_res}, null())"
         )
         gs.mapcalc(
-            f"{options['water_elevation_stddev']} = if ({tmp_size} > {size_threshold}, {tmp_water_stddev_dist_res}, null())"
+            f"{options['water_elevation_stddev']} = if ({tmp_size} > {size_threshold}, {tmp_water_elevation_stddev_zonal_res}, null())"
         )
     else:
-        gs.mapcalc(f"{options['water_elevation']} = {tmp_water_elevation_dist_res}")
-        gs.mapcalc(f"{options['water_elevation_stddev']} = {tmp_water_stddev_dist_res}")
+        gs.mapcalc(f"{options['water_elevation']} = {tmp_water_elevation_zonal_res}")
+        gs.mapcalc(
+            f"{options['water_elevation_stddev']} = {tmp_water_elevation_stddev_zonal_res}"
+        )
     gs.run_command("r.colors", map=options["water_elevation"], raster=ground)
     gs.run_command("r.colors", map=options["water_elevation_stddev"], color="reds")
     gs.raster_history(options["water_elevation"])
     gs.raster_history(options["water_elevation_stddev"])
+    if options["filled_elevation"]:
+        gs.run_command(
+            "r.patch",
+            input=[tmp_rfillstats, options["water_elevation"]],
+            output=options["filled_elevation"],
+        )
+        gs.run_command("r.colors", map=options["filled_elevation"], raster=ground)
+        gs.raster_history(options["filled_elevation"])
 
 
 if __name__ == "__main__":