diff --git a/worldeditadditions/lib/erode/erode.lua b/worldeditadditions/lib/erode/erode.lua
index 6184e6a..76bd0e5 100644
--- a/worldeditadditions/lib/erode/erode.lua
+++ b/worldeditadditions/lib/erode/erode.lua
@@ -1,6 +1,7 @@
 worldeditadditions.erode = {}
 
 dofile(worldeditadditions.modpath.."/lib/erode/snowballs.lua")
+dofile(worldeditadditions.modpath.."/lib/erode/river.lua")
 
 
 function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
diff --git a/worldeditadditions/lib/erode/river.lua b/worldeditadditions/lib/erode/river.lua
index 790d005..b6ef753 100644
--- a/worldeditadditions/lib/erode/river.lua
+++ b/worldeditadditions/lib/erode/river.lua
@@ -14,8 +14,10 @@ end
 function worldeditadditions.erode.river(heightmap_initial, heightmap, heightmap_size, region_height, params_custom)
 	local params = {
 		steps = 1,				-- Number of rounds/passes of the algorithm to run
-		remove_sides = "0,1",	-- Cells with this many adjacent horizontal neighbours will be removed
-		fill_sides = "4,3"		-- Cells with this many adjaect horizontal neighbours will be filled in
+		remove_sides = "4,3",	-- Cells with this many adjacent horizontal neighbours that are lower than the current pixel will be removed
+		fill_sides = "4,3",		-- Cells with this many adjacect horizontal neighbours that are higher than the current pixel will be filled in
+		doraise = true,			-- Whether to do raise operations or not
+		dolower = true			-- Whether to do lower operations or not
 	}
 	-- Apply the custom settings
 	wea.table_apply(params_custom, params)
@@ -27,64 +29,92 @@ function worldeditadditions.erode.river(heightmap_initial, heightmap, heightmap_
 	local filled = 0
 	local removed = 0
 	for i=1,params.steps do
+		print("[DEBUG:river] step ", i)
+		wea.format.array_2d(heightmap, heightmap_size.x)
 		local time_start = wea.get_ms_time()
 		
+		-- Store up changes to make and make them at the end of the step
+		-- This is important, because decisions 
+		local fill = { }	-- Indexes to add 1 to
+		local remove = { }	-- Indexes to take 1 away from
+		
 		for z = heightmap_size.z - 1, 0, -1 do
 			for x = heightmap_size.x - 1, 0, -1 do
 				local hi = z*heightmap_size.x + x
 				local thisheight = heightmap[hi]
+				print("[DEBUG:river] z", z, "x", x, "thisheight", thisheight)
 				
-				local sides = 0
-				local adjacent_heights = { }
-				if x > 0 then
-					table.insert(adjacent_heights, heightmap[z*heightmap_size.x + x-1])
-					if heightmap[z*heightmap_size.x + x-1] >= thisheight then
-						sides = sides + 1
-					end
-				end
-				if x < heightmap_size.x - 1 then
-					table.insert(adjacent_heights, heightmap[z*heightmap_size.x + x+1])
-					if heightmap[z*heightmap_size.x + x+1] >= thisheight then
-						sides = sides + 1
-					end
-				end
-				if z > 0 then
-					table.insert(adjacent_heights, heightmap[(z-1)*heightmap_size.x + x])
-					if heightmap[(z-1)*heightmap_size.x + x] >= thisheight then
-						sides = sides + 1
-					end
-				end
-				if z < heightmap_size.z - 1 then
-					table.insert(adjacent_heights, heightmap[(z+1)*heightmap_size.x + x])
-					if heightmap[(z+1)*heightmap_size.x + x] >= thisheight then
-						sides = sides + 1
-					end
+				local height_up = heightmap[hi]
+				local height_down = heightmap[hi]
+				local height_left = heightmap[hi]
+				local height_down = heightmap[hi]
+				
+				if x > 0 then height_left = heightmap[z*heightmap_size.x + x-1] end
+				if x < heightmap_size.x - 1 then height_right = heightmap[z*heightmap_size.x + x+1] end
+				if z > 0 then height_up = heightmap[(z-1)*heightmap_size.x + x] end
+				if z < heightmap_size.z - 1 then height_down = heightmap[(z+1)*heightmap_size.x + x] end
+				
+				-- Whether this pixel is on the edge
+				local isedge = x <= 0
+					or z <= 0
+					or x >= heightmap_size.x - 1
+					or z >= heightmap_size.z - 1
+				
+				local sides_higher = 0	-- Number of sides higher than this pixel
+				local sides_lower = 0	-- Number of sides lower than this pixel
+				if not isedge then
+					if height_down > thisheight then sides_higher = sides_higher + 1 end
+					if height_up > thisheight then sides_higher = sides_higher + 1 end
+					if height_left > thisheight then sides_higher = sides_higher + 1 end
+					if height_right > thisheight then sides_higher = sides_higher + 1 end
+					
+					if height_down < thisheight then sides_lower = sides_lower + 1 end
+					if height_up < thisheight then sides_lower = sides_lower + 1 end
+					if height_left < thisheight then sides_lower = sides_lower + 1 end
+					if height_right < thisheight then sides_lower = sides_lower + 1 end
 				end
 				
+				-- Perform an action, but only if we're not on the edge
+				-- This is important, as we can't accurately tell how many
+				-- adjacent neighbours a pixel on the edge has.
 				local action = "none"
-				for i,sidecount in ipairs(params.fill_sides) do
-					if sidecount == sides then
-						action = "fill"
-						break
-					end
-				end
-				for i,sidecount in ipairs(params.remove_sides) do
-					if sidecount == sides then
-						action = "remove"
-						break
+				if not isedge then 
+					if sides_higher > sides_lower then
+						for i,sidecount in ipairs(params.fill_sides) do
+							if sidecount == sides_higher then
+								action = "fill"
+								break
+							end
+						end
+					else
+						for i,sidecount in ipairs(params.remove_sides) do
+							if sidecount == sides_lower then
+								action = "remove"
+								break
+							end
+						end
 					end
 				end
 				
-				if action == "fill" then
-					heightmap[hi] = heightmap[hi] + 1
+				print("[DEBUG:river] sides_higher", sides_higher, "sides_lower", sides_lower, "action", action)
+				if action == "fill" and params.doraise then
+					table.insert(fill, hi)
 					filled = filled + 1
-				elseif action == "remove" then
-					heightmap[hi] = heightmap[hi] - 1
+				elseif action == "remove" and params.dolower then
+					table.insert(remove, hi)
 					removed = removed + 1
 				end
+				wea.format.array_2d(heightmap, heightmap_size.x)
 			end
 		end
 		
+		for i,hi in ipairs(fill) do
+			heightmap[hi] = heightmap[hi] + 1
+		end
+		for i,hi in ipairs(remove) do
+			heightmap[hi] = heightmap[hi] - 1
+		end
+		
 		table.insert(timings, wea.get_ms_time() - time_start)
 	end
 	
diff --git a/worldeditadditions/utils/format/array_2d.lua b/worldeditadditions/utils/format/array_2d.lua
index 8af7a8a..5e940d1 100644
--- a/worldeditadditions/utils/format/array_2d.lua
+++ b/worldeditadditions/utils/format/array_2d.lua
@@ -3,7 +3,7 @@
 -- @param	tbl		number[]	The ZERO-indexed list of numbers
 -- @param	width	number		The width of 2D array.
 function worldeditadditions.format.array_2d(tbl, width)
-	print("==== count: "..#tbl..", width:"..width.." ====")
+	print("==== count: "..(#tbl+1)..", width:"..width.." ====")
 	local display_width = 1
 	for _i,value in pairs(tbl) do
 		display_width = math.max(display_width, #tostring(value))
diff --git a/worldeditadditions_commands/commands/erode.lua b/worldeditadditions_commands/commands/erode.lua
index 807d2ce..d41df3e 100644
--- a/worldeditadditions_commands/commands/erode.lua
+++ b/worldeditadditions_commands/commands/erode.lua
@@ -4,7 +4,7 @@
 -- ██      ██   ██ ██    ██ ██   ██ ██
 -- ███████ ██   ██  ██████  ██████  ███████
 worldedit.register_command("erode", {
-	params = "[<snowballs|...> [<key_1> [<value_1>]] [<key_2> [<value_2>]] ...]",
+	params = "[<snowballs|river> [<key_1> [<value_1>]] [<key_2> [<value_2>]] ...]",
 	description = "**experimental** Runs the specified erosion algorithm over the given defined region. This may occur in 2d or 3d. Currently implemented algorithms: snowballs (default;2d hydraulic-like). Also optionally takes an arbitrary set of key - value pairs representing parameters to pass to the algorithm. See the full documentation for details.",
 	privs = { worldedit = true },
 	require_pos = 2,