diff --git a/server/src/cmd.rs b/server/src/cmd.rs
index dcd21129ce..6171238892 100644
--- a/server/src/cmd.rs
+++ b/server/src/cmd.rs
@@ -990,6 +990,7 @@ fn handle_debug_column(server: &mut Server, entity: EcsEntity, args: String, act
         let foo = || {
             // let sim_chunk = sim.get(chunk_pos)?;
             let alt = sim.get_interpolated(wpos, |chunk| chunk.alt)?;
+            let basement = sim.get_interpolated(wpos, |chunk| chunk.basement)?;
             let water_alt = sim.get_interpolated(wpos, |chunk| chunk.water_alt)?;
             let chaos = sim.get_interpolated(wpos, |chunk| chunk.chaos)?;
             let temp = sim.get_interpolated(wpos, |chunk| chunk.temp)?;
@@ -1008,6 +1009,7 @@ fn handle_debug_column(server: &mut Server, entity: EcsEntity, args: String, act
                 r#"wpos: {:?}
 alt {:?} ({:?})
 water_alt {:?} ({:?})
+basement {:?}
 river {:?}
 downhill {:?}
 chaos {:?}
@@ -1022,6 +1024,7 @@ spawn_rate {:?} "#,
                 col.alt,
                 water_alt,
                 col.water_level,
+                basement,
                 river,
                 downhill,
                 chaos,
diff --git a/world/examples/water.rs b/world/examples/water.rs
index 92ba1e3715..335994ec96 100644
--- a/world/examples/water.rs
+++ b/world/examples/water.rs
@@ -26,6 +26,8 @@ fn main() {
     let light_direction = Vec3::new(-0.8, -1.0, 0.3).normalized();
     let light_res = 3;
 
+    let mut is_basement = false;
+
     while win.is_open() {
         let mut buf = vec![0; W * H];
         const QUADRANTS: usize = 4;
@@ -34,6 +36,8 @@ fn main() {
         let mut lakes = 0u32;
         let mut oceans = 0u32;
 
+        // let water_light = (light_direction.z + 1.0) / 2.0 * 0.8 + 0.2;
+
         for i in 0..W {
             for j in 0..H {
                 let pos = focus + Vec2::new(i as i32, j as i32) * scale;
@@ -41,11 +45,12 @@ fn main() {
                 let top_right = focus + Vec2::new(i as i32 + light_res, j as i32) * scale;
                 let bottom_left = focus + Vec2::new(i as i32, j as i32 + light_res) * scale; */
 
-                let (alt, water_alt, humidity, temperature, downhill, river_kind) = sampler
+                let (alt, basement, water_alt, humidity, temperature, downhill, river_kind) = sampler
                     .get(pos)
                     .map(|sample| {
                         (
                             sample.alt,
+                            sample.basement,
                             sample.water_alt,
                             sample.humidity,
                             sample.temp,
@@ -53,7 +58,7 @@ fn main() {
                             sample.river.river_kind,
                         )
                     })
-                    .unwrap_or((CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None));
+                    .unwrap_or((CONFIG.sea_level, CONFIG.sea_level, CONFIG.sea_level, 0.0, 0.0, None, None));
                 let humidity = humidity.min(1.0).max(0.0);
                 let temperature = temperature.min(1.0).max(-1.0) * 0.5 + 0.5;
                 let downhill_pos = (downhill
@@ -63,8 +68,9 @@ fn main() {
                     + pos;
                 let downhill_alt = sampler
                     .get(downhill_pos)
-                    .map(|s| s.alt)
+                    .map(|s| if is_basement { s.basement } else { s.alt })
                     .unwrap_or(CONFIG.sea_level);
+                let alt = if is_basement { basement } else { alt };
                 /* let alt_tl = sampler.get(top_left).map(|s| s.alt)
                     .unwrap_or(CONFIG.sea_level);
                 let alt_tr = sampler.get(top_right).map(|s| s.alt)
@@ -78,7 +84,7 @@ fn main() {
                         .map(|e| e as i32));
                 let cross_alt = sampler
                     .get(cross_pos)
-                    .map(|s| s.alt)
+                    .map(|s| if is_basement { s.basement } else { s.alt })
                     .unwrap_or(CONFIG.sea_level);
                 let forward_vec = Vec3::new(
                     (downhill_pos.x - pos.x) as f64,
@@ -98,7 +104,11 @@ fn main() {
                 let water_alt = ((alt.max(water_alt) - CONFIG.sea_level) as f64 / gain as f64)
                     .min(1.0)
                     .max(0.0);
-                let alt = ((alt - CONFIG.sea_level) as f64 / gain as f64)
+                let true_alt = (alt - CONFIG.sea_level) as f64 / gain as f64;
+                let water_depth = (water_alt - true_alt)
+                    .min(1.0)
+                    .max(0.0);
+                let alt = true_alt
                     .min(1.0)
                     .max(0.0);
                 let quad =
@@ -120,10 +130,15 @@ fn main() {
                 }
 
                 buf[j * W + i] = match river_kind {
-                    Some(RiverKind::Ocean) => u32::from_le_bytes([64, 32, 0, 255]),
+                    Some(RiverKind::Ocean) => u32::from_le_bytes([
+                        ((64.0 - water_depth * 64.0) * 1.0) as u8,
+                        ((32.0 - water_depth * 32.0) * 1.0) as u8,
+                        0,
+                        255,
+                    ]),
                     Some(RiverKind::Lake { .. }) => u32::from_le_bytes([
-                        64 + (water_alt * 191.0) as u8,
-                        32 + (water_alt * 95.0) as u8,
+                        (((64.0 + water_alt * 191.0) + (- water_depth * 64.0)) * 1.0) as u8,
+                        (((32.0 + water_alt * 95.0) + (- water_depth * 32.0)) * 1.0) as u8,
                         0,
                         255,
                     ]),
@@ -183,6 +198,9 @@ fn main() {
                 );
             }
         }
+        if win.is_key_down(minifb::Key::B) {
+            is_basement ^= true;
+        }
         if win.is_key_down(minifb::Key::W) {
             focus.y -= spd * scale;
         }
diff --git a/world/src/column/mod.rs b/world/src/column/mod.rs
index 5732e57519..dd231eac49 100644
--- a/world/src/column/mod.rs
+++ b/world/src/column/mod.rs
@@ -436,9 +436,12 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
         (temp < CONFIG.snow_temp ||
          downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/
         is_rocky {
-            sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
+            sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
+            // sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
+            // sim.get_interpolated(wpos, |chunk| chunk.alt)?
         } else {
             sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
+            // sim.get_interpolated(wpos, |chunk| chunk.alt)?
         };
 
         // Find the average distance to each neighboring body of water.
@@ -581,23 +584,6 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             .get((wposf_turb.div(128.0)).into_array()) as f32)
             .mul(24.0);
 
-        let riverless_alt_delta = (sim
-            .gen_ctx
-            .small_nz
-            .get((wposf_turb.div(/*200.0*//*50.0*//*24.0*//*56.0 / (chaos as f64).max(0.05)*/50.0)).into_array()) as f32)
-            .abs()
-            .mul(3.0)
-            /* .mul(chaos.max(0.05))
-            .mul(27.0) */
-            /* + (sim
-                .gen_ctx
-                .small_nz
-                .get((wposf_turb.div(400.0)).into_array()) as f32)
-                .abs()
-                .mul((1.0 - chaos).max(0.3))
-                .mul(1.0 - humidity)
-                .mul(32.0) */;
-
         let alt_for_river = alt
             + if overlap_count == 0.0 {
                 0.0
@@ -818,6 +804,28 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
             (false, alt_for_river, downhill_water_alt, 1.0)
         };
 
+        let riverless_alt_delta = (sim
+            .gen_ctx
+            .small_nz
+            .get((wposf_turb.div(/*200.0*/50.0/*24.0*//*56.0 / (chaos as f64).max(0.05)*//*50.0*/)).into_array()) as f32)
+            .min(1.0).max(-1.0)
+            // .mul(0.5).add(0.5)
+            .abs()
+            .mul(3.0)
+            /* .mul(chaos.max(0.05))
+            .mul(27.0)
+            + (sim
+                .gen_ctx
+                .small_nz
+                .get((wposf_turb.div(400.0)).into_array()) as f32)
+                .min(1.0).max(-1.0)
+                // .mul(0.5).add(0.5)
+                .abs()
+                .mul((1.0 - chaos).max(0.3))
+                .mul(1.0 - humidity)
+                .mul(32.0) */;
+
+
         let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor);
         let alt = alt_ + riverless_alt_delta;
 
diff --git a/world/src/sim/diffusion.rs b/world/src/sim/diffusion.rs
new file mode 100644
index 0000000000..08d2373556
--- /dev/null
+++ b/world/src/sim/diffusion.rs
@@ -0,0 +1,421 @@
+use rayon::prelude::*;
+
+/// From https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
+///
+/// See https://fastscape-lem.github.io/fastscapelib-fortran
+///
+/// nx = grid size x
+///
+/// ny = grid size y
+///
+/// xl = x-dimension of the model topography in meters (double precision)
+///
+/// yl = y-dimension of the model topography in meters (double precision)
+///
+/// dt = length of the time step in years (double precision)
+///
+/// ibc = boundary conditions for grid.  For now, we assume all four boundaries are fixed height,
+/// so this parameter is ignored.
+///
+/// h = heights of cells at each cell in the grid.
+///
+/// b = basement height at each cell in the grid (see https://en.wikipedia.org/wiki/Basement_(geology)).
+///
+/// kd = bedrock transport coefficient (or diffusivity) for hillslope processes in meter squared per year
+/// (double precision) at each cell in the grid.
+///
+/// kdsed = sediment transport coefficient (or diffusivity) for hillslope processes in meter squared per year;
+/// (double precision;) note that when kdsed < 0, its value is not used, i.e., kd for sediment and bedrock have
+/// the same value, regardless of sediment thickness
+/* subroutine Diffusion ()
+
+  ! subroutine to solve the diffusion equation by ADI
+
+  use FastScapeContext
+
+  implicit none
+*/
+pub fn diffusion(nx: usize, ny: usize, xl: f64, yl: f64, dt: f64, _ibc: (),
+                 h: &mut [/*f64*/f32], b: &mut [/*f64*/f32],
+                 kd: impl Fn(usize) -> f64, kdsed: f64,
+                 ) {
+    let aij = |i: usize, j: usize| j * nx + i;
+/*
+  double precision, dimension(:), allocatable :: f,diag,sup,inf,res
+  double precision, dimension(:,:), allocatable :: zint,kdint,zintp
+  integer i,j,ij
+  double precision factxp,factxm,factyp,factym,dx,dy
+*/
+    let mut f : Vec<f64>;
+    let mut diag : Vec<f64>;
+    let mut sup : Vec<f64>;
+    let mut inf : Vec<f64>;
+    let mut res : Vec<f64>;
+    let mut zint : Vec<f64>;
+    let mut kdint : Vec<f64>;
+    let mut zintp : Vec<f64>;
+    let mut i : usize;
+    let mut j : usize;
+    let mut ij : usize;
+    let mut factxp : f64;
+    let mut factxm : f64;
+    let mut factyp : f64;
+    let mut factym : f64;
+    let mut dx : f64;
+    let mut dy : f64;
+/*
+  character cbc*4
+
+  !print*,'Diffusion'
+
+  write (cbc,'(i4)') ibc
+
+  dx=xl/(nx-1)
+  dy=yl/(ny-1)
+*/
+    // 2048*32/2048/2048
+    // 1 / 64 m
+    dx = xl / /*(nx - 1)*/nx as f64;
+    dy = yl / /*(ny - 1)*/ny as f64;
+/*
+  ! creates 2D internal arrays to store topo and kd
+
+  allocate (zint(nx,ny),kdint(nx,ny),zintp(nx,ny))
+*/
+    zint = vec![Default::default(); nx * ny];
+    kdint = vec![Default::default(); nx * ny];
+    zintp = vec![Default::default(); nx * ny];
+/*
+  do j=1,ny
+    do i=1,nx
+      ij=(j-1)*nx+i
+      zint(i,j)=h(ij)
+      kdint(i,j)=kd(ij)
+      if (kdsed.gt.0.d0 .and. (h(ij)-b(ij)).gt.1.d-6) kdint(i,j)=kdsed
+    enddo
+  enddo
+
+  zintp = zint
+*/
+    for j in 0..ny {
+        for i in 0..nx {
+            // ij = vec2_as_uniform_idx(i, j);
+            ij = aij(i, j);
+            zint[ij] = h[ij] as f64;
+            kdint[ij] = kd(ij);
+            if kdsed > 0.0 && (h[ij] - b[ij]) > 1.0e-6 {
+                kdint[ij] = kdsed;
+            }
+        }
+    }
+
+    zintp = zint.clone();
+/*
+  ! first pass along the x-axis
+
+  allocate (f(nx),diag(nx),sup(nx),inf(nx),res(nx))
+  f=0.d0
+  diag=0.d0
+  sup=0.d0
+  inf=0.d0
+  res=0.d0
+  do j=2,ny-1
+*/
+    f = vec![0.0; nx];
+    diag = vec![0.0; nx];
+    sup = vec![0.0; nx];
+    inf = vec![0.0; nx];
+    res = vec![0.0; nx];
+    for j in 1..ny-1 {
+/*
+    do i=2,nx-1
+    factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+    factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+      factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+    factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+    diag(i)=1.d0+factxp+factxm
+    sup(i)=-factxp
+    inf(i)=-factxm
+    f(i)=zintp(i,j)+factyp*zintp(i,j+1)-(factyp+factym)*zintp(i,j)+factym*zintp(i,j-1)
+    enddo
+*/
+        for i in 1..nx-1 {
+            factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[i] = 1.0 + factxp + factxm;
+            sup[i] = -factxp;
+            inf[i] = -factxm;
+            f[i] =
+                zintp[aij(i, j)] + factyp * zintp[aij(i, j+1)] -
+                (factyp + factym) *
+                zintp[aij(i, j)] + factym * zintp[aij(i, j-1)];
+        }
+/*
+! left bc
+    if (cbc(4:4).eq.'1') then
+    diag(1)=1.
+    sup(1)=0.
+    f(1)=zintp(1,j)
+    else
+    factxp=(kdint(2,j)+kdint(1,j))/2.d0*(dt/2.)/dx**2
+    factyp=(kdint(1,j+1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
+    factym=(kdint(1,j-1)+kdint(1,j))/2.d0*(dt/2.)/dy**2
+    diag(1)=1.d0+factxp
+    sup(1)=-factxp
+    f(1)=zintp(1,j)+factyp*zintp(1,j+1)-(factyp+factym)*zintp(1,j)+factym*zintp(1,j-1)
+    endif
+*/
+        if true {
+            diag[0] = 1.0;
+            sup[0] = 0.0;
+            f[0] = zintp[aij(0, j)];
+        } else {
+            // FIXME: implement reflective boundary
+        }
+/*
+! right bc
+    if (cbc(2:2).eq.'1') then
+    diag(nx)=1.
+    inf(nx)=0.
+    f(nx)=zintp(nx,j)
+    else
+    factxm=(kdint(nx-1,j)+kdint(nx,j))/2.d0*(dt/2.)/dx**2
+    factyp=(kdint(nx,j+1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
+    factym=(kdint(nx,j-1)+kdint(nx,j))/2.d0*(dt/2.)/dy**2
+    diag(nx)=1.d0+factxm
+    inf(nx)=-factxm
+    f(nx)=zintp(nx,j)+factyp*zintp(nx,j+1)-(factyp+factym)*zintp(nx,j)+factym*zintp(nx,j-1)
+    endif
+*/
+        if true {
+            diag[nx - 1] = 1.0;
+            inf[nx - 1] = 0.0;
+            f[nx - 1] = zintp[aij(nx - 1, j)];
+        } else {
+            // FIXME: implement reflective boundary
+        }
+/*
+  call tridag (inf,diag,sup,f,res,nx)
+    do i=1,nx
+    zint(i,j)=res(i)
+    enddo
+*/
+        tridag(&inf, &diag, &sup, &f, &mut res, nx);
+        for i in 0..nx {
+            zint[aij(i, j)] = res[i];
+        }
+/*
+  enddo
+  deallocate (f,diag,sup,inf,res)
+*/
+    }
+
+/*
+  ! second pass along y-axis
+
+  allocate (f(ny),diag(ny),sup(ny),inf(ny),res(ny))
+  f=0.d0
+  diag=0.d0
+  sup=0.d0
+  inf=0.d0
+  res=0.d0
+  do i=2,nx-1
+*/
+    f = vec![0.0; ny];
+    diag = vec![0.0; ny];
+    sup = vec![0.0; ny];
+    inf = vec![0.0; ny];
+    res = vec![0.0; ny];
+    for i in 1..nx-1 {
+/*
+    do j=2,ny-1
+    factxp=(kdint(i+1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+    factxm=(kdint(i-1,j)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+    factyp=(kdint(i,j+1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+    factym=(kdint(i,j-1)+kdint(i,j))/2.d0*(dt/2.)/dy**2
+    diag(j)=1.d0+factyp+factym
+    sup(j)=-factyp
+    inf(j)=-factym
+    f(j)=zint(i,j)+factxp*zint(i+1,j)-(factxp+factxm)*zint(i,j)+factxm*zint(i-1,j)
+    enddo
+*/
+        for j in 1..ny-1 {
+            factxp = (kdint[aij(i+1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factxm = (kdint[aij(i-1, j)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dx * dx);
+            factyp = (kdint[aij(i, j+1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            factym = (kdint[aij(i, j-1)] + kdint[aij(i,j)]) / 2.0 * (dt / 2.0) / (dy * dy);
+            diag[j] = 1.0 + factyp + factym;
+            sup[j] = -factyp;
+            inf[j] = -factym;
+            f[j] =
+                zint[aij(i, j)] + factxp * zint[aij(i+1, j)] -
+                (factxp + factxm) *
+                zint[aij(i, j)] + factxm * zint[aij(i-1, j)];
+        }
+/*
+! bottom bc
+    if (cbc(1:1).eq.'1') then
+    diag(1)=1.
+    sup(1)=0.
+    f(1)=zint(i,1)
+    else
+    factxp=(kdint(i+1,1)+kdint(i,j))/2.d0*(dt/2.)/dx**2
+    factxm=(kdint(i-1,1)+kdint(i,1))/2.d0*(dt/2.)/dx**2
+    factyp=(kdint(i,2)+kdint(i,1))/2.d0*(dt/2.)/dy**2
+    diag(1)=1.d0+factyp
+    sup(1)=-factyp
+    f(1)=zint(i,1)+factxp*zint(i+1,1)-(factxp+factxm)*zint(i,1)+factxm*zint(i-1,1)
+    endif
+*/
+        if true {
+            diag[0] = 1.0;
+            sup[0] = 0.0;
+            f[0] = zint[aij(i, 0)];
+        } else {
+            // FIXME: implement reflective boundary
+        }
+/*
+! top bc
+    if (cbc(3:3).eq.'1') then
+    diag(ny)=1.
+    inf(ny)=0.
+    f(ny)=zint(i,ny)
+    else
+    factxp=(kdint(i+1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
+    factxm=(kdint(i-1,ny)+kdint(i,ny))/2.d0*(dt/2.)/dx**2
+    factym=(kdint(i,ny-1)+kdint(i,ny))/2.d0*(dt/2.)/dy**2
+    diag(ny)=1.d0+factym
+    inf(ny)=-factym
+    f(ny)=zint(i,ny)+factxp*zint(i+1,ny)-(factxp+factxm)*zint(i,ny)+factxm*zint(i-1,ny)
+    endif
+*/
+        if true {
+            diag[ny - 1] = 1.0;
+            inf[ny - 1] = 0.0;
+            f[ny - 1] = zint[aij(i, ny - 1)];
+        } else {
+            // FIXME: implement reflective boundary
+        }
+/*
+  call tridag (inf,diag,sup,f,res,ny)
+    do j=1,ny
+    zintp(i,j)=res(j)
+    enddo
+*/
+        tridag(&inf, &diag, &sup, &f, &mut res, ny);
+        for j in 0..ny {
+            zintp[aij(i, j)] = res[j];
+        }
+/*
+  enddo
+  deallocate (f,diag,sup,inf,res)
+*/
+    }
+/*
+  ! stores result in 1D array
+
+  do j=1,ny
+    do i=1,nx
+      ij=(j-1)*nx+i
+      etot(ij)=etot(ij)+h(ij)-zintp(i,j)
+      erate(ij)=erate(ij)+(h(ij)-zintp(i,j))/dt
+      h(ij)=zintp(i,j)
+    enddo
+  enddo
+
+  b=min(h,b)
+*/
+    for j in 0..ny {
+        for i in 0..nx {
+            ij = aij(i, j);
+            // FIXME: Track total erosion and erosion rate.
+            h[aij(i, j)] = zintp[aij(i, j)] as f32;
+        }
+    }
+
+    b.par_iter_mut().zip(h).for_each(|(mut b, h)| {
+        *b = h.min(*b);
+    });
+/*
+  deallocate (zint,kdint,zintp)
+
+  return
+
+end subroutine Diffusion
+*/
+}
+
+/*
+!----------
+
+! subroutine to solve a tri-diagonal system of equations (from Numerical Recipes)
+
+      SUBROUTINE tridag(a,b,c,r,u,n)
+
+      implicit none
+
+      INTEGER n
+      double precision a(n),b(n),c(n),r(n),u(n)
+*/
+pub fn tridag(a: &[f64], b: &[f64], c: &[f64], r: &[f64], u: &mut [f64], n: usize) {
+/*
+      INTEGER j
+      double precision bet
+      double precision,dimension(:),allocatable::gam
+
+      allocate (gam(n))
+
+      if(b(1).eq.0.d0) stop 'in tridag'
+*/
+    let mut j : usize;
+    let mut bet : f64;
+    let mut precision : f64;
+    let mut gam: Vec<f64>;
+
+    gam = vec![Default::default(); n];
+
+    assert!(b[0] != 0.0);
+/*
+
+! first pass
+
+bet=b(1)
+u(1)=r(1)/bet
+do 11 j=2,n
+  gam(j)=c(j-1)/bet
+  bet=b(j)-a(j)*gam(j)
+  if(bet.eq.0.) then
+    print*,'tridag failed'
+    stop
+  endif
+  u(j)=(r(j)-a(j)*u(j-1))/bet
+  11    continue
+*/
+    bet = b[0];
+    u[0] = r[0] / bet;
+    for j in 1..n {
+        gam[j] = c[j - 1] / bet;
+        bet = b[j] - a[j] * gam[j];
+        assert!(bet != 0.0);
+        u[j] = (r[j] - a[j] * u[j - 1]) / bet;
+    }
+/*
+  ! second pass
+
+  do 12 j=n-1,1,-1
+    u(j)=u(j)-gam(j+1)*u(j+1)
+    12    continue
+*/
+    for j in (0..n - 1).rev() {
+        u[j] = u[j] - gam[j + 1] * u[j + 1];
+    }
+/*
+    deallocate (gam)
+
+    return
+
+    END
+*/
+}
diff --git a/world/src/sim/erosion.rs b/world/src/sim/erosion.rs
index 672ebadbe8..f7ffc57119 100644
--- a/world/src/sim/erosion.rs
+++ b/world/src/sim/erosion.rs
@@ -1,4 +1,5 @@
-use super::{downhill, neighbors, uniform_idx_as_vec2, uphill, WORLD_SIZE};
+use super::{diffusion, downhill, neighbors, uniform_idx_as_vec2, uphill, WORLD_SIZE};
+use bitvec::prelude::{bitbox, bitvec, BitBox};
 use crate::{config::CONFIG, util::RandomField};
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
 use noise::{NoiseFn, Point3};
@@ -157,6 +158,7 @@ pub fn get_rivers(
     let neighbor_coef = TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
     // NOTE: This technically makes us discontinuous, so we should be cautious about using this.
     let derivative_divisor = 1.0;
+    let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
     for &chunk_idx in newh.into_iter().rev() {
         let chunk_idx = chunk_idx as usize;
         let downhill_idx = downhill[chunk_idx];
@@ -351,7 +353,7 @@ pub fn get_rivers(
         // network).
         let downhill_water_alt = water_alt[downhill_idx];
         let neighbor_distance = neighbor_dim.magnitude();
-        let dz = (downhill_water_alt - chunk_water_alt) * CONFIG.mountain_scale;
+        let dz = (downhill_water_alt - chunk_water_alt) / height_scale as f32;// * CONFIG.mountain_scale;
         let slope = dz.abs() as f64 / neighbor_distance;
         if slope == 0.0 {
             // This is not a river--how did this even happen?
@@ -473,23 +475,28 @@ pub fn get_rivers(
 /// Precompute the maximum slope at all points.
 ///
 /// TODO: See if allocating in advance is worthwhile.
-fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f32]> {
-    const MIN_MAX_ANGLE: f64 = 6.0 / 360.0 * 2.0 * f64::consts::PI;
-    const MAX_MAX_ANGLE: f64 = 54.0 / 360.0 * 2.0 * f64::consts::PI;
+fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync)) -> Box<[f64]> {
+    const MIN_MAX_ANGLE: f64 = 6.0/*30.0*//*6.0*//*15.0*/ / 360.0 * 2.0 * f64::consts::PI;
+    const MAX_MAX_ANGLE: f64 = 54.0/*50.0*//*54.0*//*45.0*/ / 360.0 * 2.0 * f64::consts::PI;
     const MAX_ANGLE_RANGE: f64 = MAX_MAX_ANGLE - MIN_MAX_ANGLE;
+    let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
     h.par_iter()
         .enumerate()
         .map(|(posi, &z)| {
-            let wposf = (uniform_idx_as_vec2(posi)).map(|e| e as f64);
-            let wposz = z as f64 * CONFIG.mountain_scale as f64;
+            let wposf = uniform_idx_as_vec2(posi).map(|e| e as f64) * TerrainChunkSize::RECT_SIZE.map(|e| e as f64);
+            let wposz = z as f64 / height_scale;// * CONFIG.mountain_scale as f64;
             // Normalized to be between 6 and and 54 degrees.
-            let div_factor = 32.0;
+            let div_factor = /*32.0*//*16.0*//*64.0*/256.0/*1.0*//*4.0*/ * /*1.0*/4.0/* TerrainChunkSize::RECT_SIZE.x as f64 / 8.0 */;
             let rock_strength = rock_strength_nz
                 .get([
                     wposf.x, /* / div_factor*/
                     wposf.y, /* / div_factor*/
-                    wposz,
-                ])
+                    wposz * div_factor,
+                ]);
+            /* if rock_strength < -1.0 || rock_strength > 1.0 {
+                println!("Nooooo: {:?}", rock_strength);
+            } */
+            let rock_strength = rock_strength
                 .max(-1.0)
                 .min(1.0)
                 * 0.5
@@ -508,20 +515,22 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
             // We do log-odds against center, so that our log odds are 0 when x = 0.25, lower when x is
             // lower, and higher when x is higher.
             //
+            // (NOTE: below sea level, we invert it).
+            //
             // TODO: Make all this stuff configurable... but honestly, it's so complicated that I'm not
             // sure anyone would be able to usefully tweak them on a per-map basis?  Plus it's just a
             // hacky heuristic anyway.
             let center = 0.25;
-            let delta = 0.05;
-            let dmin = center - delta;
-            let dmax = center + delta;
+            let dmin = center - 0.15;//0.05;
+            let dmax = center + 0.05;//0.05;
             let log_odds = |x: f64| logit(x) - logit(center);
             let rock_strength = logistic_cdf(
                 1.0 * logit(rock_strength.min(1.0f64 - 1e-7).max(1e-7))
-                    + 1.0 * log_odds((z as f64).min(dmax).max(dmin)),
+                    + 1.0 * log_odds((wposz / CONFIG.mountain_scale as f64).abs().min(dmax).max(dmin)),
             );
+            // let rock_strength = 0.5;
             let max_slope = (rock_strength * MAX_ANGLE_RANGE + MIN_MAX_ANGLE).tan();
-            max_slope as f32
+            max_slope
         })
         .collect::<Vec<_>>()
         .into_boxed_slice()
@@ -569,27 +578,73 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
 ///    always be the *new* h[j](t + δt), while h[i] will still not have been computed yet, so we
 ///    only need to visit each node once.
 ///
+/// Afterwards, we also apply a hillslope diffusion process using an ADI (alternating direction
+/// implicit) method:
+///
+/// https://github.com/fastscape-lem/fastscapelib-fortran/blob/master/src/Diffusion.f90
+///
 /// [1] Guillaume Cordonnier, Jean Braun, Marie-Paule Cani, Bedrich Benes, Eric Galin, et al..
 ///     Large Scale Terrain Generation from Tectonic Uplift and Fluvial Erosion.
 ///     Computer Graphics Forum, Wiley, 2016, Proc. EUROGRAPHICS 2016, 35 (2), pp.165-175.
 ///     ⟨10.1111/cgf.12820⟩. ⟨hal-01262376⟩
 ///
+///
 fn erode(
     h: &mut [f32],
+    b: &mut [f32],
+    is_done: &mut BitBox,
+    done_val: bool,
     erosion_base: f32,
     max_uplift: f32,
+    kdsed: f64,
     _seed: &RandomField,
     rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
     uplift: impl Fn(usize) -> f32,
+    kd: impl Fn(usize) -> f64,
     is_ocean: impl Fn(usize) -> bool + Sync,
 ) {
     log::debug!("Done draining...");
-    let mmaxh = 1.0;
-    // Landslide constant: ideally scaled to 10 m^-2 / y^-1
-    let l = 200.0 * max_uplift as f64;
-    // Erosion constant.
-    let k = erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
-    let ((dh, indirection, newh, area), max_slope) = rayon::join(
+    let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
+    let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
+    // Since maximum uplift rate is expected to be 5.010e-4 m * y^-1, and
+    // 1.0 height units is 1.0 / height_scale m, whatever the
+    // max uplift rate is (in units / y), we can find dt by multiplying by
+    // 1.0 / height_scale m / unit and then dividing by 5.010e-4 m / y
+    // (to get dt in y / unit).  More formally:
+    //
+    // max_uplift h_unit / dt y = 5.010e-4 m / y
+    //
+    // 1 h_unit = 1.0 / height_scale m
+    //
+    //   max_uplift h_unit / dt * 1.0 / height_scale m / h_unit =
+    //   max_uplift / height_scale m / dt =
+    //   5.010e-4 m / y
+    //
+    //   max_uplift / height_scale m / dt / (5.010e-4 m / y) = 1
+    //   (max_uplift / height_scale / 5.010e-4) y = dt
+    let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
+    println!("dt={:?}", dt);
+    // Landslide constant: ideally scaled to 10e-2 m / y^-1
+    let l = /*200.0 * max_uplift as f64;*/10.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
+    // Net precipitation rate (m / year)
+    let p = 1.0 * height_scale;
+    // Stream power erosion constant (bedrock), in m^(1-2m) / year  (times dt).
+    let k_fb = // erosion_base as f64 + 2.244 / mmaxh as f64 * max_uplift as f64;
+        // 2.244*(5.010e-4)/512*5- (1.097e-5)
+        // 2.244*(5.010e-4)/2048*5- (1.097e-5)
+        // 2.244*(5.010e-4)/512- (8e-6)
+        // 2.244*(5.010e-4)/512- (2e-6)
+        // 2e-6 * dt;
+        // 8e-6 * dt
+        // 2e-5 * dt;
+        // 2.244/2048*5*32/(250000/4)*10^6
+        erosion_base as f64 + 2.244 / mmaxh as f64 * 5.0 * max_uplift as f64;
+        // 1e-5 * dt;
+        // (2.244*(5.010e-4)/512)/(2.244*(5.010e-4)/2500) = 4.88...
+        // 2.444 * 5
+    // Stream power erosion constant (sediment), in m^(1-2m) / year (times dt).
+    let k_fs = k_fb * 1.0/* 2.0*//*4.0*/;
+    let ((dh, indirection, newh, area), mut max_slopes) = rayon::join(
         || {
             let mut dh = downhill(h, |posi| is_ocean(posi));
             log::debug!("Computed downhill...");
@@ -622,8 +677,10 @@ fn erode(
             if posj == -1 {
                 panic!("Disconnected lake!");
             }
+            // max_slopes[posi] = kd(posi);
         // Egress with no outgoing flows.
         } else {
+            // *is_done.at(posi) = done_val;
             nland += 1;
             let posj = posj as usize;
             let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
@@ -633,23 +690,53 @@ fn erode(
             let neighbor_distance = (neighbor_coef * dxy).magnitude();
             // Since the area is in meters^(2m) and neighbor_distance is in m, so long as m = 0.5,
             // we have meters^(1) / meters^(1), so they should cancel out.  Otherwise, we would
-            // want to divide neighbor_distance by CONFIG.mountain_scale and area[posi] by
-            // CONFIG.mountain_scale^2, to make sure we were working in the correct units for dz
-            // (which has 1.0 height unit = CONFIG.mountain_scale meters).
-            let flux = k * (chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
+            // want to multiply neighbor_distance by height_scale and area[posi] by
+            // height_scale^2, to make sure we were working in the correct units for dz
+            // (which has height_scale height unit = 1.0 meters).
+            let old_h_i = h[posi] as f64;
+            let uplift_i = uplift(posi) as f64;
+            let k = if (old_h_i - b[posi] as f64) > 1.0e-6 {
+                // Sediment
+                k_fs
+            } else {
+                // Bedrock
+                k_fb
+            };
+            // let k = k * uplift_i / max_uplift as f64;
+            let flux = k * (p * chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
             // h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
             // NOTE: posj has already been computed since it's downhill from us.
             let h_j = h[posj] as f64;
-            let old_h_i = h[posi] as f64;
-            let mut new_h_i = (old_h_i + (uplift(posi) as f64 + flux * h_j)) / (1.0 + flux);
+            // hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
+            let mut new_h_i = (old_h_i + (uplift_i + flux * h_j)) / (1.0 + flux);
+            // Prevent erosion from dropping us below our receiver, unless we were already below it.
+            // new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
             // Find out if this is a lake bottom.
             let indirection_idx = indirection[posi];
             let is_lake_bottom = indirection_idx < 0;
-            // Test the slope.
-            let max_slope = max_slope[posi] as f64;
-            let dz = (new_h_i - h_j).max(0.0) * CONFIG.mountain_scale as f64;
-            let mag_slope = dz/*.abs()*/ / neighbor_distance;
             let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0;
+            // Test the slope.
+            let max_slope = max_slopes[posi] as f64;
+            // Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute
+            // uplift to other neighbors.
+            let (posk, h_k) = /* neighbors(posi)
+                .filter(|&posk| *is_done.at(posk) == done_val)
+                // .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk))
+                .map(|posk| (posk, h[posk] as f64))
+                // .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0)
+                .min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap())
+                .unwrap_or((posj, h_j)); */
+                (posj, h_j);
+                // .max(h_j);
+            let (posk, h_k) = if h_k < h_j {
+                (posk, h_k)
+            } else {
+                (posj, h_j)
+            };
+            let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
+            let neighbor_distance = (neighbor_coef * dxy).magnitude();
+            let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
+            let mag_slope = dz/*.abs()*/ / neighbor_distance;
             // If you're on the lake bottom and not right next to your neighbor, don't compute a
             // slope.
             if
@@ -661,26 +748,57 @@ fn erode(
                     // println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i);
                     // Thermal erosion says this can't happen, so we reduce dh_i to make the slope
                     // exactly max_slope.
-                    // max_slope = (old_h_i + dh - h_j) * CONFIG.mountain_scale / NEIGHBOR_DISTANCE
-                    // dh = max_slope * NEIGHBOR_DISTANCE / CONFIG.mountain_scale + h_j - old_h_i.
-                    let dh = max_slope * neighbor_distance / CONFIG.mountain_scale as f64;
-                    new_h_i = (h_j + dh).max(new_h_i - l * (mag_slope - max_slope));
-                    let dz = (new_h_i - h_j).max(0.0) * CONFIG.mountain_scale as f64;
+                    // max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
+                    // dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
+                    let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
+                    new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
+                    let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
                     let slope = dz/*.abs()*/ / neighbor_distance;
                     sums += slope;
+                    /* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
+                    sums += mag_slope; */
                 // let slope = dz.signum() * max_slope;
-                // new_h_i = slope * neighbor_distance / CONFIG.mountain_scale as f64 + h_j;
+                // new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
                 // sums += max_slope;
                 } else {
+                    // max_slopes[posi] = 0.0;
                     sums += mag_slope;
                     // Just use the computed rate.
                 }
                 h[posi] = new_h_i as f32;
+                // Make sure to update the basement as well!
+                b[posi] = (old_h_i + uplift_i).min(new_h_i) as f32;
                 sumh += new_h_i;
             }
         }
+        *is_done.at(posi) = done_val;
         maxh = h[posi].max(maxh);
     }
+    log::debug!(
+        "Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
+        maxh,
+        if nland == 0 {
+            f64::INFINITY
+        } else {
+            sumh / nland as f64
+        },
+        if nland == 0 {
+            f64::INFINITY
+        } else {
+            sums / nland as f64
+        },
+    );
+
+    // Apply hillslope diffusion.
+    diffusion(WORLD_SIZE.x, WORLD_SIZE.y,
+              WORLD_SIZE.x as f64 * TerrainChunkSize::RECT_SIZE.x as f64 * height_scale/* / CONFIG.mountain_scale as f64*/,
+              WORLD_SIZE.y as f64 * TerrainChunkSize::RECT_SIZE.y as f64 * height_scale/* / CONFIG.mountain_scale as f64*/,
+              dt,
+              (),
+              h, b,
+              /*|posi| max_slopes[posi]*/kd,
+              kdsed,
+    );
     log::debug!(
         "Done eroding (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
         maxh,
@@ -708,7 +826,7 @@ pub fn fill_sinks(
 ) -> Box<[f32]> {
     // NOTE: We are using the "exact" version of depression-filling, which is slower but doesn't
     // change altitudes.
-    let epsilon = /*1.0 / (1 << 7) as f32 / CONFIG.mountain_scale*/0.0;
+    let epsilon = /*1.0 / (1 << 7) as f32 * height_scale/* / CONFIG.mountain_scale */*/0.0;
     let infinity = f32::INFINITY;
     let range = 0..WORLD_SIZE.x * WORLD_SIZE.y;
     let oldh = range
@@ -1122,14 +1240,21 @@ pub fn do_erosion(
     seed: &RandomField,
     rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync),
     oldh: impl Fn(usize) -> f32 + Sync,
+    oldb: impl Fn(usize) -> f32 + Sync,
     is_ocean: impl Fn(usize) -> bool + Sync,
     uplift: impl Fn(usize) -> f32 + Sync,
-) -> Box<[f32]> {
+) -> (Box<[f32]>, Box<[f32]>) {
     let oldh_ = (0..WORLD_SIZE.x * WORLD_SIZE.y)
         .into_par_iter()
         .map(|posi| oldh(posi))
         .collect::<Vec<_>>()
         .into_boxed_slice();
+    // Topographic basement (The height of bedrock, not including sediment).
+    let mut b = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+        .into_par_iter()
+        .map(|posi| oldb(posi))
+        .collect::<Vec<_>>()
+        .into_boxed_slice();
     // TODO: Don't do this, maybe?
     // (To elaborate, maybe we should have varying uplift or compute it some other way).
     let uplift = (0..oldh_.len())
@@ -1150,18 +1275,37 @@ pub fn do_erosion(
         .max_by(|a, b| a.partial_cmp(&b).unwrap())
         .unwrap();
     log::debug!("Max uplift: {:?}", max_uplift);
+    // Height of terrain, including sediment.
     let mut h = oldh_;
+    // 0.01 / 2e-5 = 500
+    // Bedrock transport coefficients (diffusivity) in m^2 / year.  For now, we set them all to be equal
+    // on land, but in theory we probably want to at least differentiate between soil, bedrock, and
+    // sediment.
+    let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
+    let kdsed = /*1.5e-2*//*1e-4*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
+    let kd_bedrock = /*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
+    let kd = |posi: usize| if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
+    // Hillslope diffusion coefficient for sediment.
+    let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
     for i in 0..n {
         log::debug!("Erosion iteration #{:?}", i);
         erode(
             &mut h,
+            &mut b,
+            &mut is_done,
+            // The value to use to indicate that erosion is complete on a chunk.  Should toggle
+            // once per iteration, to avoid having to reset the bits, and start at true, since
+            // we initialize to 0 (false).
+            i & 1 == 0,
             erosion_base,
             max_uplift,
+            kdsed,
             seed,
             rock_strength_nz,
             |posi| uplift[posi],
+            kd,
             |posi| is_ocean(posi),
         );
     }
-    h
+    (h, b)
 }
diff --git a/world/src/sim/mod.rs b/world/src/sim/mod.rs
index 67fb439c47..0076945d1b 100644
--- a/world/src/sim/mod.rs
+++ b/world/src/sim/mod.rs
@@ -1,16 +1,19 @@
+mod diffusion;
 mod erosion;
 mod location;
 mod settlement;
 mod util;
 
 // Reexports
+pub use self::diffusion::diffusion;
 pub use self::erosion::{
     do_erosion, fill_sinks, get_drainage, get_lakes, get_rivers, RiverData, RiverKind,
 };
 pub use self::location::Location;
 pub use self::settlement::Settlement;
 pub use self::util::{
-    cdf_irwin_hall, downhill, get_oceans, local_cells, map_edge_factor, neighbors,
+    cdf_irwin_hall, downhill, get_oceans, HybridMulti as HybridMulti_, local_cells, map_edge_factor, neighbors,
+    ScaleBias,
     uniform_idx_as_vec2, uniform_noise, uphill, vec2_as_uniform_idx, InverseCdf,
 };
 
@@ -60,6 +63,7 @@ struct GenCdf {
     temp_base: InverseCdf,
     chaos: InverseCdf,
     alt: Box<[f32]>,
+    basement: Box<[f32]>,
     water_alt: Box<[f32]>,
     dh: Box<[isize]>,
     /// NOTE: Until we hit 4096 × 4096, this should suffice since integers with an absolute value
@@ -74,7 +78,7 @@ pub(crate) struct GenCtx {
     pub turb_x_nz: SuperSimplex,
     pub turb_y_nz: SuperSimplex,
     pub chaos_nz: RidgedMulti,
-    pub alt_nz: HybridMulti,
+    pub alt_nz: HybridMulti_,
     pub hill_nz: SuperSimplex,
     pub temp_nz: Fbm,
     // Humidity noise
@@ -111,28 +115,34 @@ pub struct WorldSim {
 impl WorldSim {
     pub fn generate(seed: u32) -> Self {
         let mut rng = ChaChaRng::from_seed(seed_expan::rng_state(seed));
-        let continent_scale = 5_000.0/*32768.0*/;
+        let continent_scale = 5_000.0f64/*32768.0*/.div(32.0).mul(TerrainChunkSize::RECT_SIZE.x as f64);
 
         let gen_ctx = GenCtx {
             turb_x_nz: SuperSimplex::new().set_seed(rng.gen()),
             turb_y_nz: SuperSimplex::new().set_seed(rng.gen()),
             chaos_nz: RidgedMulti::new()
-                .set_octaves(/*7*//*3*/ 7)
+                .set_octaves(/*7*//*3*/ /*7*//*3*/7)
                 .set_frequency(
                     /*RidgedMulti::DEFAULT_FREQUENCY **/ 3_000.0 * 8.0 / continent_scale,
                 )
+                // .set_persistence(RidgedMulti::DEFAULT_LACUNARITY.powf(-(1.0 - 0.5)))
                 .set_seed(rng.gen()),
             hill_nz: SuperSimplex::new().set_seed(rng.gen()),
-            alt_nz: HybridMulti::new()
-                .set_octaves(/*3*//*2*/ 8)
+            alt_nz: HybridMulti_::new()
+                .set_octaves(/*3*//*2*/ /*8*//*3*/8)
                 // 1/2048*32*1024 = 16
                 .set_frequency(
                     /*HybridMulti::DEFAULT_FREQUENCY*/
+                    // (2^8*(10000/5000/10000))*32 = per-chunk
                     (10_000.0/* * 2.0*/ / continent_scale) as f64,
                 )
                 // .set_frequency(1.0 / ((1 << 0) as f64))
                 // .set_lacunarity(1.0)
-                .set_persistence(/*0.5*//*0.5*/ 0.5)
+                // persistence = lacunarity^(-(1.0 - fractal increment))
+                .set_persistence(HybridMulti_::DEFAULT_LACUNARITY.powf(-(1.0 - /*0.75*/0.75)))
+                // .set_persistence(/*0.5*//*0.5*/0.5 + 1.0 / ((1 << 6) as f64))
+                // .set_offset(/*0.7*//*0.5*//*0.75*/0.7)
+                // .set_offset(/*0.7*//*0.5*//*0.75*/0.0)
                 .set_seed(rng.gen()),
             //temp_nz: SuperSimplex::new().set_seed(rng.gen()),
             temp_nz: Fbm::new()
@@ -178,17 +188,81 @@ impl WorldSim {
         };
 
         let river_seed = RandomField::new(rng.gen());
-        let rock_strength_nz = Fbm::new()
-            .set_octaves(8)
-            .set_persistence(/*0.9*/ 2.0)
-            .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (64.0 * 32.0))
+        let rock_lacunarity = 0.5/*HybridMulti::DEFAULT_LACUNARITY*/;
+        let rock_strength_nz = /*Fbm*/HybridMulti_/*BasicMulti*//*Fbm*/::new()
+            .set_octaves(/*6*//*5*//*4*//*5*/4)
+            // persistence = lacunarity^(-(1.0 - fractal increment))
+            .set_persistence(/*0.9*/ /*2.0*//*1.5*//*HybridMulti::DEFAULT_LACUNARITY*/rock_lacunarity.powf(-(1.0 - 0.9)))
+            // 256*32/2^4
+            // (0.5^(-(1.0-0.9)))^4/256/32*2^4*16*32
+            // (0.5^(-(1.0-0.9)))^4/256/32*2^4*256*4
+            // (0.5^(-(1.0-0.9)))^1/256/32*2^4*256*4
+            // (2^(-(1.0-0.9)))^4
+            // 16.0
+            .set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (/*64.0*/256.0/*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
+            // .set_persistence(/*0.9*/ /*2.0*/0.67)
+            // .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (2.0 * 32.0))
+            // .set_lacunarity(0.5)
             .set_seed(rng.gen());
+        // NOTE: octaves should definitely fit into i32, but we should check anyway to make
+        // sure.
+        /* assert!(rock_strength_nz.persistence > 0.0);
+        let rock_strength_scale = (1..rock_strength_nz.octaves as i32)
+            .map(|octave| rock_strength_nz.persistence.powi(octave + 1))
+            .sum::<f64>()
+            // For some reason, this is "scaled" by 3.0.
+            .mul(3.0);
+        let rock_strength_nz = ScaleBias::new(&rock_strength_nz)
+            .set_scale(1.0 / rock_strength_scale); */
 
-        let max_erosion_per_delta_t = 32.0 / CONFIG.mountain_scale as f64;
+        let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
+        let max_erosion_per_delta_t = 32.0 * height_scale;
         let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
         let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
         let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
-        let n_steps = 150; //150;//200;
+        let n_steps = 150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
+        let n_small_steps = 0;//8;//8; // 8
+
+        // fractal dimension should be between 0 and 0.9999...
+        // (lacunarity^octaves)^(-H) = persistence^(octaves)
+        // lacunarity^(octaves*-H) = persistence^(octaves)
+        // e^(-octaves*H*ln(lacunarity)) = e^(octaves * ln(persistence))
+        // -octaves * H * ln(lacunarity) = octaves * ln(persistence)
+        // -H = ln(persistence) / ln(lacunarity)
+        // H = -ln(persistence) / ln(lacunarity)
+        // ln(persistence) = -H * ln(lacunarity)
+        // persistence = lacunarity^(-H)
+        //
+        // -ln(2^(-0.25))/ln(2) = 0.25
+        //
+        // -ln(2^(-0.1))/ln(2)
+        //
+        // 0 = -ln(persistence) / ln(lacunarity)
+        // 0 = ln(persistence) => persistence = e^0 = 1
+        //
+        // 1 = -ln(persistence) / ln(lacunarity)
+        // -ln(lacunarity) = ln(persistence)
+        // e^(-ln(lacunarity)) = e^(ln(persistence))
+        // 1 / lacunarity = persistence
+        //
+        // Ergo, we should not set fractal dimension to anything  not between 1 / lacunarity and 1.
+        //
+        // dimension = -ln(0.25)/ln(2*pi/3) = 1.875
+        //
+        // (2*pi/3^1)^(-(-ln(0.25)/ln(2*pi/3))) = 0.25
+        //
+        // Default should be at most 1 / lacunarity.
+        //
+        // (2 * pi / 3)^(-ln(0.25)/ln(2*pi/3))
+        //
+        // -ln(0.25)/ln(2*pi/3) = 1.88
+        //
+        // (2 * pi / 3)^(-ln(0.25)/ln(2*pi/3))
+        //
+        // 2 * pi / 3
+        //
+        // 2.0^(2(-ln(1.5)/ln(2)))
+        // (1 / 1.5)^(2)
 
         // No NaNs in these uniform vectors, since the original noise value always returns Some.
         let ((alt_base, _), (chaos, _)) = rayon::join(
@@ -202,8 +276,13 @@ impl WorldSim {
                             .alt_nz
                             .get((wposf.div(10_000.0)).into_array())
                             .min(1.0)
-                            .max(-1.0))
-                        .sub(0.05)
+                            .max(-1.0)
+                            /* .mul(0.25)
+                            .add(0.125) */)
+                            // .add(0.5)
+                        // .sub(0.05)
+                        // .add(0.05)
+                        .add(0.075)
                         .mul(0.35), /*-0.0175*/
                     )
                 })
@@ -217,13 +296,13 @@ impl WorldSim {
                         //.add(0.0)
                         + gen_ctx
                             .hill_nz
-                            .get((wposf.div(1_500.0)).into_array())
+                            .get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(1_500.0)).into_array())
                             .min(1.0)
                             .max(-1.0)
                             .mul(1.0)
                         + gen_ctx
                             .hill_nz
-                            .get((wposf.div(400.0)).into_array())
+                            .get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(400.0)).into_array())
                             .min(1.0)
                             .max(-1.0)
                             .mul(0.3))
@@ -296,6 +375,7 @@ impl WorldSim {
                     .min(1.0)
                     .max(-1.0))
                 .abs()
+                // 0.5
                 .powf(1.35);
 
                 fn spring(x: f64, pow: f64) -> f64 {
@@ -305,14 +385,14 @@ impl WorldSim {
                 (0.0 + alt_main/*0.4*/
                     + (gen_ctx
                         .small_nz
-                        .get((wposf.div(300.0)).into_array())
+                        .get((wposf.mul(32.0).div(TerrainChunkSize::RECT_SIZE.map(|e| e as f64)).div(300.0)).into_array())
                         .min(1.0)
                         .max(-1.0))
                     .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
                     .mul(0.3)
                     .add(1.0)
                     .mul(0.4)
-                    /*0.52*/
+                    // 0.52
                     + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0).mul(0.045))
             };
 
@@ -357,7 +437,11 @@ impl WorldSim {
         });
 
         // Calculate oceans.
-        let old_height = |posi: usize| alt_old[posi].1;
+        let old_height = |posi: usize| alt_old[posi].1 * CONFIG.mountain_scale * height_scale as f32;
+        /* let is_ocean = (0..WORLD_SIZE.x * WORLD_SIZE.y)
+            .into_par_iter()
+            .map(|i| map_edge_factor(i) == 0.0)
+            .collect::<Vec<_>>(); */
         let is_ocean = get_oceans(old_height);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
 
@@ -406,6 +490,8 @@ impl WorldSim {
         let logistic_cdf = |x: f64| (x / logistic_2_base).tanh() * 0.5 + 0.5;
 
         let exp_inverse_cdf = |x: f64/*, pow: f64*/| -(-x).ln_1p()/* / ln(pow)*/;
+        // 2 / pi * ln(tan(pi/2 * p))
+        let hypsec_inverse_cdf = |x: f64| f64::consts::FRAC_2_PI * ((x * f64::consts::FRAC_PI_2).tan().ln());
         // 2^((2^10-2)/256) = 15.91...
         // -ln(1-(1-(2^(-22)*0.5)))
         // -ln(1-(1-(2^(-53)*0.5)))
@@ -423,8 +509,9 @@ impl WorldSim {
             1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64).max(f64::EPSILON as f64 * 0.5);
         let max_epsilon = (1.0 - 1.0 / (WORLD_SIZE.x as f64 * WORLD_SIZE.y as f64))
             .min(1.0 - f64::EPSILON as f64 * 0.5);
-        let alt_exp_min_uniform = exp_inverse_cdf(min_epsilon);
-        let alt_exp_max_uniform = exp_inverse_cdf(max_epsilon);
+        let inv_func = /*|x: f64| x*//*exp_inverse_cdf*/logit/*hypsec_inverse_cdf*/;
+        let alt_exp_min_uniform = /*exp_inverse_cdf*//*logit*/inv_func(min_epsilon);
+        let alt_exp_max_uniform = /*exp_inverse_cdf*//*logit*/inv_func(max_epsilon);
 
         // let erosion_pow = 2.0;
         // let n_steps = 100;//150;
@@ -437,58 +524,9 @@ impl WorldSim {
         };
         /* let erosion_factor = |x: f64| logistic_cdf(logistic_base * if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * log_odds(x))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
         ).powf(erosion_pow).mul(0.5))*/; */
-        let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(exp_inverse_cdf(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
-).powf(erosion_pow).mul(0.5))*/;
-        let alt = do_erosion(
-            0.0,
-            max_erosion_per_delta_t as f32,
-            n_steps,
-            &river_seed,
-            &rock_strength_nz,
-            |posi| {
-                if is_ocean_fn(posi) {
-                    old_height(posi)
-                } else {
-                    let wposf = (uniform_idx_as_vec2(posi)
-                        * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
-                    .map(|e| e as f64);
-                    let alt_main = {
-                        // Extension upwards from the base.  A positive number from 0 to 1 curved to be
-                        // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
-                        let alt_main = (gen_ctx
-                            .alt_nz
-                            .get((wposf.div(2_000.0)).into_array())
-                            .min(1.0)
-                            .max(-1.0))
-                        .abs()
-                        .powf(1.35);
-
-                        fn spring(x: f64, pow: f64) -> f64 {
-                            x.abs().powf(pow) * x.signum()
-                        }
-
-                        (0.0 + alt_main
-                            + (gen_ctx
-                                .small_nz
-                                .get((wposf.div(300.0)).into_array())
-                                .min(1.0)
-                                .max(-1.0))
-                            .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
-                            .mul(0.3)
-                            .add(1.0)
-                            .mul(0.4)
-                            + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
-                                .mul(0.045))
-                    };
-                    // old_height_uniform(posi) *
-                    (/*((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) **/(((6.0 / 360.0 * 2.0 * f64::consts::PI).tan()
-                        * TerrainChunkSize::RECT_SIZE.reduce_partial_min() as f64)
-                        .floor()
-                        / CONFIG.mountain_scale as f64)) as f32
-                    // 5.0 / CONFIG.mountain_scale
-                }
-            },
-            is_ocean_fn,
+        let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
+).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
+        let uplift_fn =
             |posi| {
                 if is_ocean_fn(posi) {
                     return 0.0;
@@ -521,8 +559,8 @@ impl WorldSim {
                         .mul(0.3)
                         .add(1.0)
                         .mul(0.4)
-                        + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
-                            .mul(0.045))
+                        /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
+                            .mul(0.045)*/)
                 };
                 let height =
                     ((old_height_uniform(posi) - alt_old_min_uniform) as f64
@@ -577,16 +615,98 @@ impl WorldSim {
                 } else {
                     (0.0, 0.0)
                 };
+                let height = 1.0f64;
+                // let height = 1.0 / 7.0f64;
                 let height = height
-                    .mul(max_erosion_per_delta_t * 7.0 / 8.0)
-                    .add(max_erosion_per_delta_t / 8.0)
+                    /* .mul(15.0 / 16.0)
+                    .add(1.0 / 16.0) */
+                    .mul(7.0 / 8.0)
+                    .add(1.0 / 8.0)
+                    .mul(max_erosion_per_delta_t)
                     .sub(/*1.0 / CONFIG.mountain_scale as f64*/ bump_max)
                     .add(bump_factor);
                 /* .sub(/*1.0 / CONFIG.mountain_scale as f64*/(f32::EPSILON * 0.5) as f64)
                 .add(bump_factor); */
                 height as f32
-            },
+            };
+        let alt_func = |posi| {
+            if is_ocean_fn(posi) {
+                // -max_erosion_per_delta_t as f32
+                // -1.0 / CONFIG.mountain_scale
+                // -0.75
+                // -CONFIG.sea_level / CONFIG.mountain_scale
+                // 0.0
+                old_height(posi) // 0.0
+            } else {
+                // uplift_fn(posi)
+                let wposf = (uniform_idx_as_vec2(posi)
+                    * TerrainChunkSize::RECT_SIZE.map(|e| e as i32))
+                .map(|e| e as f64);
+                let alt_main = {
+                    // Extension upwards from the base.  A positive number from 0 to 1 curved to be
+                    // maximal at 0.  Also to be multiplied by CONFIG.mountain_scale.
+                    let alt_main = (gen_ctx
+                        .alt_nz
+                        .get((wposf.div(2_000.0)).into_array())
+                        .min(1.0)
+                        .max(-1.0))
+                    .abs()
+                    .powf(1.35);
+
+                    fn spring(x: f64, pow: f64) -> f64 {
+                        x.abs().powf(pow) * x.signum()
+                    }
+
+                    (0.0 + alt_main
+                        + (gen_ctx
+                            .small_nz
+                            .get((wposf.div(300.0)).into_array())
+                            .min(1.0)
+                            .max(-1.0))
+                        .mul(alt_main.powf(0.8).max(/*0.25*/ 0.15))
+                        .mul(0.3)
+                        .add(1.0)
+                        .mul(0.4)
+                        /* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
+                            .mul(0.045)*/)
+                };
+                old_height_uniform(posi)
+                /* // 0.0
+                // -/*CONFIG.sea_level / CONFIG.mountain_scale*//* 0.75 */1.0
+                // ((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) as f32
+                // uplift_fn(posi) / max_erosion_per_delta_t as f32
+                // old_height_uniform(posi) *
+                (/*((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) **/(((6.0 / 360.0 * 2.0 * f64::consts::PI).tan()
+                    * TerrainChunkSize::RECT_SIZE.reduce_partial_min() as f64)
+                    .floor()
+                    * height_scale)) as f32
+                // 5.0 / CONFIG.mountain_scale */
+            }
+        };
+        let (alt, basement) = do_erosion(
+            0.0,
+            max_erosion_per_delta_t as f32,
+            n_steps,
+            &river_seed,
+            &rock_strength_nz,
+            |posi| alt_func(posi),
+            |posi| alt_func(posi),// if is_ocean_fn(posi) { old_height(posi) } else { 0.0 },
+            is_ocean_fn,
+            uplift_fn,
         );
+        // Quick "small scale" erosion cycle in order to lower extreme angles.
+        let (alt, basement) = do_erosion(
+            0.0,
+            (1.0 * height_scale) as f32,
+            n_small_steps,
+            &river_seed,
+            &rock_strength_nz,
+            |posi| /* if is_ocean_fn(posi) { old_height(posi) } else { alt[posi] } */alt[posi],
+            |posi| basement[posi],
+            is_ocean_fn,
+            |posi| uplift_fn(posi) * (1.0 * height_scale / max_erosion_per_delta_t) as f32,
+        );
+
         let is_ocean = get_oceans(|posi| alt[posi]);
         let is_ocean_fn = |posi: usize| is_ocean[posi];
         let mut dh = downhill(&alt, /*old_height*/ is_ocean_fn);
@@ -765,6 +885,7 @@ impl WorldSim {
             temp_base,
             chaos,
             alt,
+            basement,
             water_alt,
             dh,
             flux: flux_old,
@@ -786,7 +907,7 @@ impl WorldSim {
             rng,
         };
 
-        this.seed_elements();
+        // this.seed_elements();
 
         this
     }
@@ -1199,6 +1320,7 @@ impl WorldSim {
 pub struct SimChunk {
     pub chaos: f32,
     pub alt: f32,
+    pub basement: f32,
     pub water_alt: f32,
     pub downhill: Option<Vec2<i32>>,
     pub flux: f32,
@@ -1244,6 +1366,7 @@ impl SimChunk {
         let _map_edge_factor = map_edge_factor(posi);
         let (_, chaos) = gen_cdf.chaos[posi];
         let alt_pre = gen_cdf.alt[posi];
+        let basement_pre = gen_cdf.basement[posi];
         let water_alt_pre = gen_cdf.water_alt[posi];
         let downhill_pre = gen_cdf.dh[posi];
         let flux = gen_cdf.flux[posi];
@@ -1300,10 +1423,12 @@ impl SimChunk {
             panic!("Halp!");
         } */
 
-        let mut alt = CONFIG.sea_level.add(alt_pre.mul(CONFIG.mountain_scale));
+        let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale;
+        let mut alt = CONFIG.sea_level.add(alt_pre.div(height_scale));
+        let mut basement = CONFIG.sea_level.add(basement_pre.div(height_scale));
         let water_alt = CONFIG
             .sea_level
-            .add(water_alt_pre.mul(CONFIG.mountain_scale));
+            .add(water_alt_pre.div(height_scale));
         let downhill = if downhill_pre == -2 {
             None
         } else if downhill_pre < 0 {
@@ -1391,6 +1516,7 @@ impl SimChunk {
             chaos,
             flux,
             alt,
+            basement: basement.min(alt),
             water_alt,
             downhill,
             temp,
diff --git a/world/src/sim/util.rs b/world/src/sim/util.rs
index e7fc3aa7b7..b96ca14e99 100644
--- a/world/src/sim/util.rs
+++ b/world/src/sim/util.rs
@@ -1,9 +1,10 @@
 use super::WORLD_SIZE;
 use bitvec::prelude::{bitbox, bitvec, BitBox};
 use common::{terrain::TerrainChunkSize, vol::RectVolSize};
+use noise::{MultiFractal, NoiseFn, Perlin, Point2, Point3, Point4, Seedable};
 use num::Float;
 use rayon::prelude::*;
-use std::{f32, f64, u32};
+use std::{f32, f64, ops::Mul, u32};
 use vek::*;
 
 /// Calculates the smallest distance along an axis (x, y) from an edge of
@@ -301,3 +302,392 @@ pub fn get_oceans(oldh: impl Fn(usize) -> f32 + Sync) -> BitBox {
     }
     is_ocean
 }
+
+/// A 2-dimensional vector, for internal use.
+type Vector2<T> = [T; 2];
+/// A 3-dimensional vector, for internal use.
+type Vector3<T> = [T; 3];
+/// A 4-dimensional vector, for internal use.
+type Vector4<T> = [T; 4];
+
+#[inline]
+fn zip_with2<T, U, V, F>(a: Vector2<T>, b: Vector2<U>, f: F) -> Vector2<V>
+where
+    T: Copy,
+    U: Copy,
+    F: Fn(T, U) -> V,
+{
+    let (ax, ay) = (a[0], a[1]);
+    let (bx, by) = (b[0], b[1]);
+    [f(ax, bx), f(ay, by)]
+}
+
+#[inline]
+fn zip_with3<T, U, V, F>(a: Vector3<T>, b: Vector3<U>, f: F) -> Vector3<V>
+where
+    T: Copy,
+    U: Copy,
+    F: Fn(T, U) -> V,
+{
+    let (ax, ay, az) = (a[0], a[1], a[2]);
+    let (bx, by, bz) = (b[0], b[1], b[2]);
+    [f(ax, bx), f(ay, by), f(az, bz)]
+}
+
+#[inline]
+fn zip_with4<T, U, V, F>(a: Vector4<T>, b: Vector4<U>, f: F) -> Vector4<V>
+where
+    T: Copy,
+    U: Copy,
+    F: Fn(T, U) -> V,
+{
+    let (ax, ay, az, aw) = (a[0], a[1], a[2], a[3]);
+    let (bx, by, bz, bw) = (b[0], b[1], b[2], b[3]);
+    [f(ax, bx), f(ay, by), f(az, bz), f(aw, bw)]
+}
+
+#[inline]
+fn mul2<T>(a: Vector2<T>, b: T) -> Vector2<T>
+where
+    T: Copy + Mul<T, Output = T>,
+{
+    zip_with2(a, const2(b), Mul::mul)
+}
+
+#[inline]
+fn mul3<T>(a: Vector3<T>, b: T) -> Vector3<T>
+where
+    T: Copy + Mul<T, Output = T>,
+{
+    zip_with3(a, const3(b), Mul::mul)
+}
+
+#[inline]
+fn mul4<T>(a: Vector4<T>, b: T) -> Vector4<T>
+where
+    T: Copy + Mul<T, Output = T>,
+{
+    zip_with4(a, const4(b), Mul::mul)
+}
+
+#[inline]
+fn const2<T: Copy>(x: T) -> Vector2<T> {
+    [x, x]
+}
+
+#[inline]
+fn const3<T: Copy>(x: T) -> Vector3<T> {
+    [x, x, x]
+}
+
+#[inline]
+fn const4<T: Copy>(x: T) -> Vector4<T> {
+    [x, x, x, x]
+}
+
+fn build_sources(seed: u32, octaves: usize) -> Vec<Perlin> {
+    let mut sources = Vec::with_capacity(octaves);
+    for x in 0..octaves {
+        sources.push(Perlin::new().set_seed(seed + x as u32));
+    }
+    sources
+}
+
+/// Noise function that outputs hybrid Multifractal noise.
+///
+/// The result of this multifractal noise is that valleys in the noise should
+/// have smooth bottoms at all altitudes.
+#[derive(Clone, Debug)]
+pub struct HybridMulti {
+    /// Total number of frequency octaves to generate the noise with.
+    ///
+    /// The number of octaves control the _amount of detail_ in the noise
+    /// function. Adding more octaves increases the detail, with the drawback
+    /// of increasing the calculation time.
+    pub octaves: usize,
+
+    /// The number of cycles per unit length that the noise function outputs.
+    pub frequency: f64,
+
+    /// A multiplier that determines how quickly the frequency increases for
+    /// each successive octave in the noise function.
+    ///
+    /// The frequency of each successive octave is equal to the product of the
+    /// previous octave's frequency and the lacunarity value.
+    ///
+    /// A lacunarity of 2.0 results in the frequency doubling every octave. For
+    /// almost all cases, 2.0 is a good value to use.
+    pub lacunarity: f64,
+
+    /// A multiplier that determines how quickly the amplitudes diminish for
+    /// each successive octave in the noise function.
+    ///
+    /// The amplitude of each successive octave is equal to the product of the
+    /// previous octave's amplitude and the persistence value. Increasing the
+    /// persistence produces "rougher" noise.
+    ///
+    /// H = 1.0 - fractal increment = -ln(persistence) / ln(lacunarity).  For
+    /// a fractal increment between 0 (inclusive) and 1 (exclusive), keep
+    /// persistence between 1 / lacunarity (inclusive, for low fractal
+    /// dimension) and 1 (exclusive, for high fractal dimension).
+    pub persistence: f64,
+
+    /// An offset that is added to the output of each sample of the underlying
+    /// Perlin noise function.  Because each successive octave is weighted in
+    /// part by the previous signal's output, increasing the offset will weight
+    /// the output more heavily towards 1.0.
+    pub offset: f64,
+
+    seed: u32,
+    sources: Vec<Perlin>,
+}
+
+impl HybridMulti {
+    pub const DEFAULT_SEED: u32 = 0;
+    pub const DEFAULT_OCTAVES: usize = 6;
+    pub const DEFAULT_FREQUENCY: f64 = 2.0;
+    pub const DEFAULT_LACUNARITY: f64 = /* std::f64::consts::PI * 2.0 / 3.0 */2.0;
+    // -ln(2^(-0.25))/ln(2) = 0.25
+    // 2^(-0.25) ~ 13/16
+    pub const DEFAULT_PERSISTENCE: f64 = /* 0.25 *//* 0.5*/ 13.0 / 16.0;
+    pub const DEFAULT_OFFSET: f64 = /* 0.25 *//* 0.5*/ 0.7;
+    pub const MAX_OCTAVES: usize = 32;
+
+    pub fn new() -> Self {
+        Self {
+            seed: Self::DEFAULT_SEED,
+            octaves: Self::DEFAULT_OCTAVES,
+            frequency: Self::DEFAULT_FREQUENCY,
+            lacunarity: Self::DEFAULT_LACUNARITY,
+            persistence: Self::DEFAULT_PERSISTENCE,
+            offset: Self::DEFAULT_OFFSET,
+            sources: build_sources(Self::DEFAULT_SEED, Self::DEFAULT_OCTAVES),
+        }
+    }
+
+    pub fn set_offset(self, offset: f64) -> Self {
+        Self { offset, ..self }
+    }
+}
+
+impl Default for HybridMulti {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl MultiFractal for HybridMulti {
+    fn set_octaves(self, mut octaves: usize) -> Self {
+        if self.octaves == octaves {
+            return self;
+        }
+
+        octaves = octaves.max(1).min(Self::MAX_OCTAVES);
+        Self {
+            octaves,
+            sources: build_sources(self.seed, octaves),
+            ..self
+        }
+    }
+
+    fn set_frequency(self, frequency: f64) -> Self {
+        Self { frequency, ..self }
+    }
+
+    fn set_lacunarity(self, lacunarity: f64) -> Self {
+        Self { lacunarity, ..self }
+    }
+
+    fn set_persistence(self, persistence: f64) -> Self {
+        Self {
+            persistence,
+            ..self
+        }
+    }
+}
+
+impl Seedable for HybridMulti {
+    fn set_seed(self, seed: u32) -> Self {
+        if self.seed == seed {
+            return self;
+        }
+
+        Self {
+            seed,
+            sources: build_sources(seed, self.octaves),
+            ..self
+        }
+    }
+
+    fn seed(&self) -> u32 {
+        self.seed
+    }
+}
+
+/// 2-dimensional `HybridMulti` noise
+impl NoiseFn<Point2<f64>> for HybridMulti {
+    fn get(&self, mut point: Point2<f64>) -> f64 {
+        // First unscaled octave of function; later octaves are scaled.
+        point = mul2(point, self.frequency);
+        // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
+        let bias = 1.0;
+        let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
+        let mut scale = self.persistence;
+        let mut weight = result;
+
+        // Spectral construction inner loop, where the fractal is built.
+        for x in 1..self.octaves {
+            // Prevent divergence.
+            weight = weight.min(1.0);
+
+            // Raise the spatial frequency.
+            point = mul2(point, self.lacunarity);
+
+            // Get noise value, and scale it to the [offset - 1.0, 1.0 + offset] range.
+            let mut signal = (self.sources[x].get(point) + self.offset) * bias;
+
+            // Scale the amplitude appropriately for this frequency.
+            let exp_scale = self.persistence.powi(x as i32);
+            signal *= exp_scale;
+
+            // Add it in, weighted by previous octave's noise value.
+            result += weight * signal;
+
+            // Update the weighting value.
+            weight *= signal;
+            scale += exp_scale;
+        }
+
+        // Scale the result to the [-1,1] range
+        (result / scale) / bias - self.offset
+    }
+}
+
+/// 3-dimensional `HybridMulti` noise
+impl NoiseFn<Point3<f64>> for HybridMulti {
+    fn get(&self, mut point: Point3<f64>) -> f64 {
+        // First unscaled octave of function; later octaves are scaled.
+        point = mul3(point, self.frequency);
+        // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
+        let bias = 1.0;
+        let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
+        let mut scale = self.persistence;
+        let mut weight = result;
+
+        // Spectral construction inner loop, where the fractal is built.
+        for x in 1..self.octaves {
+            // Prevent divergence.
+            weight = weight.min(1.0);
+
+            // Raise the spatial frequency.
+            point = mul3(point, self.lacunarity);
+
+            // Get noise value, and scale it to the [0, 1.0] range.
+            let mut signal = (self.sources[x].get(point) + self.offset) * bias;
+
+            // Scale the amplitude appropriately for this frequency.
+            let exp_scale = self.persistence.powi(x as i32);
+            signal *= exp_scale;
+
+            // Add it in, weighted by previous octave's noise value.
+            result += weight * signal;
+
+            // Update the weighting value.
+            weight *= signal;
+            scale += exp_scale;
+        }
+
+        // Scale the result to the [-1,1] range
+        (result / scale) / bias - self.offset
+    }
+}
+
+/// 4-dimensional `HybridMulti` noise
+impl NoiseFn<Point4<f64>> for HybridMulti {
+    fn get(&self, mut point: Point4<f64>) -> f64 {
+        // First unscaled octave of function; later octaves are scaled.
+        point = mul4(point, self.frequency);
+        // Offset and bias to scale into [offset - 1.0, 1.0 + offset] range.
+        let bias = 1.0;
+        let mut result = (self.sources[0].get(point) + self.offset) * bias * self.persistence;
+        let mut exp_scale = self.persistence;
+        let mut scale = self.persistence;
+        let mut weight = result;
+
+        // Spectral construction inner loop, where the fractal is built.
+        for x in 1..self.octaves {
+            // Prevent divergence.
+            weight = weight.min(1.0);
+
+            // Raise the spatial frequency.
+            point = mul4(point, self.lacunarity);
+
+            // Get noise value, and scale it to the [0, 1.0] range.
+            let mut signal = (self.sources[x].get(point) + self.offset) * bias;
+
+            // Scale the amplitude appropriately for this frequency.
+            exp_scale *= self.persistence;
+            signal *= exp_scale;
+
+            // Add it in, weighted by previous octave's noise value.
+            result += weight * signal;
+
+            // Update the weighting value.
+            weight *= signal;
+            scale += exp_scale;
+        }
+
+        // Scale the result to the [-1,1] range
+        (result / scale) / bias - self.offset
+    }
+}
+
+/// Noise function that applies a scaling factor and a bias to the output value
+/// from the source function.
+///
+/// The function retrieves the output value from the source function, multiplies
+/// it with the scaling factor, adds the bias to it, then outputs the value.
+pub struct ScaleBias<'a, F: 'a> {
+    /// Outputs a value.
+    pub source: &'a F,
+
+    /// Scaling factor to apply to the output value from the source function.
+    /// The default value is 1.0.
+    pub scale: f64,
+
+    /// Bias to apply to the scaled output value from the source function.
+    /// The default value is 0.0.
+    pub bias: f64,
+}
+
+impl<'a, F> ScaleBias<'a, F> {
+    pub fn new(source: &'a F) -> Self {
+        ScaleBias {
+            source,
+            scale: 1.0,
+            bias: 0.0,
+        }
+    }
+
+    pub fn set_scale(self, scale: f64) -> Self {
+        ScaleBias { scale, ..self }
+    }
+
+    pub fn set_bias(self, bias: f64) -> Self {
+        ScaleBias { bias, ..self }
+    }
+}
+
+impl<'a, F: NoiseFn<T> + 'a, T> NoiseFn<T> for ScaleBias<'a, F> {
+    #[cfg(not(target_os = "emscripten"))]
+    fn get(&self, point: T) -> f64 {
+        (self.source.get(point)).mul_add(self.scale, self.bias)
+    }
+
+    #[cfg(target_os = "emscripten")]
+    fn get(&self, point: T) -> f64 {
+        (self.source.get(point) * self.scale) + self.bias
+    }
+}
+
+