diff --git a/terrainlib/erosion.py b/terrainlib/erosion.py
index c8b2dc2..fbed7e6 100644
--- a/terrainlib/erosion.py
+++ b/terrainlib/erosion.py
@@ -11,12 +11,6 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
     v = K * flux^m
     """
 
-    dirs = dirs.copy()
-    dirs[0,:] = 0
-    dirs[-1,:] = 0
-    dirs[:,0] = 0
-    dirs[:,-1] = 0
-
     adv_time = 1 / (K*rivers**m) # For every pixel, calculate the time an "erosion wave" will need to cross it.
     dem = np.maximum(dem, sea_level)
     dem_new = np.zeros(dem.shape)
@@ -51,7 +45,7 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
             c = remaining / adv_time[y0,x0]
             dem_new[y,x] = c*dem[y1,x1] + (1-c)*dem[y0,x0] # If between 2 pixels, perform linear interpolation.
 
-    return np.minimum(dem, dem_new)
+    return dem_new
 
 def diffusion(dem, time, d=1):
     radius = d * time**.5
@@ -80,7 +74,7 @@ class EvolutionModel:
         self.flow_uptodate = True
 
     def advection(self, time):
-        dem = advection(self.lakes, self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
+        dem = advection(np.maximum(self.dem, self.lakes), self.dirs, self.rivers, time, K=self.K, m=self.m, sea_level=self.sea_level)
         self.dem = np.minimum(dem, self.dem)
         self.flow_uptodate = False
 
diff --git a/terrainlib/rivermapper.py b/terrainlib/rivermapper.py
index 687c14d..9e09d9d 100644
--- a/terrainlib/rivermapper.py
+++ b/terrainlib/rivermapper.py
@@ -153,7 +153,7 @@ def flow(dem):
 
     waterq = accumulate_flow(dirs2)
 
-    return dirs2, np.maximum(basins[basin_id], dem), waterq
+    return dirs2, basins[basin_id], waterq
 
 def accumulate_flow(dirs):
     ndonors = np.zeros(dirs.shape, dtype=int)
diff --git a/terrainlib/view.py b/terrainlib/view.py
index d3d9c78..9b03820 100644
--- a/terrainlib/view.py
+++ b/terrainlib/view.py
@@ -22,11 +22,11 @@ if has_matplotlib:
     def view_map(dem, lakes, scale=1, title=None):
         lakes_sea = np.maximum(lakes, 0)
         water = np.maximum(lakes_sea - dem, 0)
-        max_elev = lakes_sea.max()
+        max_elev = dem.max()
         max_depth = water.max()
 
         ls = mcl.LightSource(azdeg=315, altdeg=45)
-        rgb = ls.shade(lakes_sea, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', vmin=0, vmax=max_elev)
+        rgb = ls.shade(dem, cmap=cmap1, vert_exag=1/scale, blend_mode='soft', vmin=0, vmax=max_elev)
 
         (X, Y) = dem.shape
         extent = (0, Y*scale, 0, X*scale)
@@ -69,13 +69,13 @@ else:
     def plot(*args, **kwargs):
         pass
 
-def stats(dem, lake_dem, scale=1):
+def stats(dem, lakes, scale=1):
     surface = dem.size
 
-    continent = lake_dem >= 0
+    continent = np.maximum(dem, lakes) >= 0
     continent_surface = continent.sum()
 
-    lake = continent & (lake_dem>dem)
+    lake = continent & (lakes>dem)
     lake_surface = lake.sum()
 
     print('---   General    ---')