Fix residual surface water blocking rainfall infiltration

src/processes/surface/runoff/direct_surface_runoff.jl

@@ -137,16 +137,17 @@ end
     k_unsat = hydraulic_conductivity(i, j, ground_grid.Nz, grid, fields, soil_hydrology)
     sat_top = saturation_water_ice(i, j, ground_grid.Nz, grid, fields, soil_hydrology)
 
-    if excess_water > zero(NF)
+    # A decaying pool can stay positive at roundoff scale and block fresh rain indefinitely.
+    if excess_water > eps(NF)
         # Case 1: Excess water present at the surface -> precipitation adds to excess water
         # and we set the infiltration rate to the min of hydraulic conductivity and surface_excess_water
         # First, compute rate of excess water removal (surface drainage)
         surface_drainage = compute_surface_drainage(runoff, excess_water)
         # Calculate infiltration
         infil = out.infiltration[i, j, end] = compute_infiltration(runoff, surface_drainage, sat_top, k_unsat)
     else
-        # Case 2: No excess water -> rainfall is routed directly to infiltration
-        surface_drainage = zero(NF)
+        # Case 2: No appreciable excess water -> rainfall is routed directly to infiltration
+        surface_drainage = compute_surface_drainage(runoff, excess_water)
         infil = out.infiltration[i, j, end] = compute_infiltration(runoff, influx, sat_top, k_unsat)
     end
 

test/surface/hydrology/surface_runoff_tests.jl

@@ -57,6 +57,52 @@ end
     @test R ≈ precip + surface_drainage - infil
 end
 
+@testset "fresh rainfall infiltrates after surface water decays" begin
+    grid = ColumnGrid(CPU(), Float32, Float32[-3, -2, -1, -0.5, -0.2, 0])
+    runoff = DirectSurfaceRunoff(Float32)
+    rainfall = 1.0f-7
+    excess_water = eps(Float32) / 2
+    drainage = excess_water / runoff.τ_r
+    fields = (
+        rainfall_ground = fill(rainfall, 1, 1, 1),
+        surface_excess_water = fill(excess_water, 1, 1, 1),
+        saturation_water_ice = fill(0.43f0, 1, 1, 5),
+        hydraulic_conductivity = fill(4.0f-6, 1, 1, 6),
+    )
+    out = (
+        infiltration = zeros(Float32, 1, 1, 1),
+        surface_runoff = zeros(Float32, 1, 1, 1),
+    )
+
+    Terrarium.compute_surface_runoff!(
+        out, 1, 1, grid, fields, runoff,
+        Terrarium.NoCanopyInterception(Float32),
+        SoilHydrology(Float32, NoFlow())
+    )
+
+    @test out.infiltration[1, 1, 1] ≈ rainfall
+    @test out.surface_runoff[1, 1, 1] ≈ drainage rtol = 1.0f-3
+
+    # The seasonal replay reached this subnormal pool value; its drainage rounds to zero.
+    fields.surface_excess_water[1, 1, 1] = 2.388f-42
+    Terrarium.compute_surface_runoff!(
+        out, 1, 1, grid, fields, runoff,
+        Terrarium.NoCanopyInterception(Float32),
+        SoilHydrology(Float32, NoFlow())
+    )
+    @test out.infiltration[1, 1, 1] ≈ rainfall
+    @test iszero(out.surface_runoff[1, 1, 1])
+
+    fields.surface_excess_water[1, 1, 1] = 1.0f-4
+    Terrarium.compute_surface_runoff!(
+        out, 1, 1, grid, fields, runoff,
+        Terrarium.NoCanopyInterception(Float32),
+        SoilHydrology(Float32, NoFlow())
+    )
+    @test out.infiltration[1, 1, 1] ≈ 1.0f-4 / runoff.τ_r
+    @test out.surface_runoff[1, 1, 1] ≈ rainfall
+end
+
 @testset "surface_excess_water tendency" begin
     grid = ColumnGrid(UniformSpacing(Δz = 0.1, N = 10))