Implement the vanilla Aquifer; stop flooding caves

Every air block below y=63 was turned into water. That is one line of code and
it cost the entire underground: no dry caves, no lava lakes, no air pockets, a
solid block of water from the sea floor to bedrock.

Vanilla decides fluid per position instead. Aquifer centres sit on a jittered
16x12x16 grid; each gets a fluid level and type from the floodedness and spread
noises, with centres near open sky inheriting the sea and buried ones getting a
much lower randomised level or nothing at all. A position takes its nearest
centre's fluid unless the barrier noise raises enough pressure between the two
or three nearest centres to seal it back to stone. Deep centres turn to lava.

Porting it means fixing the order of generation, not just adding a file. Vanilla
resolves stone/water/lava/air during the density pass and only then runs the
surface rules over a finished column; we did it the other way round, which is
what forced the unconditional flood in the first place. fillVanillaColumn now
asks the aquifer per position, and applySurfaceRule walks the finished column
carrying the bookkeeping SurfaceSystem carries: air resets the counters, a fluid
records its water height, and stone gets a depth from the top of its run plus
one from the bottom, found by looking ahead to the next non-stone block below.

That last one fixes stone_depth's ceiling form, which had no bottom-up depth to
work with and was testing the top-down one instead -- fourteen rules in the
overworld tree use it to dress cave roofs. The floor form is unchanged: vanilla
counts from 1 and compares against 1 + offset, we counted from 0 and compared
against offset.

The aquifer grid is built eagerly per chunk rather than lazily, because our
columns fill concurrently; every cell is a pure function of its grid coordinate
and every cell in the computed range gets consulted anyway. Cost is ~0.5% of
chunk generation, most of it absorbed by the shared preliminary-surface cache.

Inland caves go from 100% water to 3.8%, and lava exists for the first time.
cmd/gendump grows a census that would have failed loudly before, and
TestCavesAreDry guards it in the suite.

cmd/gendump/main.go

@@ -142,6 +142,75 @@ func main() {
 	}
 
 
+	// Caves are dry: the aquifer decides fluid per position, so the open volume
+	// underground is overwhelmingly air, with occasional aquifer pools and lava
+	// down low. The defect this catches is the old unconditional "flood every
+	// air block below sea level" pass, under which this number was 100%.
+	fmt.Println("\n=== Underground fluids: water fraction y=-50..40 over inland chunks, lava anywhere ===")
+	air, water, lava, solidU := 0, 0, 0, 0
+	deepLava := 0
+	inland := 0
+	for cx := int32(-12); cx <= 12; cx += 4 {
+		for cz := int32(-12); cz <= 12; cz += 4 {
+			ch := gen(cx, cz)
+			// Lava is counted everywhere; the water fraction only over land,
+			// since an ocean's water legitimately reaches its floor. Lava
+			// pockets cluster, so a narrow sample can miss them entirely.
+			land := isInland(ch)
+			if land {
+				inland++
+			}
+			for wy := world.MinY; wy <= 40; wy++ {
+				// The water fraction is measured over y=-50..40, above the band
+				// where the global fluid rule makes lava unconditional.
+				census := land && wy >= -50
+				for lx := 0; lx < 16; lx++ {
+					for lz := 0; lz < 16; lz++ {
+						switch ch.GetBlock(lx, wy, lz) {
+						case world.StateAir:
+							if census {
+								air++
+							}
+						case world.StateWater:
+							if census {
+								water++
+							}
+						case world.StateLava:
+							lava++
+							if wy < -54 {
+								deepLava++
+							}
+							if census {
+								air++ // open volume, just not water
+							}
+						default:
+							if census {
+								solidU++
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+	open := air + water
+	fmt.Printf("  chunks=%d solid=%d open=%d (air+lava=%d water=%d) | lava total=%d, of it below y=-54: %d\n",
+		inland, solidU, open, air, water, lava, deepLava)
+	switch {
+	case open == 0:
+		fmt.Println("  FAIL: no open volume underground at all")
+	default:
+		frac := float64(water) / float64(open)
+		fmt.Printf("  water is %.1f%% of the open volume\n", frac*100)
+		if frac > 0.35 {
+			fmt.Println("  FAIL: caves are flooded; the aquifer is not deciding fluid")
+		} else if lava == 0 {
+			fmt.Println("  FAIL: no lava anywhere underground")
+		} else {
+			fmt.Println("  OK: caves are dry and lava exists")
+		}
+	}
+
 	// Subsurface banding: find grass-topped land columns and print the top ~8
 	// blocks (grass cap → dirt band → stone) to confirm surfaceDepth widened the
 	// dirt band beyond a single block.
@@ -188,6 +257,28 @@ func loadName(od *worldgen.OverworldDensity, s2 worldgen.Sample2D, wx, wz int) s
 	return world.BiomeNameAt(od, wx, wz)
 }
 
+// isInland reports whether most of the chunk's columns break the surface above
+// sea level. Ocean chunks are excluded from the cave-fluid census because their
+// water legitimately reaches all the way down to the sea floor.
+func isInland(c *world.Chunk) bool {
+	aboveSea := 0
+	for lx := 0; lx < 16; lx += 2 {
+		for lz := 0; lz < 16; lz += 2 {
+			for wy := world.MinY + world.WorldHeight - 1; wy >= world.MinY; wy-- {
+				b := c.GetBlock(lx, wy, lz)
+				if b == world.StateAir {
+					continue
+				}
+				if b != world.StateWater && wy >= world.SeaLevel {
+					aboveSea++
+				}
+				break
+			}
+		}
+	}
+	return aboveSea > 48 // of 64 sampled columns
+}
+
 func crossSection(c *world.Chunk) {
 	// vertical band from y=40..136
 	for wy := 130; wy >= 40; wy-- {
@@ -205,6 +296,8 @@ func glyph(b uint16) string {
 		return "."
 	case world.StateWater:
 		return "~"
+	case world.StateLava:
+		return "!"
 	case world.StateStone:
 		return "#"
 	case world.StateDirt:

internal/world/aquifer_verify_test.go

@@ -0,0 +1,112 @@
+package world
+
+import "testing"
+
+// TestCavesAreDry is the load-bearing check for the aquifer. Before it landed,
+// every air block below sea level was turned into water unconditionally, so
+// every cave under y=63 was a solid block of water and no lava existed
+// anywhere. The aquifer decides fluid per position instead, and the visible
+// consequence is that inland caves are overwhelmingly air.
+//
+// The thresholds are deliberately loose — this is a regression guard against
+// the whole underground filling up again, not a parity check.
+func TestCavesAreDry(t *testing.T) {
+	gen := NewVanillaGenerator(12345)
+	air, water, lava := 0, 0, 0
+	inland := 0
+	// A wide grid rather than a handful of chunks: lava pockets are clustered,
+	// so a small sample can legitimately contain none.
+	for cx := int32(-12); cx <= 12; cx += 4 {
+		for cz := int32(-12); cz <= 12; cz += 4 {
+			ch := gen(cx, cz)
+			if ch == nil {
+				continue
+			}
+			// Lava counts everywhere; the water fraction only over land, since
+			// an ocean's water legitimately reaches its floor.
+			land := inlandChunk(ch)
+			if land {
+				inland++
+			}
+			for wy := MinY; wy <= 40; wy++ {
+				census := land && wy >= -50
+				for lx := 0; lx < 16; lx++ {
+					for lz := 0; lz < 16; lz++ {
+						switch ch.GetBlock(lx, wy, lz) {
+						case StateAir:
+							if census {
+								air++
+							}
+						case StateWater:
+							if census {
+								water++
+							}
+						case StateLava:
+							lava++
+							if census {
+								air++
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+	if inland == 0 {
+		t.Skip("no inland chunks in the scanned area")
+	}
+	open := air + water
+	if open == 0 {
+		t.Fatalf("no open volume underground across %d inland chunks", inland)
+	}
+	frac := float64(water) / float64(open)
+	t.Logf("inland chunks=%d open=%d water=%d (%.1f%%) lava=%d", inland, open, water, frac*100, lava)
+	if frac > 0.35 {
+		t.Errorf("water is %.1f%% of the open volume below ground; caves are flooded", frac*100)
+	}
+	if lava == 0 {
+		t.Error("no lava underground: the aquifer never picks a lava fluid type")
+	}
+}
+
+// TestNoFluidUnderBedrock guards the world floor: the aquifer runs all the way
+// down, and a fluid level reaching below y=-59 would put water or lava inside
+// the bedrock band.
+func TestNoFluidUnderBedrock(t *testing.T) {
+	gen := NewVanillaGenerator(12345)
+	for _, p := range [][2]int32{{0, 0}, {5, -7}, {-13, 21}} {
+		ch := gen(p[0], p[1])
+		for wy := MinY; wy <= MinY+5; wy++ {
+			for lx := 0; lx < 16; lx++ {
+				for lz := 0; lz < 16; lz++ {
+					switch b := ch.GetBlock(lx, wy, lz); b {
+					case StateAir, StateWater, StateLava:
+						t.Fatalf("chunk(%d,%d) block %d at (%d,%d,%d) inside the bedrock band", p[0], p[1], b, lx, wy, lz)
+					}
+				}
+			}
+		}
+	}
+}
+
+// inlandChunk reports whether most of the chunk breaks the surface above sea
+// level. Ocean chunks are excluded from the cave census: their water reaches
+// the sea floor legitimately.
+func inlandChunk(c *Chunk) bool {
+	aboveSea := 0
+	for lx := 0; lx < 16; lx += 2 {
+		for lz := 0; lz < 16; lz += 2 {
+			for wy := MinY + WorldHeight - 1; wy >= MinY; wy-- {
+				b := c.GetBlock(lx, wy, lz)
+				if b == StateAir {
+					continue
+				}
+				if b != StateWater && wy >= SeaLevel {
+					aboveSea++
+				}
+				break
+			}
+		}
+	}
+	return aboveSea > 48 // of 64 sampled columns
+}

internal/world/encode.go

@@ -26,6 +26,7 @@ const (
 	StateDirt    uint16 = 10
 	StateBedrock uint16 = 85
 	StateWater   uint16 = 86
+	StateLava    uint16 = 102
 	StateSand    uint16 = 118
 	StateGravel  uint16 = 124
 	StateOakLog  uint16 = 137

internal/world/store.go

@@ -32,7 +32,7 @@ const dataVersion26 = 4790
 // first time it ran: chunkAt prefers the store over the generator, so the
 // already-explored area around spawn keeps its old terrain and every later fix
 // looks like it did nothing in exactly the place you are standing.
-const generatorVersion = 1
+const generatorVersion = 2
 
 // generatorVersionTag is the NBT key holding generatorVersion. It is namespaced
 // because it is ours, not part of the vanilla chunk format.

internal/world/vanilla.go

@@ -1,6 +1,7 @@
 package world
 
 import (
+	"math"
 	"math/rand"
 	"sync"
 
@@ -29,12 +30,13 @@ func NewVanillaGenerator(seed int64) Generator {
 	if err != nil {
 		panic("world: loading overworld density: " + err.Error())
 	}
+	fluidPicker := worldgen.OverworldFluidPicker(od.SeaLevel)
 	return func(cx, cz int32) *Chunk {
-		return generateVanilla(od, seed, cx, cz)
+		return generateVanilla(od, fluidPicker, seed, cx, cz)
 	}
 }
 
-func generateVanilla(od *worldgen.OverworldDensity, seed int64, cx, cz int32) *Chunk {
+func generateVanilla(od *worldgen.OverworldDensity, fluidPicker worldgen.FluidPicker, seed int64, cx, cz int32) *Chunk {
 	c := NewChunk(cx, cz, BiomePlains) // per-cell biomes override below
 	baseX, baseZ := int(cx)*16, int(cz)*16
 
@@ -81,6 +83,14 @@ func generateVanilla(od *worldgen.OverworldDensity, seed int64, cx, cz int32) *C
 	// to parse, surface fill falls back to the biome-blind heuristics.
 	surfaceRule, ruleErr := od.SurfaceRule()
 
+	// The aquifer decides fluid per position while the column is laid down. Its
+	// cell grid spans the chunk plus a margin, so it is built once per chunk and
+	// shared, read-only, by the parallel column fill.
+	var aq *worldgen.Aquifer
+	if od.AquifersEnabled {
+		aq = worldgen.NewAquifer(od, int(cx), int(cz), fluidPicker)
+	}
+
 	var columns [16][16][WorldHeight]uint16
 	var surfTop [16][16]int // top solid index, -1 if none
 	var grass [16][16]bool  // grassy land surface (tree-plantable)
@@ -94,7 +104,7 @@ func generateVanilla(od *worldgen.OverworldDensity, seed int64, cx, cz int32) *C
 				if ruleErr == nil {
 					rule = surfaceRule
 				}
-				surfTop[lx][lz], grass[lx][lz] = fillVanillaColumn(od, grids, interp, &columns[lx][lz], baseX+lx, baseZ+lz, lx, lz, seed, rule, biomeName[lx][lz])
+				surfTop[lx][lz], grass[lx][lz] = fillVanillaColumn(od, aq, fluidPicker, grids, interp, &columns[lx][lz], baseX+lx, baseZ+lz, lx, lz, seed, rule, biomeName[lx][lz])
 			}
 		}(lx)
 	}
@@ -147,17 +157,20 @@ func fillBiomes3D(c *Chunk, od *worldgen.OverworldDensity, s2D [16][16]worldgen.
 }
 
 // fillVanillaColumn lays the blocks for one column and returns the top solid
-// index and whether the surface is grassy land (suitable for trees). When a
+// index and whether the surface is grassy land (suitable for trees).
-// surface rule tree is provided, surface blocks are decided by it (vanilla
+//
-// behaviour: biome/depth/steepness/water/y-driven); otherwise the legacy
+// The order matches vanilla: the density pass decides stone-or-not, the aquifer
-// beach/grass/dirt heuristics are used as a fallback.
+// turns every non-stone position into air, water or lava (and can also seal a
-func fillVanillaColumn(od *worldgen.OverworldDensity, grids []cornerGrid, interp []float64, out *[WorldHeight]uint16, wx, wz, lx, lz int, seed int64, rule worldgen.SurfaceRule, biomeName string) (int, bool) {
+// position back to stone where the barrier noise says the rock holds), and only
+// then does the surface rule tree walk the finished column. Doing it the other
+// way round is what forced the old unconditional "flood everything under sea
+// level" pass, which left every cave below y=63 underwater.
+func fillVanillaColumn(od *worldgen.OverworldDensity, aq *worldgen.Aquifer, fluidPicker worldgen.FluidPicker, grids []cornerGrid, interp []float64, out *[WorldHeight]uint16, wx, wz, lx, lz int, seed int64, rule worldgen.SurfaceRule, biomeName string) (int, bool) {
 	cx0 := lx / cellWidth
 	cz0 := lz / cellWidth
 	fx := float64(lx%cellWidth) / cellWidth
 	fz := float64(lz%cellWidth) / cellWidth
 
-	var solid [WorldHeight]bool
 	top := -1
 	for i := 0; i < WorldHeight; i++ {
 		cy0 := i / cellHeight
@@ -165,9 +178,12 @@ func fillVanillaColumn(od *worldgen.OverworldDensity, grids []cornerGrid, interp
 		for n := range grids {
 			interp[n] = trilerp(&grids[n], cx0, cy0, cz0, fx, fy, fz)
 		}
-		ctx := worldgen.FunctionContext{X: float64(wx), Y: float64(MinY + i), Z: float64(wz)}.WithInterp(interp)
+		y := MinY + i
-		if od.Final.Compute(ctx) > 0 {
+		ctx := worldgen.FunctionContext{X: float64(wx), Y: float64(y), Z: float64(wz)}.WithInterp(interp)
-			solid[i] = true
+		density := od.Final.Compute(ctx)
+		state, isDefaultBlock := substance(aq, fluidPicker, wx, y, wz, density)
+		out[i] = state
+		if isDefaultBlock {
 			top = i
 		}
 	}
@@ -183,28 +199,56 @@ func fillVanillaColumn(od *worldgen.OverworldDensity, grids []cornerGrid, interp
 	rng := newColumnRand(wx, wz, int(seed))
 
 	if rule != nil {
-		applySurfaceRule(out, solid, top, wx, wz, SeaLevel, MinY, biomeName, rule, rng)
+		applySurfaceRule(out, wx, wz, SeaLevel, MinY, biomeName, rule, rng, top)
 	} else {
-		fillLegacySurface(out, solid, top, beach, deepWater, topY, rng)
+		fillLegacySurface(out, top, beach, deepWater, rng)
-	}
-	// Water fills air below sea level regardless of rule path.
-	for i := 0; i < WorldHeight; i++ {
-		if out[i] == StateAir && MinY+i < SeaLevel {
-			out[i] = StateWater
-		}
 	}
 	return top, top >= 0 && !beach && !deepWater && topY >= SeaLevel
 }
 
-// applySurfaceRule walks the column top-to-surface applying the rule tree. For
+// substance resolves one position to the block the terrain pass leaves behind:
-// each solid block it builds a SurfaceContext and lets the rule decide; the
+// the default block where the density is solid, otherwise whatever the aquifer
-// stone depth counts how far below the surface the block sits. Air blocks
+// puts there — air, water or lava. The second result says which of the two
-// above the surface are left for the water fill.
+// happened, so the caller can track the top solid block without re-testing.
+func substance(aq *worldgen.Aquifer, fluidPicker worldgen.FluidPicker, x, y, z int, density float64) (state uint16, isDefaultBlock bool) {
+	if aq == nil {
+		// aquifers_enabled=false: Aquifer.createDisabled, the global fluid rule
+		// with no cells and no barriers.
+		if density > 0 {
+			return StateStone, true
+		}
+		return fluidPicker(x, y, z).At(y), false
+	}
+	if s, ok := aq.ComputeSubstance(x, y, z, density); ok {
+		return s, false
+	}
+	return StateStone, true
+}
+
+// applySurfaceRule walks the finished column from the top down, applying the
+// rule tree to every default-block position, and mirrors SurfaceSystem's
+// bookkeeping as it goes:
+//
+//   - air resets both the stone depth and the water height;
+//   - a fluid records the height of the first (topmost) block of its run;
+//   - stone carries a depth counted down from the top of its run, and a depth
+//     counted up from the bottom, found by looking ahead to the next non-stone
+//     block below.
+//
+// The rule only replaces the default block, so anything the aquifer placed —
+// water in an ocean, lava in a deep pocket — survives untouched.
 //
 // One *rand.Rand is created per column (not per block) — bandlands/gradient
 // consume from it sequentially, which is correct because vanilla seeds those
 // per-column too. This avoids ~98k rand.New allocations per chunk.
-func applySurfaceRule(out *[WorldHeight]uint16, solid [WorldHeight]bool, top int, wx, wz, seaLevel, minY int, biomeName string, rule worldgen.SurfaceRule, rng chunkRand) {
+func applySurfaceRule(out *[WorldHeight]uint16, wx, wz, seaLevel, minY int, biomeName string, rule worldgen.SurfaceRule, rng chunkRand, topSolid int) {
+	top := -1
+	for i := WorldHeight - 1; i >= 0; i-- {
+		if out[i] != StateAir {
+			top = i
+			break
+		}
+	}
 	if top < 0 {
 		return
 	}
@@ -225,52 +269,82 @@ func applySurfaceRule(out *[WorldHeight]uint16, solid [WorldHeight]bool, top int
 		MinY:               minY,
 		SurfaceNoise:       surfaceNoise,
 		SurfaceDepth:       0,
-		PreliminarySurface: minY + top,
+		PreliminarySurface: minY + topSolid,
 		Rng:                colRng,
 	}
+	stoneDepthAbove := 0
+	waterHeight := math.MinInt
+	nextCeilingStoneY := math.MaxInt
 	for i := top; i >= 0; i-- {
-		if !solid[i] {
+		y := minY + i
+		old := out[i]
+		if old == StateAir {
+			stoneDepthAbove = 0
+			waterHeight = math.MinInt
+			continue
+		}
+		if isFluidState(old) {
+			if waterHeight == math.MinInt {
+				waterHeight = y + 1
+			}
+			continue
+		}
+		if nextCeilingStoneY >= y {
+			// Look ahead to the first non-stone block below; the scan runs one
+			// past the world floor, which reads as air, so it always terminates.
+			nextCeilingStoneY = worldgen.WayBelowMinY
+			for j := i - 1; j >= -1; j-- {
+				if j >= 0 && isStoneState(out[j]) {
+					continue
+				}
+				nextCeilingStoneY = minY + j + 1
+				break
+			}
+		}
+		stoneDepthAbove++
+		sctx.Y = y
+		sctx.StoneDepthAbove = stoneDepthAbove
+		sctx.StoneDepthBelow = y - nextCeilingStoneY + 1
+		sctx.WaterHeight = waterHeight
+		if old != StateStone {
 			continue
 		}
-		sctx.Y = minY + i
-		sctx.StoneDepthAbove = top - i
-		// Solid blocks default to stone; the rule tree overrides only the
-		// surface layers it matches (grass/sand/terracotta/etc). Blocks where
-		// the rule does not match (depth > surface band) keep stone, matching
-		// vanilla: surface rules replace only the top few blocks, the column is
-		// otherwise stone down to bedrock.
-		out[i] = StateStone
 		if state, ok := rule.Apply(sctx); ok && state != 0 {
 			out[i] = state
 		}
 	}
 }
 
+// isFluidState reports whether a raw terrain block is a fluid (SurfaceSystem
+// branches on getFluidState().isEmpty()). Only the aquifer's own fluids can
+// appear here, since the rule pass runs before decoration.
+func isFluidState(s uint16) bool { return s == StateWater || s == StateLava }
+
+// isStoneState is SurfaceSystem.isStone: solid, non-fluid, non-air.
+func isStoneState(s uint16) bool { return s != StateAir && !isFluidState(s) }
+
 // fillLegacySurface is the biome-blind heuristic used when no surface rule is
-// available (parse failure). It mirrors the pre-surface-rule block switch.
+// available (parse failure). It dresses the stone the terrain and aquifer
-func fillLegacySurface(out *[WorldHeight]uint16, solid [WorldHeight]bool, top int, beach, deepWater bool, topY int, rng chunkRand) {
+// passes already laid down, leaving their air and fluids alone.
+func fillLegacySurface(out *[WorldHeight]uint16, top int, beach, deepWater bool, rng chunkRand) {
 	for i := 0; i < WorldHeight; i++ {
 		y := MinY + i
+		if !isStoneState(out[i]) {
+			continue
+		}
 		switch {
 		case y <= MinY:
 			out[i] = StateBedrock
-		case y <= MinY+4 && solid[i] && bedrockAt(&rng, y-MinY):
+		case y <= MinY+4 && bedrockAt(&rng, y-MinY):
 			out[i] = StateBedrock
-		case solid[i]:
+		case beach && i > top-4:
-			switch {
+			out[i] = StateSand
-			case beach && i > top-4:
+		case deepWater && i == top:
-				out[i] = StateSand
+			out[i] = StateGravel
-			case deepWater && i == top:
+		case i == top && y >= SeaLevel:
-				out[i] = StateGravel
+			out[i] = StateGrass
-			case i == top && y >= SeaLevel:
+		case i > top-4:
-				out[i] = StateGrass
+			out[i] = StateDirt
-			case i > top-4:
-				out[i] = StateDirt
-			default:
-				out[i] = StateStone
-			}
-		case y < SeaLevel:
-			out[i] = StateWater
 		}
 	}
 }

internal/worldgen/aquifer.go

@@ -0,0 +1,446 @@
+package worldgen
+
+import (
+	"math"
+	"sync"
+)
+
+// aquifer.go ports net.minecraft.world.level.levelgen.Aquifer.NoiseBasedAquifer.
+//
+// The aquifer is what decides, for every position the density function leaves
+// empty, whether it becomes air, water or lava. Without it a generator has to
+// guess — the usual guess being "water everywhere below sea level", which
+// drowns every cave under y=63 and leaves no lava lakes anywhere.
+//
+// Vanilla scatters aquifer centres on a 16×12×16 grid, jittered by a positional
+// RNG. Each centre gets a FluidStatus: a fluid level and a fluid type. A
+// position takes the fluid of its nearest centre, unless the barrier noise
+// raises enough "pressure" between the two or three nearest centres to seal the
+// position off as stone instead. Centres near the open sky inherit the global
+// sea level, so oceans and lakes still fill normally; centres buried deep get a
+// randomised, usually much lower level, which is why caves are dry.
+
+// Block-state network IDs the aquifer places. The worldgen package deliberately
+// does not import the world package; these match blockids.go.
+const (
+	blockAir   uint16 = 0
+	blockWater uint16 = 86
+	blockLava  uint16 = 102
+)
+
+// Aquifer grid geometry (Aquifer.NoiseBasedAquifer constants).
+const (
+	aquiferXSpacing = 16
+	aquiferYSpacing = 12
+	aquiferZSpacing = 16
+	aquiferXRange   = 10
+	aquiferYRange   = 9
+	aquiferZRange   = 10
+
+	// WayBelowMinY is DimensionType.WAY_BELOW_MIN_Y (MIN_Y << 4, MIN_Y=-2032):
+	// the "this aquifer holds nothing" sentinel fluid level.
+	WayBelowMinY = -32512
+)
+
+// deepDark is OverworldBiomeBuilder.isDeepDarkRegion's thresholds, kept at the
+// exact double values the float constants widen to.
+const (
+	deepDarkErosionMax = -0.22499999403953552
+	deepDarkDepthMin   = 0.8999999761581421
+)
+
+// FluidStatus is a fluid level plus the fluid filling up to it (Aquifer.FluidStatus).
+type FluidStatus struct {
+	Level int
+	Type  uint16
+}
+
+// At returns the fluid at blockY, or air above the level.
+func (f FluidStatus) At(blockY int) uint16 {
+	if blockY < f.Level {
+		return f.Type
+	}
+	return blockAir
+}
+
+// FluidPicker is the dimension-wide fluid rule (Aquifer.FluidPicker): what a
+// position would hold if there were no aquifer at all.
+type FluidPicker func(x, y, z int) FluidStatus
+
+// OverworldFluidPicker is NoiseBasedChunkGenerator.createFluidPicker: lava
+// below y=-54, sea water above it.
+func OverworldFluidPicker(seaLevel int) FluidPicker {
+	lava := FluidStatus{Level: -54, Type: blockLava}
+	sea := FluidStatus{Level: seaLevel, Type: blockWater}
+	lavaBelow := min(-54, seaLevel)
+	return func(_, y, _ int) FluidStatus {
+		if y < lavaBelow {
+			return lava
+		}
+		return sea
+	}
+}
+
+// surfaceSamplingOffsets is SURFACE_SAMPLING_OFFSETS_IN_CHUNKS: the thirteen
+// chunk offsets an aquifer centre probes to work out whether it is under open
+// sky or buried. The set is lopsided towards -X on purpose — it is vanilla's.
+var surfaceSamplingOffsets = [13][2]int{
+	{0, 0}, {-2, -1}, {-1, -1}, {0, -1}, {1, -1}, {-3, 0}, {-2, 0},
+	{-1, 0}, {1, 0}, {-2, 1}, {-1, 1}, {0, 1}, {1, 1},
+}
+
+// Aquifer resolves fluid for one chunk. Its cell grid is computed up front so
+// the chunk's columns can be filled in parallel without locking.
+type Aquifer struct {
+	od     *OverworldDensity
+	global FluidPicker
+
+	minGridX, minGridY, minGridZ    int
+	gridSizeX, gridSizeY, gridSizeZ int
+
+	locations []aquiferPos
+	status    []FluidStatus
+
+	// skipSamplingAboveY is the height above which the grid is irrelevant and
+	// the global fluid rule answers directly.
+	skipSamplingAboveY int
+}
+
+type aquiferPos struct{ x, y, z int }
+
+// NewAquifer builds the aquifer covering the given chunk.
+//
+// Vanilla fills the cell grid lazily as columns are generated; we fill it
+// eagerly because our columns are generated concurrently. That is not a
+// fidelity change: every cell's centre and status is a pure function of its
+// grid coordinate, and every cell in the range computed here is consulted by
+// some position in the chunk anyway.
+func NewAquifer(od *OverworldDensity, chunkX, chunkZ int, picker FluidPicker) *Aquifer {
+	minBlockX, minBlockZ := chunkX*16, chunkZ*16
+	maxBlockX, maxBlockZ := minBlockX+15, minBlockZ+15
+
+	a := &Aquifer{od: od, global: picker}
+	a.minGridX = aquiferGridX(minBlockX - 5)
+	maxGridX := aquiferGridX(maxBlockX-5) + 1
+	a.gridSizeX = maxGridX - a.minGridX + 1
+	a.minGridY = aquiferGridY(od.MinY+1) - 1
+	maxGridY := aquiferGridY(od.MinY+od.Height+1) + 1
+	a.gridSizeY = maxGridY - a.minGridY + 1
+	a.minGridZ = aquiferGridZ(minBlockZ - 5)
+	maxGridZ := aquiferGridZ(maxBlockZ-5) + 1
+	a.gridSizeZ = maxGridZ - a.minGridZ + 1
+
+	n := a.gridSizeX * a.gridSizeY * a.gridSizeZ
+	a.locations = make([]aquiferPos, n)
+	a.status = make([]FluidStatus, n)
+
+	maxAdjusted := adjustSurfaceLevel(od.MaxPreliminarySurfaceLevel(
+		fromAquiferGridX(a.minGridX, 0), fromAquiferGridZ(a.minGridZ, 0),
+		fromAquiferGridX(maxGridX, 9), fromAquiferGridZ(maxGridZ, 9)))
+	a.skipSamplingAboveY = fromAquiferGridY(aquiferGridY(maxAdjusted+12)+1, 11) - 1
+
+	// Cells above the highest consulted anchor are never read: computeSubstance
+	// returns the global fluid before touching the grid once y climbs past
+	// skipSamplingAboveY, and the anchor search reaches at most one cell higher.
+	topUsedGridY := min(aquiferGridY(a.skipSamplingAboveY+1)+1, maxGridY)
+
+	var wg sync.WaitGroup
+	for gy := a.minGridY; gy <= topUsedGridY; gy++ {
+		wg.Add(1)
+		go func(gy int) {
+			defer wg.Done()
+			for gz := a.minGridZ; gz < a.minGridZ+a.gridSizeZ; gz++ {
+				for gx := a.minGridX; gx < a.minGridX+a.gridSizeX; gx++ {
+					i := a.index(gx, gy, gz)
+					r := od.AquiferRandom.At(gx, gy, gz)
+					pos := aquiferPos{
+						x: fromAquiferGridX(gx, int(r.NextIntN(aquiferXRange))),
+						y: fromAquiferGridY(gy, int(r.NextIntN(aquiferYRange))),
+						z: fromAquiferGridZ(gz, int(r.NextIntN(aquiferZRange))),
+					}
+					a.locations[i] = pos
+					a.status[i] = a.computeFluid(pos.x, pos.y, pos.z)
+				}
+			}
+		}(gy)
+	}
+	wg.Wait()
+	return a
+}
+
+func (a *Aquifer) index(gridX, gridY, gridZ int) int {
+	x := gridX - a.minGridX
+	y := gridY - a.minGridY
+	z := gridZ - a.minGridZ
+	return (y*a.gridSizeZ+z)*a.gridSizeX + x
+}
+
+// ComputeSubstance decides what fills (x,y,z) given the final density there.
+// ok=false means the position stays the settings' default block (stone);
+// otherwise the returned state is the fluid — which may be air.
+//
+// Vanilla additionally tracks shouldScheduleFluidUpdate here, to mark positions
+// where two neighbouring aquifers disagree so the fluid flows on first tick. We
+// have no fluid ticking yet and the flag never affects the block placed, so it
+// is left out; the fourth-nearest centre, which only feeds that flag, is not
+// tracked either.
+func (a *Aquifer) ComputeSubstance(x, y, z int, density float64) (uint16, bool) {
+	if density > 0 {
+		return 0, false
+	}
+	global := a.global(x, y, z)
+	if y > a.skipSamplingAboveY {
+		return global.At(y), true
+	}
+	if global.At(y) == blockLava {
+		return blockLava, true
+	}
+
+	xAnchor := aquiferGridX(x - 5)
+	yAnchor := aquiferGridY(y + 1)
+	zAnchor := aquiferGridZ(z - 5)
+	dist1, dist2, dist3 := math.MaxInt32, math.MaxInt32, math.MaxInt32
+	idx1, idx2, idx3 := 0, 0, 0
+	for dx := 0; dx <= 1; dx++ {
+		for dy := -1; dy <= 1; dy++ {
+			for dz := 0; dz <= 1; dz++ {
+				i := a.index(xAnchor+dx, yAnchor+dy, zAnchor+dz)
+				p := a.locations[i]
+				ox, oy, oz := p.x-x, p.y-y, p.z-z
+				d := ox*ox + oy*oy + oz*oz
+				switch {
+				case dist1 >= d:
+					idx3, idx2, idx1 = idx2, idx1, i
+					dist3, dist2, dist1 = dist2, dist1, d
+				case dist2 >= d:
+					idx3, idx2 = idx2, i
+					dist3, dist2 = dist2, d
+				case dist3 >= d:
+					idx3, dist3 = i, d
+				}
+			}
+		}
+	}
+
+	closest1 := a.status[idx1]
+	sim12 := aquiferSimilarity(dist1, dist2)
+	fluid := closest1.At(y)
+	if sim12 <= 0 {
+		return fluid, true
+	}
+	// Water sitting directly on the global lava level always wins: it is what
+	// makes the lava-lake shorelines steam rather than vanish.
+	if fluid == blockWater && a.global(x, y-1, z).At(y-1) == blockLava {
+		return fluid, true
+	}
+
+	barrierNoise := math.NaN()
+	closest2 := a.status[idx2]
+	if density+sim12*a.calculatePressure(x, y, z, &barrierNoise, closest1, closest2) > 0 {
+		return 0, false
+	}
+	closest3 := a.status[idx3]
+	if sim13 := aquiferSimilarity(dist1, dist3); sim13 > 0 {
+		if density+sim12*sim13*a.calculatePressure(x, y, z, &barrierNoise, closest1, closest3) > 0 {
+			return 0, false
+		}
+	}
+	if sim23 := aquiferSimilarity(dist2, dist3); sim23 > 0 {
+		if density+sim12*sim23*a.calculatePressure(x, y, z, &barrierNoise, closest2, closest3) > 0 {
+			return 0, false
+		}
+	}
+	return fluid, true
+}
+
+// aquiferSimilarity falls from 1 to 0 as the second distance pulls away from
+// the first; at or below 0 the nearest centre wins outright and no barrier is
+// evaluated.
+func aquiferSimilarity(distSqr1, distSqr2 int) float64 {
+	return 1.0 - float64(distSqr2-distSqr1)/25.0
+}
+
+// calculatePressure is the barrier between two aquifers: how hard the rock
+// between them resists being carved open. barrierNoise memoises the noise
+// sample across the (up to three) pressure evaluations at one position, exactly
+// as vanilla's MutableDouble does.
+func (a *Aquifer) calculatePressure(x, y, z int, barrierNoise *float64, s1, s2 FluidStatus) float64 {
+	type1 := s1.At(y)
+	type2 := s2.At(y)
+	if (type1 == blockLava && type2 == blockWater) || (type1 == blockWater && type2 == blockLava) {
+		return 2.0
+	}
+	fluidYDiff := s1.Level - s2.Level
+	if fluidYDiff < 0 {
+		fluidYDiff = -fluidYDiff
+	}
+	if fluidYDiff == 0 {
+		return 0.0
+	}
+	averageFluidY := 0.5 * float64(s1.Level+s2.Level)
+	howFarAboveAverage := float64(y) + 0.5 - averageFluidY
+	baseValue := float64(fluidYDiff) / 2.0
+	// Distance from the barrier's edge towards its middle; the biases below are
+	// vanilla's, and they are asymmetric: rock reaches much further down from a
+	// fluid surface than up from it.
+	distanceFromEdge := baseValue - math.Abs(howFarAboveAverage)
+	var gradient float64
+	if howFarAboveAverage > 0 {
+		if centerPoint := 0.0 + distanceFromEdge; centerPoint > 0 {
+			gradient = centerPoint / 1.5
+		} else {
+			gradient = centerPoint / 2.5
+		}
+	} else {
+		if centerPoint := 3.0 + distanceFromEdge; centerPoint > 0 {
+			gradient = centerPoint / 3.0
+		} else {
+			gradient = centerPoint / 10.0
+		}
+	}
+	var noiseValue float64
+	if gradient >= -2.0 && gradient <= 2.0 {
+		if math.IsNaN(*barrierNoise) {
+			*barrierNoise = a.od.Barrier.Compute(FunctionContext{X: float64(x), Y: float64(y), Z: float64(z)})
+		}
+		noiseValue = *barrierNoise
+	}
+	return 2.0 * (noiseValue + gradient)
+}
+
+// computeFluid decides one aquifer centre's fluid level and type.
+func (a *Aquifer) computeFluid(x, y, z int) FluidStatus {
+	global := a.global(x, y, z)
+	lowestPreliminarySurface := math.MaxInt32
+	topOfCell := y + aquiferYSpacing
+	bottomOfCell := y - aquiferYSpacing
+	surfaceAtCentreIsUnderFluid := false
+	for _, off := range surfaceSamplingOffsets {
+		sampleX := x + off[0]*16
+		sampleZ := z + off[1]*16
+		preliminary := a.od.PreliminarySurfaceLevelAt(sampleX, sampleZ)
+		adjusted := adjustSurfaceLevel(preliminary)
+		start := off[0] == 0 && off[1] == 0
+		// Wholly below the terrain: an ordinary underground aquifer, whose
+		// level the noise decides.
+		if start && bottomOfCell > adjusted {
+			return global
+		}
+		pokesAboveSurface := topOfCell > adjusted
+		if pokesAboveSurface || start {
+			if atSurface := a.global(sampleX, adjusted, sampleZ); atSurface.At(adjusted) != blockAir {
+				if start {
+					surfaceAtCentreIsUnderFluid = true
+				}
+				// Breaking the surface under an ocean: take the ocean's level,
+				// so sea floors do not dry out.
+				if pokesAboveSurface {
+					return atSurface
+				}
+			}
+		}
+		lowestPreliminarySurface = min(lowestPreliminarySurface, preliminary)
+	}
+	level := a.computeSurfaceLevel(x, y, z, global, lowestPreliminarySurface, surfaceAtCentreIsUnderFluid)
+	return FluidStatus{Level: level, Type: a.computeFluidType(x, y, z, global, level)}
+}
+
+func adjustSurfaceLevel(preliminarySurfaceLevel int) int { return preliminarySurfaceLevel + 8 }
+
+// computeSurfaceLevel picks the aquifer's fluid level: the global one when the
+// floodedness noise says "fully flooded", a randomised low one when it says
+// "partially", and nothing at all otherwise — which is what leaves caves dry.
+func (a *Aquifer) computeSurfaceLevel(x, y, z int, global FluidStatus, lowestPreliminarySurface int, surfaceAtCentreIsUnderFluid bool) int {
+	ctx := FunctionContext{X: float64(x), Y: float64(y), Z: float64(z)}
+	var partiallyFloodedness, fullyFloodedness float64
+	if a.isDeepDarkRegion(ctx) {
+		// The deep dark is never flooded.
+		partiallyFloodedness, fullyFloodedness = -1.0, -1.0
+	} else {
+		distanceBelowSurface := lowestPreliminarySurface + 8 - y
+		floodednessFactor := 0.0
+		if surfaceAtCentreIsUnderFluid {
+			floodednessFactor = clampedMap(float64(distanceBelowSurface), 0.0, 64.0, 1.0, 0.0)
+		}
+		floodednessNoise := clamp(a.od.FluidLevelFloodedness.Compute(ctx), -1.0, 1.0)
+		fullyFloodedThreshold := mapRange(floodednessFactor, 1.0, 0.0, -0.3, 0.8)
+		partiallyFloodedThreshold := mapRange(floodednessFactor, 1.0, 0.0, -0.8, 0.4)
+		partiallyFloodedness = floodednessNoise - partiallyFloodedThreshold
+		fullyFloodedness = floodednessNoise - fullyFloodedThreshold
+	}
+	switch {
+	case fullyFloodedness > 0:
+		return global.Level
+	case partiallyFloodedness > 0:
+		return a.computeRandomizedFluidSurfaceLevel(x, y, z, lowestPreliminarySurface)
+	default:
+		return WayBelowMinY
+	}
+}
+
+// computeRandomizedFluidSurfaceLevel puts the water table somewhere in the
+// middle of a 40-block-tall cell, nudged by the spread noise and quantised to
+// three blocks so neighbouring cells share levels often enough to connect.
+func (a *Aquifer) computeRandomizedFluidSurfaceLevel(x, y, z, lowestPreliminarySurface int) int {
+	const cellWidth, cellHeight, maxSpread = 16, 40, 10
+	cellX := floorDivInt(x, cellWidth)
+	cellY := floorDivInt(y, cellHeight)
+	cellZ := floorDivInt(z, cellWidth)
+	middleY := cellY*cellHeight + cellHeight/2
+	spread := a.od.FluidLevelSpread.Compute(FunctionContext{X: float64(cellX), Y: float64(cellY), Z: float64(cellZ)}) * maxSpread
+	return min(lowestPreliminarySurface, middleY+quantizeToMultiple(spread, 3))
+}
+
+// computeFluidType turns deep aquifers into lava lakes.
+func (a *Aquifer) computeFluidType(x, y, z int, global FluidStatus, fluidSurfaceLevel int) uint16 {
+	if fluidSurfaceLevel > -10 || fluidSurfaceLevel == WayBelowMinY || global.Type == blockLava {
+		return global.Type
+	}
+	const cellWidth, cellHeight = 64, 40
+	lavaNoise := a.od.Lava.Compute(FunctionContext{
+		X: float64(floorDivInt(x, cellWidth)),
+		Y: float64(floorDivInt(y, cellHeight)),
+		Z: float64(floorDivInt(z, cellWidth)),
+	})
+	if math.Abs(lavaNoise) > 0.3 {
+		return blockLava
+	}
+	return global.Type
+}
+
+// isDeepDarkRegion is OverworldBiomeBuilder.isDeepDarkRegion.
+func (a *Aquifer) isDeepDarkRegion(ctx FunctionContext) bool {
+	if a.od.Erosion == nil || a.od.Depth == nil {
+		return false
+	}
+	return a.od.Erosion.Compute(ctx) < deepDarkErosionMax && a.od.Depth.Compute(ctx) > deepDarkDepthMin
+}
+
+// ---- grid arithmetic ---------------------------------------------------
+
+func aquiferGridX(blockCoord int) int       { return blockCoord >> 4 }
+func aquiferGridZ(blockCoord int) int       { return blockCoord >> 4 }
+func aquiferGridY(blockCoord int) int       { return floorDivInt(blockCoord, aquiferYSpacing) }
+func fromAquiferGridX(grid, offset int) int { return grid<<4 + offset }
+func fromAquiferGridZ(grid, offset int) int { return grid<<4 + offset }
+func fromAquiferGridY(grid, offset int) int { return grid*aquiferYSpacing + offset }
+
+// floorDivInt is Math.floorDiv: division rounding towards negative infinity.
+func floorDivInt(a, b int) int {
+	q := a / b
+	if a%b != 0 && (a < 0) != (b < 0) {
+		q--
+	}
+	return q
+}
+
+// quantizeToMultiple is Mth.quantize: round down to a multiple of factor.
+func quantizeToMultiple(value float64, factor int) int {
+	return int(math.Floor(value/float64(factor))) * factor
+}
+
+// mapRange is Mth.map: an unclamped linear remap (clampedMap is the clamped one).
+func mapRange(value, from0, to0, from1, to1 float64) float64 {
+	t := (value - from0) / (to0 - from0)
+	return from1 + t*(to1-from1)
+}

internal/worldgen/surface.go

@@ -22,11 +22,17 @@ import (
 type SurfaceContext struct {
 	// X, Y, Z are the block's world coordinates.
 	X, Y, Z int
-	// StoneDepthAbove counts solid blocks at or above Y in this column down to
+	// StoneDepthAbove counts solid blocks from the top of the current stone run
-	// the surface; it is the vanilla "stone_depth" the stone_depth condition
+	// down to and including Y — 1 for the block directly under air or fluid.
-	// compares against (with offset/surface_depth adjustments applied by the
+	// StoneDepthBelow counts the other way, 1 for the block directly above the
-	// test).
+	// cave roof under it. Together they are the vanilla "stone_depth" the
+	// stone_depth condition compares against, floor and ceiling respectively.
 	StoneDepthAbove int
+	StoneDepthBelow int
+	// WaterHeight is one above the lowest fluid block of the run of fluid
+	// directly above Y, or math.MinInt when no fluid sits above Y with no air
+	// in between. It is what the water condition measures against.
+	WaterHeight int
 	// SeaLevel is the world sea level (63 for the overworld).
 	SeaLevel int
 	// BiomeName is the resolved surface biome (e.g. "minecraft:desert").
@@ -210,13 +216,11 @@ func (t yAboveTest) Test(ctx *SurfaceContext) bool {
 	return ctx.Y >= threshold
 }
 
-// stoneDepthTest passes based on the block's depth relative to the surface
+// stoneDepthTest passes when the block is within `offset` of the surface it
-// floor/ceiling. surface_type "floor" counts blocks from the surface downward
+// names: "floor" measures down from the top of the stone run (the ground you
-// and passes when that depth is at or below offset (i.e. near/at the surface);
+// walk on), "ceiling" measures up from its bottom (the roof of whatever cave or
-// "ceiling" passes when the block is the surface cap — the topmost block whose
+// ocean sits underneath). The overworld tree uses ceiling with offset 0 to dress
-// depth-above-surface is 0, i.e. air sits directly on it. This matches vanilla:
+// cave roofs — fourteen times, more than any other stone_depth form.
-// desert's "ceiling → sandstone, else sand" puts sandstone just below the sand
-// cap, not on top.
 type stoneDepthTest struct {
 	surfaceType     string // "floor" or "ceiling"
 	offset          int
@@ -226,16 +230,18 @@ type stoneDepthTest struct {
 
 func (t stoneDepthTest) Test(ctx *SurfaceContext) bool {
 	depth := ctx.StoneDepthAbove
-	if t.addSurfaceDepth {
-		depth += ctx.SurfaceDepth
-	}
 	if t.surfaceType == "ceiling" {
-		// Ceiling: the surface cap. Passes when depth-above-surface equals the
+		depth = ctx.StoneDepthBelow
-		// offset (0 for the topmost block). Used to special-case the block
-		// directly under air.
-		return depth == t.offset
 	}
-	return depth <= t.offset
+	surfaceDepth := 0
+	if t.addSurfaceDepth {
+		surfaceDepth = ctx.SurfaceDepth
+	}
+	// Vanilla widens the band by map(surface_secondary noise, -1..1, 0..range).
+	// That noise is not sampled yet, so the secondary term stays 0; the two
+	// rules that use it also set add_surface_depth, and both currently reduce to
+	// the same single-block band either way.
+	return depth <= 1+t.offset+surfaceDepth
 }
 
 // noiseThresholdTest passes when the named surface noise is within [min,max].