diff --git a/sote/game/scenes/world-gen.lua b/sote/game/scenes/world-gen.lua
index 238a22ff..10bafd07 100644
--- a/sote/game/scenes/world-gen.lua
+++ b/sote/game/scenes/world-gen.lua
@@ -23,31 +23,25 @@ local function map_tiles_to_hex()
 		local lat, lon = tile:latlon()
 		local q, r, face = hex.latlon_to_hex_coords(lat, lon, wg.world.size)
 
-		wg.world:cache_tile_coord(tile.tile_id, q, r, face)
+		wg.world:cache_square_ti_by_hex_coords(q, r, face)
 	end
 end
 
 local function load_mapping_from_file(file)
 	for row in require("game.file-utils").csv_rows(file) do
-		local tile_id = tonumber(row[1])
-		local q = tonumber(row[2])
-		local r = tonumber(row[3])
-		local face = tonumber(row[4])
-
-		wg.world:cache_tile_coord(tile_id, q, r, face)
+		wg.world:cache_square_ti_by_hex_ti(tonumber(row[2]))
 	end
 end
 
-local function cache_tile_coord()
+local function cache_cube_world_tiles()
 	print("Caching tile coordinates...")
 
 	if debug.map_tiles_from_file then
 		-- it's faster to load the pre-calculated coordinates from a file than to calculate them on the fly
-		load_mapping_from_file("d:\\temp\\hex_mapping.csv")
+		load_mapping_from_file("d:\\temp\\hex_mapping2.csv")
 	else
 		map_tiles_to_hex()
 	end
-	-- wg.world:map_hex_coords()
 
 	print("Done caching tile coordinates")
 end
@@ -119,7 +113,7 @@ function wg.generate_coro()
 			coroutine.yield()
 
 			wg.message = "Caching tile coordinates"
-			prof.run_with_profiling(function() cache_tile_coord() end, "[scenes.world-gen]", "cache_tile_coord")
+			prof.run_with_profiling(function() cache_cube_world_tiles() end, "[scenes.world-gen]", "cache_cube_world_tiles")
 			coroutine.yield()
 		end
 
diff --git a/sote/libsote/debug-control-panel.lua b/sote/libsote/debug-control-panel.lua
index 342343b7..4a496a15 100644
--- a/sote/libsote/debug-control-panel.lua
+++ b/sote/libsote/debug-control-panel.lua
@@ -1,8 +1,9 @@
 local dcp = {}
 
-dcp.align_to_sote_coords = false -- this will align hex world storage to match the order from original sote, very useful when debugging/validating port
-dcp.map_tiles_from_file = false -- this will load cube world tile IDs mapping to hex coordinates from a file, as it's faster than computing them from lat lon
+dcp.map_tiles_from_file = false -- this will load cube world tile IDs mapping to hex from a file, as it's faster than computing them from lat lon
 dcp.use_sote_climate_data = false -- climate model was ported from sote, but with some changes; this will enable import of original sote climate data from a csv file, to aid in debugging/validating port
+dcp.use_sote_ice_data = false -- this will enable import of original sote ice data from a csv file, to aid in debugging/validating port
+dcp.align_to_sote_coords = dcp.use_sote_climate_data or dcp.use_sote_ice_data -- this will align hex world storage to match the order from original sote, very useful when debugging/validating port
 
 dcp.save_maps = false -- this will export maps to PNG
 dcp.maps_selection = {
@@ -11,7 +12,9 @@ dcp.maps_selection = {
 	climate = false,
 	waterflow = false,
 	waterbodies = false,
-	debug = false
+	watersheds = false,
+	debug1 = false,
+	debug2 = false
 }
 
 -- seed = 58738 -- climate integration was done on this one
@@ -20,7 +23,7 @@ dcp.maps_selection = {
 -- seed = 6618 -- tiny islands?
 -- seed = 49597 -- interesting looking one, huge northern ice cap (with lua climate model)
 -- seed = 91170 -- huge lake
-dcp.fixed_seed = nil
+-- dcp.fixed_seed = nil
 -- dcp.fixed_seed = 12177
 
 return dcp
\ No newline at end of file
diff --git a/sote/libsote/debug-loggers.lua b/sote/libsote/debug-loggers.lua
index bb8a6ce0..5d4db8a2 100644
--- a/sote/libsote/debug-loggers.lua
+++ b/sote/libsote/debug-loggers.lua
@@ -44,16 +44,6 @@ function logger:log(message, do_flush)
 	end
 end
 
-local loggers = {}
-
-local latlon_logger = nil
-local neighbors_logger = nil
-local waterflow_logger = nil
-local parent_material_logger = nil
-local glacial_logger = nil
-local climate_logger = nil
-local lakes_logger = nil
-
 local function get_logger(logger_instance, logname, path, unique)
 	if logger_instance == nil then
 		logger_instance = logger:new(logname, path, unique)
@@ -62,32 +52,46 @@ local function get_logger(logger_instance, logname, path, unique)
 	return logger_instance
 end
 
+local loggers = {}
+
+local latlon_logger = nil
 function loggers.get_latlon_logger(path)
 	return get_logger(latlon_logger, "latlon", path)
 end
 
+local neighbors_logger = nil
 function loggers.get_neighbors_logger(path)
 	return get_logger(neighbors_logger, "neighbours", path)
 end
 
+local waterflow_logger = nil
 function loggers.get_waterflow_logger(path)
 	return get_logger(waterflow_logger, "waterflow", path)
 end
 
+local parent_material_logger = nil
 function loggers.get_parent_material_logger(path)
 	return get_logger(parent_material_logger, "parent_material", path)
 end
 
+local glacial_logger = nil
 function loggers.get_glacial_logger(path)
 	return get_logger(glacial_logger, "glacial", path)
 end
 
+local climate_logger = nil
 function loggers.get_climate_logger(path)
 	return get_logger(climate_logger, "climate", path)
 end
 
+local lakes_logger = nil
 function loggers.get_lakes_logger(path)
 	return get_logger(lakes_logger, "lakes", path)
 end
 
+local rivers_logger = nil
+function loggers.get_rivers_logger(path)
+	return get_logger(rivers_logger, "rivers", path)
+end
+
 return loggers
diff --git a/sote/libsote/glaciation/glacial-formation.lua b/sote/libsote/glaciation/glacial-formation.lua
index 4ee3b741..ea0367ea 100644
--- a/sote/libsote/glaciation/glacial-formation.lua
+++ b/sote/libsote/glaciation/glacial-formation.lua
@@ -731,6 +731,7 @@ function gf.run(world_obj)
 
 	if enable_debug then
 		world:adjust_debug_channels(2)
+		world:reset_debug_all()
 	end
 
 	world:fill_ffi_array(glacial_seed, false)
diff --git a/sote/libsote/heap-sort.lua b/sote/libsote/heap-sort.lua
index a4ae6354..d23c1005 100644
--- a/sote/libsote/heap-sort.lua
+++ b/sote/libsote/heap-sort.lua
@@ -64,10 +64,10 @@ end
 local ffi = require("ffi")
 
 ---@param get_primary fun(i:number):any
----@param get_secondary fun(i:number):any
+---@param get_secondary nil|fun(i:number):any
 ---@param n number
 ---@param desc_primary boolean
----@param desc_secondary boolean
+---@param desc_secondary boolean|nil
 function hs.heap_sort_indices(get_primary, get_secondary, n, desc_primary, desc_secondary)
 	desc_primary = desc_primary or false
 	desc_secondary = desc_secondary or false
diff --git a/sote/libsote/hydrology/calculate-waterflow.lua b/sote/libsote/hydrology/calculate-waterflow.lua
index e380ba09..6be663a7 100644
--- a/sote/libsote/hydrology/calculate-waterflow.lua
+++ b/sote/libsote/hydrology/calculate-waterflow.lua
@@ -19,8 +19,6 @@ end
 
 local function clear_current_elevation_on_lakes(world)
 	world:for_each_waterbody(function(wb)
-		if not wb:is_valid() then return end
-
 		if wb.type == wb.TYPES.freshwater_lake or wb.type == wb.TYPES.saltwater_lake then
 			wb.tmp_float_1 = 0
 		end
diff --git a/sote/libsote/hydrology/gen-dynamic-lakes.lua b/sote/libsote/hydrology/gen-dynamic-lakes.lua
index 00c45b5b..655c0f6e 100644
--- a/sote/libsote/hydrology/gen-dynamic-lakes.lua
+++ b/sote/libsote/hydrology/gen-dynamic-lakes.lua
@@ -2,12 +2,13 @@ local dl = {}
 
 local world
 
+local waterbody = require "libsote.hydrology.waterbody"
 local open_issues = require "libsote.hydrology.open-issues"
 
 -- local logger = require("libsote.debug-loggers").get_lakes_logger("d:/temp")
+
 local prof = require "libsote.profiling-helper"
 local prof_prefix = "[gen-dynamic-lakes]"
-
 local function run_with_profiling(func, log_txt)
 	prof.run_with_profiling(func, prof_prefix, log_txt)
 end
@@ -112,7 +113,7 @@ local function water_flow_from_tile_to_tile()
 		if world:get_waterbody_by_tile(ti) then goto continue1 end
 
 		--* If elevation difference is 0, it can be inferred that we have no tiles lower than the target tile, therefore it should construct a lake.
-		local new_wb = world:create_new_waterbody_from_tile(ti)
+		local new_wb = world:create_waterbody_from_tile(ti, waterbody.TYPES.saltwater_lake)
 
 		--* Mark the tile as water
 		world.is_land[ti] = false
@@ -125,7 +126,6 @@ local function water_flow_from_tile_to_tile()
 		new_wb:set_lowest_shore_tile(world)
 		new_wb.tmp_float_1 = new_wb.tmp_float_1 + water_to_give
 		new_wb.water_level = true_elevation_for_waterflow
-		new_wb.type = new_wb.TYPES.saltwater_lake
 
 		-- logger:log("\tlake " .. new_wb.id .. " created at " .. ti)
 
@@ -139,7 +139,7 @@ local function add_lowest_shore_tile_to_waterbody(wb, lsti)
 	wb.tmp_float_1 = wb.tmp_float_1 + world.tmp_float_2[lsti] / lake_divisor --* If a tile gets eaten by a lake, we automatically move any water that was in the tile to the lake
 	world.tmp_float_2[lsti] = 0
 
-	world:add_tile_to_waterbody(wb, lsti)
+	world:add_tile_to_waterbody(lsti, wb)
 
 	local true_elevation_for_waterflow = world:true_elevation_for_waterflow(lsti)
 	world:for_each_neighbor(lsti, function(nti)
@@ -210,13 +210,12 @@ local function manage_expansion_and_drainage(wb, water_to_disburse)
 
 	--* Drain lake into shore tile with low neighbor that is not the same waterbody
 	if has_lower_neigh then
-		-- local log_str = "\t\t\tlake " .. wb.id .. " draining into tile " .. lowest_shore_ti
 		local lsti_wb = world:get_waterbody_by_tile(lowest_shore_ti)
-		if lsti_wb then
-			-- logger:log(log_str .. " which belongs to lake " .. lsti_wb.id)
+		-- if lsti_wb then
+		-- 	logger:log("\t\t\tlake " .. wb.id .. " draining into tile " .. lowest_shore_ti .. " which belongs to lake " .. lsti_wb.id)
 		-- else
-		-- 	logger:log(log_str)
-		end
+		-- 	logger:log("\t\t\tlake " .. wb.id .. " draining into tile " .. lowest_shore_ti)
+		-- end
 
 		wb.lake_open = true
 		wb.type = wb.TYPES.freshwater_lake
@@ -241,8 +240,6 @@ local function resize_lakes()
 	--* Loop through all waterbodies. Check for active, and check for tempWater.
 
 	world:for_each_waterbody(function(wb)
-		-- if not wb:is_valid() then return end -- do we really need this check?
-
 		if wb.tmp_float_1 <= 0 or wb.type == wb.TYPES.ocean then return end
 		--* Only resize lakes and seas
 
@@ -333,14 +330,16 @@ function dl.run(world_obj)
 	run_with_profiling(function() water_flow_phase() end, "water_flow_phase")
 
 	world:for_each_waterbody(function(wb)
-		if not wb:is_valid() then return end
 		-- logger:log("lake " .. wb.id .. " (" .. wb:size() .. ", " .. wb.water_level .. ")")
-		for _, ti in ipairs(wb.tiles) do
+
+		wb:for_each_tile(function(ti)
 			world.is_land[ti] = false
-		end
+		end)
 	end)
 end
 
+return dl
+
 --* River Plan ///
 --* Possibly generate wetlands first, so that when rivers flow through them, they can then have an "out tile" similar to a lake. The difference is that you don't
 --* have "standing water" like a lake. Instead, its more like a broad, slow moving, shallow river.
@@ -359,5 +358,3 @@ end
 --* Check all endoric waterbodies, including oceans, seas, and saltwater lakes
 --* Check shoreline of those waterbodies and check for tiles with watermovement reaching a particular threshold.
 --* Build a list and follow the waterflow backward from low elevation to high. Stop once the river forks.
-
-return dl
\ No newline at end of file
diff --git a/sote/libsote/hydrology/gen-initial-waterbodies.lua b/sote/libsote/hydrology/gen-initial-waterbodies.lua
index a2aadf7a..cc1f6432 100644
--- a/sote/libsote/hydrology/gen-initial-waterbodies.lua
+++ b/sote/libsote/hydrology/gen-initial-waterbodies.lua
@@ -1,5 +1,6 @@
 local giw = {}
 
+local waterbody = require "libsote.hydrology.waterbody"
 local queue = require("engine.queue"):new()
 
 local waterbodies_created = 0
@@ -11,7 +12,7 @@ local function process(tile_index, world)
 
 	-- "no ice" check is skipped for now
 
-	local new_wb = world:create_new_waterbody_from_tile(tile_index)
+	local new_wb = world:create_waterbody_from_tile(tile_index, waterbody.TYPES.ocean) -- everything seems to start out as oceans?
 	waterbodies_created = waterbodies_created + 1
 
 	queue:enqueue(tile_index)
@@ -22,7 +23,7 @@ local function process(tile_index, world)
 		world:for_each_neighbor(ti, function(nti)
 			if world.is_land[nti] or world:is_tile_waterbody_valid(nti) then return end
 
-			world:add_tile_to_waterbody(new_wb, nti)
+			world:add_tile_to_waterbody(nti, new_wb)
 
 			queue:enqueue(nti)
 		end)
diff --git a/sote/libsote/hydrology/gen-rivers.lua b/sote/libsote/hydrology/gen-rivers.lua
new file mode 100644
index 00000000..173a5a12
--- /dev/null
+++ b/sote/libsote/hydrology/gen-rivers.lua
@@ -0,0 +1,544 @@
+local rivers = {}
+
+local world
+
+local waterbody = require "libsote.hydrology.waterbody"
+
+local enable_debug = false
+-- local logger = require("libsote.debug-loggers").get_rivers_logger("d:/temp")
+
+local prof = require "libsote.profiling-helper"
+local prof_prefix = "[gen-dynamic-lakes]"
+local function run_with_profiling(func, log_txt)
+	prof.run_with_profiling(func, prof_prefix, log_txt)
+end
+
+local initial_candidates = {}
+local sorted_candidates = {}
+
+local stored_bodies = {} -- key: tile index, value: freshwater lake waterbody
+local watershed = {} -- key: tile index, value: river waterbody
+local true_lake
+local true_river
+local fork_count
+
+-- local function set_waterbodies_to_debug(channel)
+-- 	world:for_each_waterbody(function(wb)
+-- 		wb:for_each_tile(function(ti)
+-- 			if wb.type == wb.TYPES.river then
+-- 				world:set_debug_rgba(channel, ti, 173, 216, 230, 255)
+-- 			elseif wb.type == wb.TYPES.freshwater_lake then
+-- 				world:set_debug_rgba(channel, ti, 0, 255, 0, 255)
+-- 			elseif wb.type == wb.TYPES.saltwater_lake then
+-- 				world:set_debug_rgba(channel, ti, 0, 255, 255, 255)
+-- 			elseif wb.type == wb.TYPES.ocean then
+-- 				world:set_debug_rgba(channel, ti, 0, 0, 255, 255)
+-- 			end
+-- 		end)
+-- 	end)
+-- end
+
+local function construct_start_locations()
+	--* Here we are iterating along the coast of each endoreic lake and ocean to find the start tile of rivers
+	world:for_each_waterbody(function(wb)
+		if wb:is_lake_or_ocean() then
+			wb:for_each_tile_in_perimeter(function(ti)
+				if world.water_movement[ti] >= 6000 then
+					table.insert(initial_candidates, ti)
+				end
+			end)
+		end
+
+		--* Setting all tiles inside of a waterbody to 0 watermovement since they are now submerged
+		wb:for_each_tile(function(ti)
+			world.water_movement[ti] = 0
+		end)
+	end)
+end
+
+local function process_drainage_basins(coast_ti)
+	if world:is_tile_waterbody_valid(coast_ti) then return end
+
+	local wb = world:create_waterbody_from_tile(coast_ti, waterbody.TYPES.river)
+
+	local old_layer = {}
+	local new_layer = {}
+
+	table.insert(old_layer, coast_ti)
+	local num_tiles = 1
+
+	while num_tiles > 0 do
+		--* At some point we will need to run a check which evaluated whether we have run into a freshwater lake
+		for _, ti in ipairs(old_layer) do
+			local true_elev = world:true_elevation_for_waterflow(ti)
+
+			world:for_each_neighbor(ti, function(nti)
+				local nwb = world:get_waterbody_by_tile(nti)
+
+				if not world:is_tile_waterbody_valid(nti) then
+					if world.water_movement[nti] > 2000 and world:true_elevation_for_waterflow(nti) > true_elev then
+						world:add_tile_to_waterbody(nti, wb)
+						table.insert(new_layer, nti)
+					end
+				elseif nwb.id ~= wb.id and nwb.type == nwb.TYPES.freshwater_lake then --* Freshwater lakes get tiles IDs, but don't get added to the list of tiles in the drainage basin
+					stored_bodies[nti] = nwb
+					world:reassign_tile_to_waterbody(nti, wb)
+					nwb.basin = wb
+					table.insert(new_layer, nti)
+				else
+					--* ?????
+				end
+			end)
+		end
+
+		old_layer = {}
+		for _, ti in ipairs(new_layer) do
+			table.insert(old_layer, ti)
+		end
+		new_layer = {}
+		num_tiles = #old_layer
+	end
+end
+
+local function construct_drainage_basins()
+	for i = 0, #initial_candidates - 1 do
+		local coast_ti = initial_candidates[sorted_candidates[i] + 1]
+		process_drainage_basins(coast_ti)
+	end
+end
+
+local function find_path_using_waterbodies(ti, wb)
+	local true_elev = world:true_elevation_for_waterflow(ti)
+	local lowest_elev = 100000
+	local lowest_nti = -1
+
+	world:for_each_neighbor(ti, function(nti) --* Here we count candidates and determine who has the lowest elevation
+		local nwb = world:get_waterbody_by_tile(nti)
+		if not nwb or not nwb:is_valid() then return end
+
+		local actual_neigh_elev = world:true_elevation_for_waterflow(nti)
+
+		local is_freshwater_lake_with_standing_water = false
+		local swb = stored_bodies[nti]
+		if swb and swb.water_level > 0 then --* we need to consider the level of the standing water body if we happen to bump into it
+			is_freshwater_lake_with_standing_water = true
+			actual_neigh_elev = swb.water_level
+		end
+
+		if actual_neigh_elev >= true_elev then return end
+
+		if is_freshwater_lake_with_standing_water then
+			lowest_nti = swb.lowest_shore_tile
+		elseif nwb.id == wb.id and world.water_movement[nti] >= 2000 and lowest_elev > actual_neigh_elev then
+			lowest_elev = actual_neigh_elev
+			lowest_nti = nti
+		end
+	end)
+
+	return lowest_nti ~= -1, lowest_nti
+end
+
+local function tag_and_prep_all_tributaries()
+	initial_candidates = {}
+	sorted_candidates = {}
+
+	world:for_each_tile(function(ti)
+		if world.water_movement[ti] < 6000 or world.ice[ti] > 0 then return end
+
+		local wb = world:get_waterbody_by_tile(ti)
+		if not wb or not wb:is_valid() then return end
+
+		local ellibigle_candidate = true
+		world:for_each_neighbor(ti, function(nti) --* Determine elligible candidates for headwaters
+			local nwb = world:get_waterbody_by_tile(nti)
+			if not nwb or not nwb:is_valid() then return end
+
+			if nwb.id == wb.id and world.water_movement[nti] >= 6000 and world:true_elevation_for_waterflow(nti) > world:true_elevation_for_waterflow(ti) then
+				ellibigle_candidate = false
+			end
+		end)
+
+		if ellibigle_candidate then
+			table.insert(initial_candidates, ti)
+		end
+	end)
+
+	--* Now we can just have a second loop which goes through all freshwater lakes with a drain tile and we just add that drain tile to the candidate list
+	--* and allow it to flow toward the ocean like all the others. That should rectify all of our downstream fork counts on all tributaries.
+
+	world:for_each_waterbody(function(wb)
+		if wb.type ~= wb.TYPES.freshwater_lake then return end
+
+		local drain_ti = wb.lowest_shore_tile
+		table.insert(initial_candidates, drain_ti)
+		true_lake[drain_ti] = true
+	end)
+
+	--* Sort candidates by elevation
+	sorted_candidates = require("libsote.heap-sort").heap_sort_indices(
+		function(i) return world:true_elevation_for_waterflow(initial_candidates[i + 1]) end,
+		nil,
+		#initial_candidates,
+		true
+	)
+
+	--* Now we construct a path down from each headwater tributary and we mark every tile as we go down. This method will be used to differentiate different tributaries from one another.
+	for i = 0, #initial_candidates - 1 do
+		local ti = initial_candidates[sorted_candidates[i] + 1]
+
+		local wb = world:get_waterbody_by_tile(ti)
+		if not wb or not wb:is_valid() then goto continue1 end
+
+		--* We need lake drain tiles to be checked as candidates and not over-written
+		if fork_count[ti] > 0 and not true_lake[ti] then goto continue1 end
+
+		local found_path = true
+
+		while found_path do
+			fork_count[ti] = fork_count[ti] + 1
+			found_path, ti = find_path_using_waterbodies(ti, wb)
+		end
+
+		::continue1::
+	end
+end
+
+local function kill_old_basins()
+	world:for_each_waterbody(function(wb) --* Kill old drainage basin rivers and convert their members to a logic variable
+		if wb.type == wb.TYPES.river then --* Kill old rivers
+			wb:for_each_tile(function(ti)
+				watershed[ti] = wb
+			end)
+			world:kill_waterbody(wb)
+		elseif wb.type == wb.TYPES.ocean or wb.type == wb.TYPES.saltwater_lake then --* Prep standing water body variables for next phase
+			wb:for_each_tile(function(ti)
+				fork_count[ti] = 1000000
+				true_river[ti] = true
+			end)
+		else
+			--* ???
+		end
+	end)
+end
+
+local function find_path_using_watershed(ti, wb)
+	local true_elev = world:true_elevation_for_waterflow(ti)
+	local lowest_elevation = 100000
+	local lowest_nti = -1
+
+	world:for_each_neighbor(ti, function(nti) --* Here we count candidates and determine who has the lowest elevation. Lowest elevation neighbor will be next tile in the path
+		local actual_neigh_elev = world:true_elevation_for_waterflow(nti) --* Will function as either the elevation of the tile or the water level of the tile (if it is a lake)
+
+		local is_freshwater_lake_with_standing_water = false
+		local swb = stored_bodies[nti]
+		if swb and swb.water_level > 0 then --* Then we've bumped into a lake and we need to push the water to the drain tile of the lake
+			is_freshwater_lake_with_standing_water = true
+			actual_neigh_elev = swb.water_level
+		end
+
+		if actual_neigh_elev >= true_elev then return end
+
+		local nwb = swb or watershed[nti] or world:get_waterbody_by_tile(nti)
+		if not nwb then return end
+
+		if is_freshwater_lake_with_standing_water then
+			--* We need to terminate expansion and start the next tributary on the drain tile
+			lowest_nti = swb.lowest_shore_tile
+		elseif nwb.type == nwb.TYPES.saltwater_lake or nwb.type == nwb.TYPES.ocean then -- found river end
+			lowest_elevation = actual_neigh_elev
+			lowest_nti = nti
+		elseif nwb.id == wb.id and world.water_movement[nti] >= 2000 and lowest_elevation > actual_neigh_elev then
+			lowest_elevation = actual_neigh_elev
+			lowest_nti = nti
+		end
+	end)
+
+	return lowest_nti
+end
+
+local function process_tributary(ti, wb, members)
+	true_river[ti] = true --* Set as true river so it can never be counted again as a waterbody member
+	table.insert(members, ti)
+
+	local lowest_nti = find_path_using_watershed(ti, wb)
+	if lowest_nti == -1 then return false, -1 end
+
+	--* If greater than -1, it implies we actually have a lower neighbor, so therefore we need to check if it is a true river
+
+	--* We also need to check to see if it is part of the same tributary as well.
+	if fork_count[lowest_nti] == fork_count[ti] then return true, lowest_nti end
+
+	--* We need to terminate expansion for this tributary, and construct a waterbody for the tributary using the list we created earlier
+	local new_tributary_wb = world:create_waterbody(waterbody.TYPES.river)
+	new_tributary_wb.tmp_float_1 = 0
+	new_tributary_wb.water_level = 0
+	for _, trib_ti in ipairs(members) do
+		world:add_tile_to_waterbody(trib_ti, new_tributary_wb)
+	end
+	while #members > 0 do table.remove(members) end
+
+	if true_river[lowest_nti] then return false, -1 end
+
+	return true, lowest_nti
+end
+
+local function split_river_up_into_tributaries()
+	--* I need to connect all my rivers, lakes, and oceans now...
+
+	for i = 0, #initial_candidates - 1 do
+		local ti = initial_candidates[sorted_candidates[i] + 1]
+
+		 --* If the tile has been already turned into a true river, don't bother continuing.
+		if true_river[ti] then goto continue2 end
+
+		local wb = watershed[ti]
+		if not wb then goto continue2 end
+
+		local tributary_members = {} --* This list will be turned into the members of the river tributary.
+		local found_path = true
+
+		while found_path do --* Continue to the ocean or until we hit a tile which has already been turned into a true river
+			found_path, ti = process_tributary(ti, wb, tributary_members)
+		end
+
+		::continue2::
+	end
+end
+
+local function reassign_proper_tile_waterbody_to_lakes()
+	world:for_each_waterbody(function(wb)
+		if wb.type == wb.TYPES.freshwater_lake then
+			wb:for_each_tile(function(ti)
+				world:reassign_tile_to_waterbody(ti, wb)
+				watershed[ti] = wb
+			end)
+		end
+
+		if wb.type == wb.TYPES.river then
+			local members_under_ice = 0
+			local total_members = wb:size()
+			wb:for_each_tile(function(ti)
+				local ice = world.ice[ti]
+				if ice == 0 then return end
+
+				if ice > 1000 then members_under_ice = members_under_ice + 2
+				elseif ice > 500 then members_under_ice = members_under_ice + 1.75
+				elseif ice > 200 then members_under_ice = members_under_ice + 1.5
+				elseif ice > 100 then members_under_ice = members_under_ice + 1.25
+				elseif ice > 50 then members_under_ice = members_under_ice + 1
+				elseif ice > 25 then members_under_ice = members_under_ice + 0.75
+				else members_under_ice = members_under_ice + 0.51
+				end
+			end)
+			if members_under_ice / total_members > 0.5 then
+				world:kill_waterbody(wb)
+			end
+		end
+	end)
+end
+
+local function connect_all_waterbodies()
+	world:for_each_waterbody(function(wb)
+		if wb.type == wb.TYPES.river then
+			--* Check first tile first to determine whether there are waterbody sources feeding into the current waterbody
+			local first_ti = wb.tiles[1]
+			local first_tile_elev = world:true_elevation_for_waterflow(first_ti)
+			world:for_each_neighbor(first_ti, function(nti)
+				local nwb = world:get_waterbody_by_tile(nti)
+				if not nwb or not nwb:is_valid() then return end
+
+				if world:true_elevation_for_waterflow(nti) > first_tile_elev then
+					wb:add_source(nwb)
+				end
+			end)
+			wb.lake_open = #wb.source > 0 and true or false --* Open means we're getting water from another waterbody and not just the ambient environment
+
+			--* Check last tile to determine the waterbodies that are being fed by the current waterbody. EVERY river should feed somewhere
+			local lowest_elevation = 100000
+			local lowest_wb = nil
+			local last_ti = wb.tiles[#wb.tiles]
+			world:for_each_neighbor(last_ti, function(nti)
+				local nwb = world:get_waterbody_by_tile(nti)
+				if not nwb or not nwb:is_valid() then return end
+
+				local elev_to_check = world:true_elevation_for_waterflow(nti)
+
+				if nwb:is_lake_or_ocean() then
+					elev_to_check = nwb.water_level
+				end
+
+				if elev_to_check < lowest_elevation then
+					lowest_elevation = elev_to_check
+					lowest_wb = nwb
+				end
+			end)
+			if lowest_wb == nil then error("River " .. wb.id .. " does not feed into a waterbody") end
+			wb.drain = lowest_wb
+
+			if wb.drain.type == wb.TYPES.freshwater_lake then
+				local lowest_shore_tile_wb = world:get_waterbody_by_tile(wb.drain.lowest_shore_tile)
+				wb.drain = lowest_shore_tile_wb and lowest_shore_tile_wb:is_valid() and lowest_shore_tile_wb or nil
+			end
+		end
+
+		if wb:is_lake_or_ocean() then
+			--* Receive water from sources
+			wb:for_each_tile_in_perimeter(function(ti)
+				--* Check for higher elevation than waterlevel... check for waterbody ID to make sure it is not zero and not different.
+				if world:true_elevation_for_waterflow(ti) <= wb.water_level then return end
+
+				--* If criteria is met, add as source
+				local nwb = world:get_waterbody_by_tile(ti)
+				if not nwb or not nwb:is_valid() then return end
+
+				wb:add_source(nwb)
+			end)
+		end
+
+		if wb.type == wb.TYPES.freshwater_lake then
+			local lowest_shore_tile_wb = world:get_waterbody_by_tile(wb.lowest_shore_tile)
+			if lowest_shore_tile_wb and lowest_shore_tile_wb:is_valid() then --* If the drain tile is a waterbody...
+				wb.drain = lowest_shore_tile_wb
+				lowest_shore_tile_wb:add_source(wb)
+			else
+				wb.lake_open = false
+				wb.drain = nil
+			end
+
+			if lowest_shore_tile_wb and lowest_shore_tile_wb.id == wb.id then
+				wb.lake_open = false
+				wb.drain = nil
+			end
+		end
+	end)
+end
+
+local function assigning_drainage_basin_value_to_rivers()
+	world:for_each_waterbody(function(wb)
+		if wb.type ~= wb.TYPES.river then return end
+		wb.basin = watershed[wb.tiles[1]]
+	end)
+end
+
+local function construct_wetlands()
+	world:fill_ffi_array(true_river, false)
+
+	--* Identify Wetlands
+	world:for_each_tile(function(ti)
+		local wb = world:get_waterbody_by_tile(ti)
+		if wb and wb:is_valid() then return end
+
+		if world.water_movement[ti] <= 2000 then return end
+
+		local points_to_check = 0
+		world:for_each_neighbor(ti, function(nti)
+			local nwb = world:get_waterbody_by_tile(nti)
+			if nwb and nwb:is_valid() or world.water_movement[nti] > 2000 then
+				points_to_check = points_to_check + 1
+			end
+		end)
+
+		if points_to_check >= 3 then
+			true_river[ti] = true
+		end
+	end)
+
+	--* Construct wetlands as actual waterbodies
+	world:for_each_tile(function(ti)
+		if not true_river[ti] or world.ice[ti] > 0 then return end
+
+		local wb = world:get_waterbody_by_tile(ti)
+		if wb and wb:is_valid() then return end
+
+		--* If elligible wetland but not already assigned to a waterbody
+		wb = world:create_waterbody_from_tile(ti, waterbody.TYPES.wetland)
+		wb.lake_open = true
+
+		local old_layer = {}
+		local new_layer = {}
+
+		table.insert(old_layer, ti)
+		local tiles_to_check = 1
+
+		while tiles_to_check > 0 do
+			for _, expansion_ti in ipairs(old_layer) do
+				world:for_each_neighbor(expansion_ti, function(nti)
+					if not true_river[nti] or world.ice[nti] > 0 then return end
+
+					local nwb = world:get_waterbody_by_tile(nti)
+					if nwb and nwb:is_valid() then return end
+
+					world:add_tile_to_waterbody(nti, wb)
+					table.insert(new_layer, nti)
+				end)
+			end
+
+			old_layer = {}
+			for _, new_ti in ipairs(new_layer) do
+				table.insert(old_layer, new_ti)
+			end
+			new_layer = {}
+			tiles_to_check = #old_layer
+		end
+	end)
+end
+
+function rivers.run(world_obj)
+	world = world_obj
+	true_lake = world.tmp_bool_1
+	true_river = world.tmp_bool_2
+	fork_count = world.tmp_int_1
+
+	world:fill_ffi_array(true_lake, false)
+	world:fill_ffi_array(true_river, false)
+	world:fill_ffi_array(fork_count, 0)
+
+	if enable_debug then
+		world:adjust_debug_channels(1)
+		world:reset_debug_all()
+	end
+
+	run_with_profiling(function() construct_start_locations() end, "construct_start_locations")
+	run_with_profiling(function()
+		sorted_candidates = require("libsote.heap-sort").heap_sort_indices(
+			function(i) return world:true_elevation_for_waterflow(initial_candidates[i + 1]) end,
+			nil,
+			#initial_candidates,
+			false
+		)
+	end, "sort_lowest_elevation_to_highest")
+	run_with_profiling(function() construct_drainage_basins() end, "construct_drainage_basins")
+	run_with_profiling(function() tag_and_prep_all_tributaries() end, "tag_and_prep_all_tributaries")
+	run_with_profiling(function() kill_old_basins() end, "kill_old_basins")
+	run_with_profiling(function() split_river_up_into_tributaries() end, "split_river_up_into_tributaries")
+	run_with_profiling(function() reassign_proper_tile_waterbody_to_lakes() end, "reassign_proper_tile_waterbody_to_lakes")
+	run_with_profiling(function() connect_all_waterbodies() end, "connect_all_waterbodies")
+	run_with_profiling(function() assigning_drainage_basin_value_to_rivers() end, "assigning_drainage_basin_value_to_rivers")
+	run_with_profiling(function() construct_wetlands() end, "construct_wetlands")
+end
+
+return rivers
+
+--* Make wetlands or create variables for rivers *--
+--*
+--* Major types of wetlands we need to represent:
+--*  Riverine
+--*  Coastal
+--*  Lake skirt
+
+
+
+--*  Nominal Plan? *--
+--* Okay, so we may want to generate basic soil data from bedrock weathering. We'll need a bit of global random variation layers for each sediment grain size.
+--* Sand, Silt, Clay... we may want to generate them first from the bedrock based on weathering rates, then transport them by water and *possibly* wind?
+--* Need weathering sources... need glacial and wind blown sources (particularly with silt), and need alluvial sources.
+--* 
+--* Theoretically we can generate the values of both of the first two in either order and aggregate them both... but the alluvial source should pull soil from upstream and
+--* deposit it down stream so we need the first two done.
+--*
+--* Generating soil texture material from bedrock is going to depend on local weathering conditions, but also the age of the rock. We'll also need a variability matrix similar
+--* to our randQuality matrix, but possibly for each type of texture? We can have a clay / sand slider variable and then a silt abundance variable. Then we use our last
+--* 3 logic variables for the purpose of transporting texture.
+--* We'll need to transport organics and mineral nutrients in a later job so that we're not juggling too many variables in the same job.
+--*
+--* Step 1: Generating variability matrix. Go back to C++ code for some inspiration.
diff --git a/sote/libsote/hydrology/open-issues.lua b/sote/libsote/hydrology/open-issues.lua
index 24ef2f5d..8360093d 100644
--- a/sote/libsote/hydrology/open-issues.lua
+++ b/sote/libsote/hydrology/open-issues.lua
@@ -24,6 +24,8 @@ function oi.waterflow_for_rank_7()
 end
 
 -- odd thing to do, wiping out water movement that was calculated in 'calculate-waterflow'
+-- 2024.09.28: might be ok, according to original authors;
+--             water movement might have been calculated and used in intermediary step(s) (like gen-parent-material), then wiped out and re-calculated
 function oi.set_water_movement_for_lakes(world, ti)
 	world.water_movement[ti] = 0
 end
diff --git a/sote/libsote/hydrology/waterbody.lua b/sote/libsote/hydrology/waterbody.lua
index db8eef45..57ba4763 100644
--- a/sote/libsote/hydrology/waterbody.lua
+++ b/sote/libsote/hydrology/waterbody.lua
@@ -17,10 +17,13 @@ function waterbody:new(id)
 	obj.id = id or 0
 	obj.tiles = {}
 	obj.type = waterbody.TYPES.invalid
+	obj.basin = nil
 	obj.water_level = 0
 	obj.perimeter = {}
 	obj.lowest_shore_tile = nil
 	obj.lake_open = false
+	obj.source = {}
+	obj.drain = nil
 	obj.tmp_float_1 = 0
 
 	setmetatable(obj, self)
@@ -32,12 +35,15 @@ end
 function waterbody:kill()
 	self.id = 0
 	self.tiles = {}
-	self.type = waterbody.TYPES.invalid
-	self.water_level = 0
+	-- self.type = waterbody.TYPES.invalid
+	self.basin = nil
+	-- self.water_level = 0
 	self.perimeter = {}
-	self.lowest_shore_tile = nil
-	self.lake_open = false
-	self.tmp_float_1 = 0
+	-- self.lowest_shore_tile = nil
+	-- self.lake_open = false
+	self.source = {}
+	self.drain = nil
+	-- self.tmp_float_1 = 0
 end
 
 ---@return number
@@ -50,6 +56,11 @@ function waterbody:add_tile(ti)
 	table.insert(self.tiles, ti)
 end
 
+---@param wb table
+function waterbody:add_source(wb)
+	table.insert(self.source, wb)
+end
+
 ---@param callback fun(tile_index:number)
 function waterbody:for_each_tile(callback)
 	for _, ti in ipairs(self.tiles) do
@@ -57,6 +68,13 @@ function waterbody:for_each_tile(callback)
 	end
 end
 
+---@param callback fun(tile_index:number)
+function waterbody:for_each_tile_in_perimeter(callback)
+	for ti, _ in pairs(self.perimeter) do
+		callback(ti)
+	end
+end
+
 ---@param ti number
 function waterbody:add_to_perimeter(ti)
 	self.perimeter[ti] = true
@@ -100,4 +118,13 @@ function waterbody:is_valid()
 	return self.id > 0
 end
 
+---@return boolean
+function waterbody:is_lake_or_ocean()
+	return self.type == waterbody.TYPES.freshwater_lake or self.type == waterbody.TYPES.saltwater_lake or self.type == waterbody.TYPES.ocean
+end
+
+function waterbody:is_salty()
+	return self.type == waterbody.TYPES.saltwater_lake or self.type == waterbody.TYPES.ocean
+end
+
 return waterbody
\ No newline at end of file
diff --git a/sote/libsote/world-generator.lua b/sote/libsote/world-generator.lua
index 70b422a5..7a906c98 100644
--- a/sote/libsote/world-generator.lua
+++ b/sote/libsote/world-generator.lua
@@ -40,10 +40,9 @@ end
 local fu = require "game.file-utils"
 
 local function override_climate_data()
-	local climate_generator = fu.csv_rows("d:\\temp\\sote\\12177\\sote_climate_data_by_elev.csv")
-	-- local logger = require("libsote.debug-loggers").get_climate_logger("d:/temp")
+	local climate_generator = fu.csv_rows("d:\\temp\\sote\\12177\\sote_climate_data.csv")
 
-	wg.world:for_each_tile_by_elevation_for_waterflow(function(ti, _)
+	wg.world:for_each_tile(function(ti, _)
 		local row = climate_generator()
 		if row == nil then
 			error("Not enough rows in climate data")
@@ -57,9 +56,6 @@ local function override_climate_data()
 		wg.world.jul_humidity[ti] = tonumber(row[8])
 		wg.world.jan_wind_speed[ti] = tonumber(row[9])
 		wg.world.jul_wind_speed[ti] = tonumber(row[10])
-
-		-- local log_str = row[1] .. "," .. row[2] .. " --- " .. wg.world.colatitude[ti] .. "," .. wg.world.minus_longitude[ti] .. " --- " .. wg.world:true_elevation_for_waterflow(ti) .. " <-> " .. tonumber(row[11])
-		-- logger:log(log_str)
 	end)
 end
 
@@ -127,6 +123,31 @@ local function gen_phase_02()
 	initial_waterflow()
 	glaciers()
 	run_with_profiling(function() require "libsote.hydrology.gen-dynamic-lakes".run(wg.world) end, "gen-dynamic-lakes")
+	run_with_profiling(function() require "libsote.hydrology.gen-rivers".run(wg.world) end, "gen-rivers")
+
+	local ocean_count = 0
+	local freshwater_lake_count = 0
+	local saltwater_lake_count = 0
+	local river_count = 0
+	local wetland_count = 0
+	wg.world:for_each_waterbody(function(wb)
+		if wb.type == wb.TYPES.ocean then
+			ocean_count = ocean_count + 1
+		elseif wb.type == wb.TYPES.freshwater_lake then
+			freshwater_lake_count = freshwater_lake_count + 1
+		elseif wb.type == wb.TYPES.saltwater_lake then
+			saltwater_lake_count = saltwater_lake_count + 1
+		elseif wb.type == wb.TYPES.river then
+			river_count = river_count + 1
+		elseif wb.type == wb.TYPES.wetland then
+			wetland_count = wetland_count + 1
+		end
+	end)
+	print("Ocean count: " .. ocean_count)
+	print("Freshwater lake count: " .. freshwater_lake_count)
+	print("Saltwater lake count: " .. saltwater_lake_count)
+	print("River count: " .. river_count)
+	print("Wetland count: " .. wetland_count)
 end
 
 local libsote_cpp = require "libsote.libsote"
diff --git a/sote/libsote/world-loader.lua b/sote/libsote/world-loader.lua
index a646c831..24045a9a 100644
--- a/sote/libsote/world-loader.lua
+++ b/sote/libsote/world-loader.lua
@@ -140,14 +140,21 @@ local function get_waterbody_color(wb)
 	end
 end
 
+local function get_empty_color(world, ti)
+	if world.is_land[ti] then
+		return 0, 0, 0
+	end
+
+	return 30, 30, 30
+end
+
 local rock_layers = require "libsote.rock-layers"
 
 function wl.load_maps_from(world)
 	local start = love.timer.getTime()
 
 	for _, tile in pairs(WORLD.tiles) do
-		local q, r, face = world:get_tile_coord(tile.tile_id)
-		local ti = world:get_tile_index(q, r, face)
+		local ti = world:get_tile_coord(tile.tile_id)
 
 		local is_land = world.is_land[ti]
 
@@ -255,6 +262,7 @@ function wl.dump_maps_from(world)
 	local image_jan_rainfall_data
 	local image_jan_waterflow_data
 	local image_waterbodies_data
+	local image_watersheds_data
 	local image_debug_data_1
 	local image_debug_data_2
 
@@ -273,8 +281,13 @@ function wl.dump_maps_from(world)
 	if debug_ms.waterbodies then
 		image_waterbodies_data = love.image.newImageData(width, height)
 	end
-	if debug_ms.debug then
+	if debug_ms.watersheds then
+		image_watersheds_data = love.image.newImageData(width, height)
+	end
+	if debug_ms.debug1 then
 		image_debug_data_1 = love.image.newImageData(width, height)
+	end
+	if debug_ms.debug2 then
 		image_debug_data_2 = love.image.newImageData(width, height)
 	end
 
@@ -359,14 +372,21 @@ function wl.dump_maps_from(world)
 				col_r, col_g, col_b = get_waterbody_color(wb)
 				image_waterbodies_data:setPixel(x, y, col_r / 255, col_g / 255, col_b / 255, 1)
 			end
+			-- watersheds ----------------------------------------------------
+			if debug_ms.watersheds then
+				if wb and wb:is_salty()  then
+					col_r, col_g, col_b = get_empty_color(world, ti)
+				else
+					col_r, col_g, col_b = get_waterbody_color(wb and wb.basin or nil)
+				end
+				image_watersheds_data:setPixel(x, y, col_r / 255, col_g / 255, col_b / 255, 1)
+			end
 
 			-- debug ---------------------------------------------------------
-			if debug_ms.debug then
+			if debug_ms.debug1 then
 				col_r, col_g, col_b, _ = world:get_debug_rgba_by_tile(ti, 1)
 				image_debug_data_1:setPixel(x, y, col_r / 255, col_g / 255, col_b / 255, 1)
 
-				col_r, col_g, col_b, _ = world:get_debug_rgba_by_tile(ti, 2)
-				image_debug_data_2:setPixel(x, y, col_r / 255, col_g / 255, col_b / 255, 1)
 				-- local r_blend, g_blend, b_blend, a_blend = world:get_debug_rgba(world.num_debug_channels, q, r, face)
 				-- for channel = world.num_debug_channels - 1, 1, -1 do
 				-- 	local cr, cg, cb, ca = world:get_debug_rgba(channel, q, r, face)
@@ -374,6 +394,10 @@ function wl.dump_maps_from(world)
 				-- end
 				-- image_debug_data:setPixel(x, y, r_blend / 255, g_blend / 255, b_blend / 255, a_blend)
 			end
+			if debug_ms.debug2 then
+				col_r, col_g, col_b, _ = world:get_debug_rgba_by_tile(ti, 2)
+				image_debug_data_2:setPixel(x, y, col_r / 255, col_g / 255, col_b / 255, 1)
+			end
 		end
 	end
 
@@ -397,11 +421,16 @@ function wl.dump_maps_from(world)
 		local waterbodies_file_data = image_waterbodies_data:encode('png')
 		love.filesystem.write(world.seed .. '_waterbodies.png', waterbodies_file_data)
 	end
-	if debug_ms.debug then
+	if debug_ms.watersheds then
+		local watersheds_file_data = image_watersheds_data:encode('png')
+		love.filesystem.write(world.seed .. '_watersheds.png', watersheds_file_data)
+	end
+	if debug_ms.debug1 then
 		local debug_file_data_1 = image_debug_data_1:encode('png')
-		local debug_file_data_2 = image_debug_data_2:encode('png')
-
 		love.filesystem.write(world.seed .. '_debug_1.png', debug_file_data_1)
+	end
+	if debug_ms.debug2 then
+		local debug_file_data_2 = image_debug_data_2:encode('png')
 		love.filesystem.write(world.seed .. '_debug_2.png', debug_file_data_2)
 	end
 end
diff --git a/sote/libsote/world.lua b/sote/libsote/world.lua
index 0ab02c3a..66422e63 100644
--- a/sote/libsote/world.lua
+++ b/sote/libsote/world.lua
@@ -1,12 +1,5 @@
 local world = {}
 
--- local transform = {
---     -0.86615, 0,       -0.49979, 0,
---     -0.17829, 0.93420,  0.30899, 0,
---      0.46690, 0.35674, -0.80916, 0,
---      0,       0,        0,       1
--- }
-
 local ffi = require("ffi")
 local ffi_mem_tally = 0
 
@@ -53,10 +46,11 @@ function world:new(world_size, seed)
 	obj.tile_count = obj.size * obj.size * 30 + 2
 	print("[world allocation] tile count: " .. obj.tile_count)
 	obj.coord = {}
-	obj.coord_by_tile_id = {} -- tile_id -> { q, r, face }, this is using cube world tile IDs
+	obj.square_to_hex = {} -- square tile id -> hex tile id, so map cube world tiles to hex world tiles
 	obj.coord_by_ti = {}      -- tile index -> { q, r, face }, this is using hex world 0-based tile indices
 	obj.climate_cells = {}
 	obj.waterbodies = {}
+	obj.killed_waterbodies = {}
 
 	obj.neighbors               = allocate_array("neighbors",               obj.tile_count * 6, "int32_t")
 	obj.waterbody_id_by_tile    = allocate_array("waterbody_id_by_tile",    obj.tile_count,     "uint32_t")
@@ -272,13 +266,16 @@ function world:get_hex_coords(ti)
 	return coord[1], coord[2], coord[3]
 end
 
-function world:cache_tile_coord(tile_id, q, r, face)
-	self.coord_by_tile_id[tile_id] = { q, r, face }
+function world:cache_square_ti_by_hex_ti(hex_ti)
+	table.insert(self.square_to_hex, hex_ti)
 end
 
-function world:get_tile_coord(tile_id)
-	local coord = self.coord_by_tile_id[tile_id]
-	return coord[1], coord[2], coord[3]
+function world:cache_square_ti_by_hex_coords(q, r, face)
+	table.insert(self.square_to_hex, self.coord[self:_key_from_coord(q, r, face)])
+end
+
+function world:get_tile_coord(square_ti)
+	return self.square_to_hex[square_ti]
 end
 
 function world:get_raw_colatitude(q, r, face)
@@ -468,24 +465,60 @@ end
 
 local waterbody = require "libsote.hydrology.waterbody"
 
----@param wb table waterbody
 ---@param ti number 0-based tile index
-function world:add_tile_to_waterbody(wb, ti)
+---@param wb table waterbody
+function world:add_tile_to_waterbody(ti, wb)
 	wb:add_tile(ti)
 	self.waterbody_id_by_tile[ti] = wb.id
 end
 
 ---@param ti number 0-based tile index
+---@param wb table waterbody
+function world:reassign_tile_to_waterbody(ti, wb)
+	self.waterbody_id_by_tile[ti] = wb.id
+end
+
+---@param wb_type waterbody_type
 ---@return table waterbody
-function world:create_new_waterbody_from_tile(ti)
-	-- no reuse of killed waterbodies for now
+function world:create_waterbody(wb_type)
+	if #self.killed_waterbodies > 0 then
+		local id = table.remove(self.killed_waterbodies, 1)
+		local wb = self.waterbodies[id]
+		wb.id = id
+		wb.type = wb_type
+		return wb
+	end
+
 	local id = #self.waterbodies + 1
 	local new_wb = waterbody:new(id)
-	self:add_tile_to_waterbody(new_wb, ti)
-	self.waterbodies[id] = new_wb
+	new_wb.type = wb_type
+	table.insert(self.waterbodies, new_wb)
 	return new_wb
 end
 
+---@param ti number 0-based tile index
+---@param wb_type waterbody_type
+---@return table waterbody
+function world:create_waterbody_from_tile(ti, wb_type)
+	local wb = self:create_waterbody(wb_type)
+	self:add_tile_to_waterbody(ti, wb)
+	return wb
+end
+
+---@param wb table waterbody
+function world:kill_waterbody(wb)
+	wb:for_each_tile(function(ti)
+		self.waterbody_id_by_tile[ti] = 0
+	end)
+
+	table.insert(self.killed_waterbodies, wb.id)
+	table.sort(self.killed_waterbodies)
+	wb:kill()
+end
+
+---@param q number
+---@param r number
+---@param face number
 function world:get_waterbody(q, r, face)
 	return self.waterbodies[self.waterbody_id_by_tile[self.coord[self:_key_from_coord(q, r, face)]]]
 end
@@ -513,7 +546,7 @@ end
 ---@param callback fun(waterbody:table)
 function world:for_each_waterbody(callback)
 	for _, wb in ipairs(self.waterbodies) do
-		if wb and wb:is_valid() then callback(wb) end
+		if wb:is_valid() then callback(wb) end
 	end
 end
 
@@ -521,14 +554,15 @@ end
 ---@param from_wb table waterbody to merge from
 function world:merge_waterbodies(to_wb, from_wb)
 	for _, ti in ipairs(from_wb.tiles) do
-		self:add_tile_to_waterbody(to_wb, ti)
+		self:add_tile_to_waterbody(ti, to_wb)
 	end
 	for ti, _ in pairs(from_wb.perimeter) do
 		to_wb:add_to_perimeter(ti)
 	end
 
+	table.insert(self.killed_waterbodies, from_wb.id)
+	table.sort(self.killed_waterbodies)
 	from_wb:kill()
-	-- no reuse of killed waterbodies for now, they are just left in the array
 end
 
 ---------------------------------------------------------------------------------------------------
@@ -686,7 +720,8 @@ end
 function world:sort_by_elevation_for_waterflow()
 	self.tiles_by_elevation_for_waterflow = require("libsote.heap-sort").heap_sort_indices(
 		function(i) return self:true_elevation_for_waterflow(i) end,
-		function(i) return self.colatitude[i] end,
+		nil,
+		-- function(i) return self.colatitude[i] end,
 		self.tile_count,
 		true, false)
 end