From 0393575f0846ce395f63b7918df9faa9728d8b0a Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 13 Aug 2025 17:21:31 +0200
Subject: [PATCH 01/24] allow for infiltration of overlandflow

---
 Wflow/src/sbm.jl               |   2 +-
 Wflow/src/sbm_model.jl         |   6 +-
 Wflow/src/soil/soil.jl         | 112 ++++++++++++++++++++++++++++++++-
 Wflow/src/soil/soil_process.jl |  67 ++++++++++++++++++++
 4 files changed, 180 insertions(+), 7 deletions(-)

diff --git a/Wflow/src/sbm.jl b/Wflow/src/sbm.jl
index fe5291c7f..e50c29067 100644
--- a/Wflow/src/sbm.jl
+++ b/Wflow/src/sbm.jl
@@ -116,7 +116,7 @@ function update!(
         (; interception, runoff, demand, allocation),
     )
 
-    update!(soil, atmospheric_forcing, (; snow, runoff, demand), config, dt)
+    update!(soil, domain, atmospheric_forcing, (; snow, runoff, demand), config, dt)
     @. soil.variables.actevap += interception.variables.interception_rate
     return nothing
 end
diff --git a/Wflow/src/sbm_model.jl b/Wflow/src/sbm_model.jl
index 2c34ba80a..99d2ba1f0 100644
--- a/Wflow/src/sbm_model.jl
+++ b/Wflow/src/sbm_model.jl
@@ -107,12 +107,12 @@ Update SBM model after subsurface flow for a single timestep. This function is a
 accessible through BMI, to couple the SBM model to an external groundwater model.
 """
 function update_after_subsurfaceflow!(model::AbstractModel{<:SbmModel})
-    (; routing, land) = model
+    (; routing, land, domain, config) = model
     (; soil, runoff, demand) = land
-    (; subsurface_flow) = routing
+    (; subsurface_flow, overland_flow) = routing
 
     # update SBM soil model (runoff, ustorelayerdepth and satwaterdepth)
-    update!(soil, (; runoff, demand, subsurface_flow))
+    update!(soil, (; runoff, demand, subsurface_flow, overland_flow), domain, config)
 
     surface_routing!(model)
 
diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index cc53c396c..66bdb1ed4 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -40,6 +40,12 @@ abstract type AbstractSoilModel end
     infiltsoilpath::Vector{Float64}
     # Infiltration excess water [mm Δt⁻¹]
     infiltexcess::Vector{Float64}
+    # Infiltration from surface water [mm Δt⁻¹]
+    infilt_surfacewater::Vector{Float64}
+    # Contribution from surface water [mm Δt⁻¹]
+    infilt_available_surfacewater::Vector{Float64}
+    # Total water available for infiltration [mm Δt⁻¹]
+    infilt_total_available::Vector{Float64}
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{Float64}
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -184,6 +190,9 @@ function SbmSoilVariables(n::Int, parameters::SbmSoilParameters)
         actinfiltpath = fill(MISSING_VALUE, n),
         infiltsoilpath = fill(MISSING_VALUE, n),
         infiltexcess = fill(MISSING_VALUE, n),
+        infilt_surfacewater = fill(0.0, n),
+        infilt_available_surfacewater = fill(0.0, n),
+        infilt_total_available = fill(0.0, n),
         excesswater = fill(MISSING_VALUE, n),
         exfiltsatwater = fill(MISSING_VALUE, n),
         exfiltustore = fill(MISSING_VALUE, n),
@@ -680,6 +689,83 @@ function infiltration_reduction_factor!(
     return nothing
 end
 
+function update_available_for_infiltration!(
+    model::SbmSoilModel,
+    domain::Domain,
+    runoff::AbstractRunoffModel,
+    do_surface_water_infiltration::Bool,
+)
+    v = model.variables
+    (; water_flux_surface) = model.boundary_conditions
+    (; waterdepth_land) = runoff.boundary_conditions
+    (; river_fraction) = domain.land.parameters
+
+    n = length(v.infilt_total_available)
+    threaded_foreach(1:n; basesize = 1000) do i
+        v.infilt_available_surfacewater[i] = 0.0
+        if do_surface_water_infiltration
+            v.infilt_available_surfacewater[i] =
+                waterdepth_land[i] * (1.0 - river_fraction[i]) * 0.95
+            water_flux_surface[i] += v.infilt_available_surfacewater[i]
+        end
+        v.infilt_total_available[i] = water_flux_surface[i]
+    end
+
+    return nothing
+end
+
+function correct_infiltration!(model::SbmSoilModel)
+    v = model.variables
+    (; water_flux_surface) = model.boundary_conditions
+
+    n = length(v.actinfilt)
+    threaded_foreach(1:n; basesize = 1000) do i
+        v.infilt_surfacewater[i],
+        v.actinfilt[i],
+        v.infiltexcess[i],
+        v.excesswater[i],
+        water_flux_surface[i] = correct_infiltration(
+            v.infilt_total_available[i],
+            v.infilt_available_surfacewater[i],
+            water_flux_surface[i],
+            v.actinfilt[i],
+            v.infiltexcess[i],
+        )
+    end
+end
+
+function correct_overland_flow_level!(
+    model::SbmSoilModel,
+    # overland_flow::Union{KinWaveOverlandFlow, LocalInertialOverlandFlow},
+    overland_flow,
+    domain::Domain,
+    config::Config,
+)
+    v = model.variables
+
+    do_surface_water_infiltration =
+        get(config.model, "reinfiltration_surfacewater", false)::Bool
+
+    if do_surface_water_infiltration
+        (; surface_flow_width) = domain.land.parameters
+        n = length(surface_flow_width)
+        threaded_foreach(1:n; basesize = 1000) do i
+            q, h = correct_overland_flow_level(
+                overland_flow.variables.h[i],
+                v.infilt_surfacewater[i],
+                domain.land.parameters.river_fraction[i],
+                surface_flow_width[i],
+                overland_flow.parameters.alpha[i],
+                overland_flow.parameters.beta,
+            )
+            if !isnothing(q)
+                overland_flow.variables.flow.q[i] = q
+                overland_flow.variables.h[i] = h
+            end
+        end
+    end
+end
+
 """
     infiltration!(model::SbmSoilModel)
 
@@ -1054,6 +1140,7 @@ transpiration, capillary flux and leakage) for a single timestep.
 """
 function update!(
     model::SbmSoilModel,
+    domain::Domain,
     atmospheric_forcing::AtmosphericForcing,
     external_models::NamedTuple,
     config::Config,
@@ -1062,6 +1149,8 @@ function update!(
     soil_infiltration_reduction =
         get(config.model, "soil_infiltration_reduction__flag", false)::Bool
     modelsnow = get(config.model, "snow__flag", false)::Bool
+    do_surface_water_infiltration =
+        get(config.model, "reinfiltration_surfacewater", false)::Bool
 
     (; snow, runoff, demand) = external_models
     (; temperature) = atmospheric_forcing
@@ -1077,6 +1166,10 @@ function update!(
     # infiltration
     soil_temperature!(model, snow, temperature)
     infiltration_reduction_factor!(model; modelsnow, soil_infiltration_reduction)
+
+    # update available for infiltration in case surface water infiltration is enabled
+    update_available_for_infiltration!(model, domain, runoff, do_surface_water_infiltration)
+
     infiltration!(model)
     # unsaturated zone flow
     unsaturated_zone_flow!(model)
@@ -1085,7 +1178,11 @@ function update!(
     transpiration!(model, dt)
     # actual infiltration and excess water
     actual_infiltration!(model)
-    @. v.excesswater = water_flux_surface - v.actinfilt - v.infiltexcess
+
+    # Correct fluxes in case of reinfiltration, also to ensure correct soil and path
+    # infiltration, and excesswater
+    correct_infiltration!(model)
+
     actual_infiltration_soil_path!(model)
     @. v.excesswatersoil =
         max(water_flux_surface * (1.0 - p.pathfrac) - v.actinfiltsoil, 0.0)
@@ -1118,8 +1215,13 @@ the unsaturated store `exfiltustore`, land `runoff` and `net_runoff`, the satura
 `exfiltsatwater` are updated. Addionally, volumetric water content per soil layer and for
 the root zone are updated.
 """
-function update!(model::SbmSoilModel, external_models::NamedTuple)
-    (; runoff, demand, subsurface_flow) = external_models
+function update!(
+    model::SbmSoilModel,
+    external_models::NamedTuple,
+    domain::Domain,
+    config::Config,
+)
+    (; runoff, demand, subsurface_flow, overland_flow) = external_models
     (; runoff_land, ae_openw_l) = runoff.variables
     p = model.parameters
     v = model.variables
@@ -1223,6 +1325,10 @@ function update!(model::SbmSoilModel, external_models::NamedTuple)
     # and the h_max parameter of a paddy field)
     update_runoff!(demand.paddy, v.runoff)
     @. v.net_runoff = v.runoff - ae_openw_l
+
+    # correct overland flow water levels in case of reinfiltration
+    correct_overland_flow_level!(model, overland_flow, domain, config)
+
     return nothing
 end
 
diff --git a/Wflow/src/soil/soil_process.jl b/Wflow/src/soil/soil_process.jl
index f1314bbff..cf43c0786 100644
--- a/Wflow/src/soil/soil_process.jl
+++ b/Wflow/src/soil/soil_process.jl
@@ -292,4 +292,71 @@ function actual_infiltration_soil_path(
     end
 
     return actinfiltsoil, actinfiltpath
+end
+
+function correct_infiltration(
+    infilt_total_available,
+    infilt_available_surfacewater,
+    water_flux_surface,
+    actinfilt,
+    infiltexcess,
+)
+    # Determine ratio of water that has infiltrated
+    infilt_ratio = infilt_total_available == 0.0 ? 0.0 : actinfilt / infilt_total_available
+    # Use this ratio to determine the contribution from overland flow
+    infilt_surfacewater = max(0.0, infilt_available_surfacewater * infilt_ratio)
+    # Determine the correction factor to apply to the relevant fluxes
+    correction_surfacewater =
+        infilt_total_available == 0.0 ? 1.0 :
+        1.0 - (infilt_available_surfacewater / infilt_total_available)
+
+    # Correct fluxes
+    actinfilt *= correction_surfacewater
+    infiltexcess *= correction_surfacewater
+
+    # subtract contribution from overland flow to ensure correct fluxes
+    water_flux_surface -= infilt_available_surfacewater
+    excesswater = water_flux_surface - actinfilt - infiltexcess
+
+    return infilt_surfacewater, actinfilt, infiltexcess, excesswater, water_flux_surface
+end
+
+function correct_overland_flow_level(
+    overlandflow_depth,
+    infilt_surfacewater,
+    river_fraction,
+    surface_flow_width,
+    alpha,
+    beta,
+)
+    if infilt_surfacewater > 0.0
+        # Get original h_land in mm
+        original_h_land = overlandflow_depth * 1000.0
+
+        # Correct values for river fraction to ensure correct water accounting
+        infiltrated_surfacewater = (infilt_surfacewater / (1.0 - river_fraction))
+        # Calculate new h_land in m
+        h = (original_h_land - infiltrated_surfacewater) / 1000.0
+
+        diff = overlandflow_depth - h
+        if diff < 0.0 && infilt_surfacewater > 0.0
+            println(
+                "Corrected h_land: ",
+                overlandflow_depth - h,
+                " ",
+                infiltrated_surfacewater,
+            )
+        end
+
+        q = ifelse(
+            surface_flow_width > 0.0 && h > 0.0 && alpha > 0.0 && beta != 0.0,
+            # Compute cross-sectional area from h
+            pow.(max.((h * surface_flow_width) / alpha, 1e-10), 1.0 / beta),
+            0.0,  # Set q to 0.0 if conditions are not met
+        )
+    else
+        q = nothing
+        h = overlandflow_depth
+    end
+    return q, h
 end
\ No newline at end of file

From 4cc452ceda4f9c0956c81c665a3ee636b4b2c173 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 2 Sep 2025 14:46:40 +0200
Subject: [PATCH 02/24] fix gwf model

---
 Wflow/src/sbm_gwf_model.jl     |  9 +++++----
 Wflow/src/soil/soil_process.jl | 10 ----------
 2 files changed, 5 insertions(+), 14 deletions(-)

diff --git a/Wflow/src/sbm_gwf_model.jl b/Wflow/src/sbm_gwf_model.jl
index e6752d4ac..0595bf242 100644
--- a/Wflow/src/sbm_gwf_model.jl
+++ b/Wflow/src/sbm_gwf_model.jl
@@ -64,9 +64,10 @@ end
 function update!(model::AbstractModel{<:SbmGwfModel})
     (; routing, land, domain, clock, config) = model
     (; soil, runoff, demand) = land
+    (; subsurface_flow, overland_flow) = routing
 
     do_water_demand = haskey(config.model, "water_demand")
-    (; aquifer, boundaries) = routing.subsurface_flow
+    (; aquifer, boundaries) = subsurface_flow
     dt = tosecond(clock.dt)
 
     update!(land, routing, domain, config, dt)
@@ -89,7 +90,7 @@ function update!(model::AbstractModel{<:SbmGwfModel})
         )
     end
 
-    Q = zeros(routing.subsurface_flow.connectivity.ncell)
+    Q = zeros(subsurface_flow.connectivity.ncell)
     # exchange of recharge between SBM soil model and groundwater flow domain
     # recharge rate groundwater is required in units [m d⁻¹]
     @. boundaries.recharge.variables.rate =
@@ -99,10 +100,10 @@ function update!(model::AbstractModel{<:SbmGwfModel})
             land.allocation.variables.act_groundwater_abst / 1000.0 * (1.0 / dt_sbm)
     end
     # update groundwater domain
-    update!(routing.subsurface_flow, Q, dt_sbm, conductivity_profile)
+    update!(subsurface_flow, Q, dt_sbm, conductivity_profile)
 
     # update SBM soil model (runoff, ustorelayerdepth and satwaterdepth)
-    update!(soil, (; runoff, demand, subsurface_flow = routing.subsurface_flow))
+    update!(soil, (; runoff, demand, subsurface_flow, overland_flow), domain, config)
 
     surface_routing!(model)
 
diff --git a/Wflow/src/soil/soil_process.jl b/Wflow/src/soil/soil_process.jl
index cf43c0786..9eadf4019 100644
--- a/Wflow/src/soil/soil_process.jl
+++ b/Wflow/src/soil/soil_process.jl
@@ -338,16 +338,6 @@ function correct_overland_flow_level(
         # Calculate new h_land in m
         h = (original_h_land - infiltrated_surfacewater) / 1000.0
 
-        diff = overlandflow_depth - h
-        if diff < 0.0 && infilt_surfacewater > 0.0
-            println(
-                "Corrected h_land: ",
-                overlandflow_depth - h,
-                " ",
-                infiltrated_surfacewater,
-            )
-        end
-
         q = ifelse(
             surface_flow_width > 0.0 && h > 0.0 && alpha > 0.0 && beta != 0.0,
             # Compute cross-sectional area from h

From c2108d9fadfb44a7d7e393a4b6caff4ea58a7a0c Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 16 Sep 2025 15:41:25 +0200
Subject: [PATCH 03/24] overlandflow correction improvement

---
 Wflow/src/routing/surface_kinwave.jl        | 37 +++++++++++++++++-
 Wflow/src/routing/surface_local_inertial.jl | 42 ++++++++++++++++++++-
 Wflow/src/soil/soil.jl                      | 32 ----------------
 3 files changed, 77 insertions(+), 34 deletions(-)

diff --git a/Wflow/src/routing/surface_kinwave.jl b/Wflow/src/routing/surface_kinwave.jl
index 387578d29..22453ccaa 100644
--- a/Wflow/src/routing/surface_kinwave.jl
+++ b/Wflow/src/routing/surface_kinwave.jl
@@ -525,7 +525,7 @@ end
 
 """
 Compute a stable timestep size for the kinematice wave method for a river or overland flow
-model using a nonlinear scheme (Chow et al., 1988). 
+model using a nonlinear scheme (Chow et al., 1988).
 
 A stable time step is computed for each vector element based on the Courant timestep size
 criterion. A quantile of the vector is computed based on probability `p` to remove potential
@@ -650,3 +650,38 @@ get_inflow_reservoir(::KinWaveRiverFlow, model::LateralSSF) =
 # Exclude subsurface flow from `GroundwaterFlow`.
 get_inflow_reservoir(::AbstractRiverFlowModel, model::GroundwaterFlow) =
     zeros(model.connectivity.ncell)
+
+"""
+Update overland flow water level and discharge for KinWaveOverlandFlow model based on
+surface water infiltration.
+"""
+function correct_overland_flow_level!(
+    model::SbmSoilModel,
+    overland_flow::KinWaveOverlandFlow,
+    domain::Domain,
+    config::Config,
+)
+    v = model.variables
+
+    do_surface_water_infiltration =
+        get(config.model, "reinfiltration_surfacewater", false)::Bool
+
+    if do_surface_water_infiltration
+        (; surface_flow_width) = domain.land.parameters
+        n = length(surface_flow_width)
+        threaded_foreach(1:n; basesize = 1000) do i
+            q, h = correct_overland_flow_level(
+                overland_flow.variables.h[i],
+                v.infilt_surfacewater[i],
+                domain.land.parameters.river_fraction[i],
+                surface_flow_width[i],
+                overland_flow.parameters.alpha[i],
+                overland_flow.parameters.beta,
+            )
+            if !isnothing(q)
+                overland_flow.variables.flow.q[i] = q
+                overland_flow.variables.h[i] = h
+            end
+        end
+    end
+end
diff --git a/Wflow/src/routing/surface_local_inertial.jl b/Wflow/src/routing/surface_local_inertial.jl
index 4df7279f8..75a90be05 100644
--- a/Wflow/src/routing/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface_local_inertial.jl
@@ -617,11 +617,16 @@ end
 @with_kw struct LocalInertialOverlandFlowBC
     runoff::Vector{Float64}           # runoff from hydrological model [m³ s⁻¹]
     inflow_reservoir::Vector{Float64} # inflow to reservoir from hydrological model [m³ s⁻¹]
+    infiltration_volume::Vector{Float64} # amount of infiltration from surface water [m³]
 end
 
 "Struct to store shallow water overland flow model boundary conditions"
 function LocalInertialOverlandFlowBC(n::Int)
-    bc = LocalInertialOverlandFlowBC(; runoff = zeros(n), inflow_reservoir = zeros(n))
+    bc = LocalInertialOverlandFlowBC(;
+        runoff = zeros(n),
+        inflow_reservoir = zeros(n),
+        infiltration_volume = zeros(n),
+    )
     return bc
 end
 
@@ -962,6 +967,11 @@ function local_inertial_update_water_depth!(
                     land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] + land_bc.runoff[i] -
                     river_bc.abstraction[inds_river[i]]
                 ) * dt
+            # Apply surface water infiltration correction for river cells
+            if land_bc.infiltration_volume[i] > 0.0
+                land_v.storage[i] =
+                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
+            end
             if land_v.storage[i] < 0.0
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
@@ -1012,6 +1022,11 @@ function local_inertial_update_water_depth!(
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
             end
+            # Apply surface water infiltration correction if available
+            if land_bc.infiltration_volume[i] > 0.0
+                land_v.storage[i] =
+                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
+            end
             land_v.h[i] = land_v.storage[i] / (x_length[i] * y_length[i])
         end
         # average variables (here accumulated for model timestep Δt)
@@ -1300,3 +1315,28 @@ function FloodPlain(
     floodplain = FloodPlain(; parameters, variables)
     return floodplain
 end
+
+"""
+Update overland flow water level and storage for LocalInertialOverlandFlow model based on
+surface water infiltration. For local inertial flow, the infiltration is applied by updating
+the boundary conditions which are then used in the water depth update function.
+"""
+function correct_overland_flow_level!(
+    model::SbmSoilModel,
+    overland_flow::LocalInertialOverlandFlow,
+    domain::Domain,
+    config::Config,
+)
+    (; infilt_surfacewater) = model.variables
+    (; area) = domain.land.parameters
+
+    do_surface_water_infiltration =
+        get(config.model, "reinfiltration_surfacewater", false)::Bool
+
+    if do_surface_water_infiltration
+        # Update the boundary condition for surface water infiltration
+        # This will be used in local_inertial_update_water_depth!
+        overland_flow.boundary_conditions.infiltration_volume .=
+            infilt_surfacewater .* area .* 0.001
+    end
+end
diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 66bdb1ed4..b31929d9a 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -734,38 +734,6 @@ function correct_infiltration!(model::SbmSoilModel)
     end
 end
 
-function correct_overland_flow_level!(
-    model::SbmSoilModel,
-    # overland_flow::Union{KinWaveOverlandFlow, LocalInertialOverlandFlow},
-    overland_flow,
-    domain::Domain,
-    config::Config,
-)
-    v = model.variables
-
-    do_surface_water_infiltration =
-        get(config.model, "reinfiltration_surfacewater", false)::Bool
-
-    if do_surface_water_infiltration
-        (; surface_flow_width) = domain.land.parameters
-        n = length(surface_flow_width)
-        threaded_foreach(1:n; basesize = 1000) do i
-            q, h = correct_overland_flow_level(
-                overland_flow.variables.h[i],
-                v.infilt_surfacewater[i],
-                domain.land.parameters.river_fraction[i],
-                surface_flow_width[i],
-                overland_flow.parameters.alpha[i],
-                overland_flow.parameters.beta,
-            )
-            if !isnothing(q)
-                overland_flow.variables.flow.q[i] = q
-                overland_flow.variables.h[i] = h
-            end
-        end
-    end
-end
-
 """
     infiltration!(model::SbmSoilModel)
 

From e141d87584ebbaeb0459db6e8a4e4222b0b3fa31 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:07:06 +0100
Subject: [PATCH 04/24] Update config_structure.jl

---
 Wflow/src/config_structure.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Wflow/src/config_structure.jl b/Wflow/src/config_structure.jl
index b18cf4bde..9060ec8c7 100644
--- a/Wflow/src/config_structure.jl
+++ b/Wflow/src/config_structure.jl
@@ -70,6 +70,7 @@ end
     snow_gravitational_transport__flag::Bool = false
     glacier__flag::Bool = false
     soil_infiltration_reduction__flag::Bool = false
+    reinfiltration_surfacewater__flag::Bool = false
     soil_layer__thickness::Vector{Int} = [100, 300, 800]
     saturated_hydraulic_conductivity_profile::VerticalConductivityProfile.T =
         VerticalConductivityProfile.exponential

From 60585671f6a1bc192f808527cfd5b02e4bc72104 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:22:35 +0100
Subject: [PATCH 05/24] fix flags

---
 Wflow/src/routing/surface_kinwave.jl        | 5 +----
 Wflow/src/routing/surface_local_inertial.jl | 6 +-----
 2 files changed, 2 insertions(+), 9 deletions(-)

diff --git a/Wflow/src/routing/surface_kinwave.jl b/Wflow/src/routing/surface_kinwave.jl
index d1d54b61f..eecdbae52 100644
--- a/Wflow/src/routing/surface_kinwave.jl
+++ b/Wflow/src/routing/surface_kinwave.jl
@@ -652,10 +652,7 @@ function correct_overland_flow_level!(
 )
     v = model.variables
 
-    do_surface_water_infiltration =
-        get(config.model, "reinfiltration_surfacewater", false)::Bool
-
-    if do_surface_water_infiltration
+    if config.model.reinfiltration_surfacewater__flag
         (; surface_flow_width) = domain.land.parameters
         n = length(surface_flow_width)
         threaded_foreach(1:n; basesize = 1000) do i
diff --git a/Wflow/src/routing/surface_local_inertial.jl b/Wflow/src/routing/surface_local_inertial.jl
index 8fcdab27e..b25f95a3a 100644
--- a/Wflow/src/routing/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface_local_inertial.jl
@@ -646,7 +646,6 @@ end
     infiltration_volume::Vector{Float64} = zeros(n) # amount of infiltration from surface water [m³]
 end
 
-
 "Local inertial overland flow model using the local inertial method"
 @with_kw struct LocalInertialOverlandFlow <: AbstractOverlandFlowModel
     timestepping::TimeStepping
@@ -1348,10 +1347,7 @@ function correct_overland_flow_level!(
     (; infilt_surfacewater) = model.variables
     (; area) = domain.land.parameters
 
-    do_surface_water_infiltration =
-        get(config.model, "reinfiltration_surfacewater", false)::Bool
-
-    if do_surface_water_infiltration
+    if config.model.reinfiltration_surfacewater__flag
         # Update the boundary condition for surface water infiltration
         # This will be used in local_inertial_update_water_depth!
         overland_flow.boundary_conditions.infiltration_volume .=

From 78ff48d8ce0c04b43388397cddfb8ad54fe2660c Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:27:51 +0100
Subject: [PATCH 06/24] precommit fix

---
 Wflow/src/soil/soil_process.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Wflow/src/soil/soil_process.jl b/Wflow/src/soil/soil_process.jl
index 05572217e..da4cf3f05 100644
--- a/Wflow/src/soil/soil_process.jl
+++ b/Wflow/src/soil/soil_process.jl
@@ -323,4 +323,4 @@ function correct_overland_flow_level(
         h = overlandflow_depth
     end
     return q, h
-end
\ No newline at end of file
+end

From 9dcf097510a8e7cb46c88251789327f4f8a90eb7 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 13 Feb 2026 14:32:11 +0100
Subject: [PATCH 07/24] remove todo statement

---
 Wflow/src/soil/soil.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index f87078544..52b2f2406 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -1177,7 +1177,6 @@ function update!(
     config::Config,
     dt::Float64,
 )
-
     (; snow, runoff, demand) = external_models
     (; temperature) = atmospheric_forcing
     (; water_flux_surface) = model.boundary_conditions
@@ -1195,12 +1194,11 @@ function update!(
     )
 
     # update available for infiltration in case surface water infiltration is enabled
-    #TODO fix flag
     update_available_for_infiltration!(
         model,
         domain,
         runoff,
-        config.model.reinfiltration_surfacewater__flag
+        config.model.reinfiltration_surfacewater__flag,
     )
 
     infiltration!(model)

From 2aeb128e3bc901586fc84af25d7c2347f45bb2ed Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 6 Mar 2026 16:12:23 +0100
Subject: [PATCH 08/24] fix merge error

---
 Wflow/src/soil/soil.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 5d8c5f53e..6ffba948f 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -1393,7 +1393,7 @@ function update_soil_water_storage!(
     @. v.net_runoff = v.runoff - ae_openw_l
 
     # correct overland flow water levels in case of reinfiltration
-    correct_overland_flow_level!(model, overland_flow, domain, config)
+    correct_overland_flow_level!(soil_model, overland_flow, domain, config)
 
     return nothing
 end

From dfd4afb13258eae0b369ba9b69acda7e73cc6586 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 6 Mar 2026 16:35:25 +0100
Subject: [PATCH 09/24] add unit test

---
 Wflow/test/land_process.jl | 84 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 84 insertions(+)

diff --git a/Wflow/test/land_process.jl b/Wflow/test/land_process.jl
index dd39a5e4b..af1149298 100644
--- a/Wflow/test/land_process.jl
+++ b/Wflow/test/land_process.jl
@@ -445,3 +445,87 @@ end
     @test actinfiltsoil == 0.0
     @test actinfiltpath == 0.0
 end
+
+@testitem "unit: correct_infiltration" begin
+    potential_infiltration = 10.0
+    infilt_available_surfacewater = 2.0
+    water_flux_surface = 10.0
+    actual_infiltration = 6.0
+    infiltexcess = 1.0
+
+    infilt_surfacewater,
+    actual_infiltration,
+    infiltexcess,
+    excesswater,
+    water_flux_surface = Wflow.correct_infiltration(
+        potential_infiltration,
+        infilt_available_surfacewater,
+        water_flux_surface,
+        actual_infiltration,
+        infiltexcess,
+    )
+
+    @test infilt_surfacewater ≈ 1.2
+    @test actual_infiltration ≈ 4.8
+    @test infiltexcess ≈ 0.8
+    @test water_flux_surface ≈ 8.0
+    @test excesswater ≈ 2.4
+
+    potential_infiltration = 0.0
+    infilt_available_surfacewater = 1.0
+    water_flux_surface = 3.0
+    actual_infiltration = 0.0
+    infiltexcess = 0.0
+
+    infilt_surfacewater,
+    actual_infiltration,
+    infiltexcess,
+    excesswater,
+    water_flux_surface = Wflow.correct_infiltration(
+        potential_infiltration,
+        infilt_available_surfacewater,
+        water_flux_surface,
+        actual_infiltration,
+        infiltexcess,
+    )
+
+    @test infilt_surfacewater == 0.0
+    @test actual_infiltration == 0.0
+    @test infiltexcess == 0.0
+    @test water_flux_surface ≈ 2.0
+    @test excesswater ≈ 2.0
+end
+
+@testitem "unit: correct_overland_flow_level" begin
+    overlandflow_depth = 0.02
+    infilt_surfacewater = 5.0
+    river_fraction = 0.2
+    surface_flow_width = 10.0
+    alpha = 2.0
+    beta = 0.5
+
+    q, h = Wflow.correct_overland_flow_level(
+        overlandflow_depth,
+        infilt_surfacewater,
+        river_fraction,
+        surface_flow_width,
+        alpha,
+        beta,
+    )
+
+    @test h ≈ 0.01375
+    @test q ≈ 0.0047265625
+
+    infilt_surfacewater = 0.0
+    q, h = Wflow.correct_overland_flow_level(
+        overlandflow_depth,
+        infilt_surfacewater,
+        river_fraction,
+        surface_flow_width,
+        alpha,
+        beta,
+    )
+
+    @test q === nothing
+    @test h == overlandflow_depth
+end

From 84cdf2228ba3bd485d747390da5231240808f5b0 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Mar 2026 11:26:16 +0100
Subject: [PATCH 10/24] fix merge confict error

---
 Wflow/src/sbm_model.jl | 31 +------------------------------
 1 file changed, 1 insertion(+), 30 deletions(-)

diff --git a/Wflow/src/sbm_model.jl b/Wflow/src/sbm_model.jl
index 65cee0030..6685e4cc2 100644
--- a/Wflow/src/sbm_model.jl
+++ b/Wflow/src/sbm_model.jl
@@ -61,6 +61,7 @@ function update_model!(model::AbstractModel{<:SbmModel})
     (; routing, land, domain, clock, config) = model
     (; soil, runoff, demand) = land
     (; kv_profile) = land.soil.parameters
+    (; subsurface_flow, overland_flow) = routing
     (; boundary_conditions) = routing.subsurface_flow
     dt = tosecond(clock.dt)
 
@@ -86,36 +87,6 @@ function update_model!(model::AbstractModel{<:SbmModel})
         clock.dt / BASETIMESTEP,
     )
 
-    update_after_subsurface_flow!(model)
-    update_total_water_storage!(model)
-    return nothing
-end
-
-"""
-    update_until_recharge!model::AbstractModel{<:SbmModel})
-
-Update SBM model until recharge for a single timestep. This function is also accessible
-through BMI, to couple the SBM model to an external groundwater model.
-"""
-function update_until_recharge!(model::AbstractModel{<:SbmModel})
-    (; routing, land, domain, clock, config) = model
-    dt = tosecond(clock.dt)
-    update_land_hydrology_model!(land, routing, domain, config, dt)
-    return nothing
-end
-
-"""
-    update_after_subsurfaceflow!(model::AbstractModel{<:SbmModel})
-
-Update SBM model after subsurface flow for a single timestep. This function is also
-accessible through BMI, to couple the SBM model to an external groundwater model.
-"""
-function update_after_subsurface_flow!(model::AbstractModel{<:SbmModel})
-    (; routing, land, domain, config) = model
-
-    (; soil, runoff, demand) = land
-    (; subsurface_flow, overland_flow) = routing
-
     # update SBM soil model (runoff, ustorelayerdepth and satwaterdepth)
     update_soil_water_storage!(
         soil,

From 49b9f633fc4bf38d29c71c541cf78f0a6afd80ad Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Mar 2026 11:33:03 +0100
Subject: [PATCH 11/24] another merge fix

---
 Wflow/src/routing/surface/surface_kinwave.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Wflow/src/routing/surface/surface_kinwave.jl b/Wflow/src/routing/surface/surface_kinwave.jl
index 97cf06a04..a8afaff0a 100644
--- a/Wflow/src/routing/surface/surface_kinwave.jl
+++ b/Wflow/src/routing/surface/surface_kinwave.jl
@@ -656,7 +656,7 @@ surface water infiltration.
 """
 function correct_overland_flow_level!(
     model::SbmSoilModel,
-    overland_flow::KinWaveOverlandFlow,
+    overland_flow::KinWaveOverlandFlowModel,
     domain::Domain,
     config::Config,
 )

From 23567ef2f1498011df389679e48bbd638f1b980f Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Mar 2026 11:35:25 +0100
Subject: [PATCH 12/24] Update surface_local_inertial.jl

---
 Wflow/src/routing/surface/surface_local_inertial.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Wflow/src/routing/surface/surface_local_inertial.jl b/Wflow/src/routing/surface/surface_local_inertial.jl
index a7eacfedb..18209fe3e 100644
--- a/Wflow/src/routing/surface/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface/surface_local_inertial.jl
@@ -1372,7 +1372,7 @@ the boundary conditions which are then used in the water depth update function.
 """
 function correct_overland_flow_level!(
     model::SbmSoilModel,
-    overland_flow::LocalInertialOverlandFlow,
+    overland_flow::LocalInertialOverlandFlowModel,
     domain::Domain,
     config::Config,
 )

From 92cfb0015e269af11ef2577d425fa6f844b27b4f Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Mar 2026 09:23:46 +0100
Subject: [PATCH 13/24] fix test

---
 Wflow/test/soil.jl | 42 ++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 40 insertions(+), 2 deletions(-)

diff --git a/Wflow/test/soil.jl b/Wflow/test/soil.jl
index a22ba811b..44558ba6f 100644
--- a/Wflow/test/soil.jl
+++ b/Wflow/test/soil.jl
@@ -274,9 +274,47 @@ end
     )
 
     subsurface_flow = (; variables = (; exfiltwater = [0.0]))
-    external_models = (; runoff, demand, subsurface_flow)
+    overland_flow = Wflow.KinWaveOverlandFlowModel(;
+        timestepping = Wflow.TimeStepping(),
+        boundary_conditions = Wflow.LandFlowBC(; n),
+        parameters = Wflow.ManningFlowParameters(;
+            n,
+            beta = 0.6,
+            slope = [0.01],
+            mannings_n = [0.072],
+            alpha_pow = (2.0 / 3.0) * 0.6,
+        ),
+        variables = Wflow.OverLandFlowVariables(; n),
+    )
+    external_models = (; runoff, demand, subsurface_flow, overland_flow)
+
+    domain = Wflow.Domain(;
+        land = Wflow.DomainLand(;
+            parameters = Wflow.LandParameters(;
+                surface_flow_width = [700.0],
+                river_fraction = [0.1],
+            ),
+        ),
+    )
+    config = Wflow.Config(;
+        model = Wflow.ModelSection(;
+            type = Wflow.ModelType.sbm,
+            reinfiltration_surfacewater__flag = false,
+        ),
+        input = Wflow.InputSection(;
+            path_forcing = "",
+            path_static = "",
+            basin__local_drain_direction = "",
+            river_location__mask = "",
+            subbasin_location__count = "",
+            forcing = Wflow.InputEntries(),
+            static = Wflow.InputEntries(),
+            location_maps = Wflow.PropertyDict(Dict{String, Any}()),
+        ),
+        path = "",
+    )
 
-    Wflow.update_soil_water_storage!(soil_model, external_models)
+    Wflow.update_soil_water_storage!(soil_model, external_models, domain, config)
 
     @test demand.paddy.variables.h[1] ≈ 0.047112555773
     @test soil_model.variables.runoff[1] ≈ 0.0

From abecbff6811d605f86fa77fbbc692b3a6755a345 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 2 Apr 2026 09:54:48 +0200
Subject: [PATCH 14/24] add docs

---
 docs/model_docs/land/landhydrology_sbm.qmd | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index f50b362d2..1f16b2228 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -960,6 +960,24 @@ let
 end
 ```
 
+#### Infiltration from surface water
+
+Infiltration from surface water (overland flow) is supported when enabling the following
+setting in the TOML file:
+
+```toml
+[model]
+reinfiltration_surfacewater__flag = true
+```
+
+This can be done with both `kinematic_wave` and `local_inertial` land routing settings. When
+enabling this proces, the amount of overland flow is added as water that can potentially
+infiltrate. During infiltration, the amount of infiltrated water originating from overland
+flow is tracked using a simple ratio
+(`available_from_overland_flow/total_potential_infiltration`). After the infiltration
+process, the volume of water in the overland flow is corrected to remove this water from the
+overland flow.
+
 ### Capillary rise
 
 Capillary rise is determined when an unsaturated zone occurs in the soil column, using the

From b8767379f37d201a7499ac73303ee0ae862de0ab Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 2 Apr 2026 09:57:40 +0200
Subject: [PATCH 15/24] fix typo

---
 docs/model_docs/land/landhydrology_sbm.qmd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index 1f16b2228..6eb308550 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -971,7 +971,7 @@ reinfiltration_surfacewater__flag = true
 ```
 
 This can be done with both `kinematic_wave` and `local_inertial` land routing settings. When
-enabling this proces, the amount of overland flow is added as water that can potentially
+enabling this process, the amount of overland flow is added as water that can potentially
 infiltrate. During infiltration, the amount of infiltrated water originating from overland
 flow is tracked using a simple ratio
 (`available_from_overland_flow/total_potential_infiltration`). After the infiltration

From 8828dc0cd74acb09729631274cc606b5a7ab129b Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 10 Apr 2026 16:51:26 +0200
Subject: [PATCH 16/24] update names; updated and added tests

---
 Wflow/src/soil/soil.jl         | 33 ++++++++-------
 Wflow/src/soil/soil_process.jl | 24 +++++++----
 Wflow/test/land_process.jl     | 76 ++++++++++++++++++++++++----------
 Wflow/test/run_sbm.jl          | 72 ++++++++++++++++++++++++++++++++
 4 files changed, 159 insertions(+), 46 deletions(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index f446e3d62..02a45e437 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -45,10 +45,10 @@ abstract type AbstractSoilModel end
     infiltexcess::Vector{Float64} = fill(MISSING_VALUE, n)
     # Infiltration from surface water [mm Δt⁻¹]
     infilt_surfacewater::Vector{Float64} = fill(0.0, n)
-    # Contribution from surface water [mm Δt⁻¹]
-    infilt_available_surfacewater::Vector{Float64} = fill(0.0, n)
+    # Potential infiltration originating from surface water [mm Δt⁻¹]
+    potential_infiltration_surfacewater::Vector{Float64} = fill(0.0, n)
     # Total water available for infiltration [mm Δt⁻¹]
-    infilt_total_available::Vector{Float64} = fill(0.0, n)
+    potential_infiltration::Vector{Float64} = fill(0.0, n)
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{Float64} = fill(MISSING_VALUE, n)
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -693,12 +693,13 @@ function update_bc_soil_model!(
     evaporation!(demand.paddy, potential_soilevaporation)
     potential_soilevaporation .= potential_soilevaporation .- get_evaporation(demand.paddy)
 
-    water_flux_surface .= max.(
-        runoff.boundary_conditions.water_flux_surface .+
-        get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
-        runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
-        0.0,
-    )
+    water_flux_surface .=
+        max.(
+            runoff.boundary_conditions.water_flux_surface .+
+            get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
+            runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
+            0.0,
+        )
     return nothing
 end
 
@@ -758,15 +759,15 @@ function update_available_for_infiltration!(
     (; waterdepth_land) = runoff.boundary_conditions
     (; river_fraction) = domain.land.parameters
 
-    n = length(v.infilt_total_available)
+    n = length(v.potential_infiltration)
     threaded_foreach(1:n; basesize = 1000) do i
-        v.infilt_available_surfacewater[i] = 0.0
+        v.potential_infiltration_surfacewater[i] = 0.0
         if do_surface_water_infiltration
-            v.infilt_available_surfacewater[i] =
+            v.potential_infiltration_surfacewater[i] =
                 waterdepth_land[i] * (1.0 - river_fraction[i]) * 0.95
-            water_flux_surface[i] += v.infilt_available_surfacewater[i]
+            water_flux_surface[i] += v.potential_infiltration_surfacewater[i]
         end
-        v.infilt_total_available[i] = water_flux_surface[i]
+        v.potential_infiltration[i] = water_flux_surface[i]
     end
 
     return nothing
@@ -783,8 +784,8 @@ function correct_infiltration!(model::SbmSoilModel)
         v.infiltexcess[i],
         v.excesswater[i],
         water_flux_surface[i] = correct_infiltration(
-            v.infilt_total_available[i],
-            v.infilt_available_surfacewater[i],
+            v.potential_infiltration[i],
+            v.potential_infiltration_surfacewater[i],
             water_flux_surface[i],
             v.actinfilt[i],
             v.infiltexcess[i],
diff --git a/Wflow/src/soil/soil_process.jl b/Wflow/src/soil/soil_process.jl
index 54e057f00..5b132ae24 100644
--- a/Wflow/src/soil/soil_process.jl
+++ b/Wflow/src/soil/soil_process.jl
@@ -268,28 +268,38 @@ function actual_infiltration_soil_path(
     return actinfiltsoil, actinfiltpath
 end
 
+"""
+Correct infiltration fluxes by separating the surface water contribution from the total
+infiltration. The correction factor is based on the ratio of
+`potential_infiltration_surfacewater` to `potential_infiltration`, and is applied to
+`actinfilt` and `infiltexcess` to remove the surface water component. The surface water flux
+`water_flux_surface` is adjusted accordingly, and the remaining excess water is computed.
+
+Returns `infilt_surfacewater`, corrected `actinfilt`, corrected `infiltexcess`, `excesswater`,
+and corrected `water_flux_surface`.
+"""
 function correct_infiltration(
-    infilt_total_available,
-    infilt_available_surfacewater,
+    potential_infiltration,
+    potential_infiltration_surfacewater,
     water_flux_surface,
     actinfilt,
     infiltexcess,
 )
     # Determine ratio of water that has infiltrated
-    infilt_ratio = infilt_total_available == 0.0 ? 0.0 : actinfilt / infilt_total_available
+    infilt_ratio = potential_infiltration == 0.0 ? 0.0 : actinfilt / potential_infiltration
     # Use this ratio to determine the contribution from overland flow
-    infilt_surfacewater = max(0.0, infilt_available_surfacewater * infilt_ratio)
+    infilt_surfacewater = max(0.0, potential_infiltration_surfacewater * infilt_ratio)
     # Determine the correction factor to apply to the relevant fluxes
     correction_surfacewater =
-        infilt_total_available == 0.0 ? 1.0 :
-        1.0 - (infilt_available_surfacewater / infilt_total_available)
+        potential_infiltration == 0.0 ? 1.0 :
+        1.0 - (potential_infiltration_surfacewater / potential_infiltration)
 
     # Correct fluxes
     actinfilt *= correction_surfacewater
     infiltexcess *= correction_surfacewater
 
     # subtract contribution from overland flow to ensure correct fluxes
-    water_flux_surface -= infilt_available_surfacewater
+    water_flux_surface -= potential_infiltration_surfacewater
     excesswater = water_flux_surface - actinfilt - infiltexcess
 
     return infilt_surfacewater, actinfilt, infiltexcess, excesswater, water_flux_surface
diff --git a/Wflow/test/land_process.jl b/Wflow/test/land_process.jl
index af1149298..9534f39e0 100644
--- a/Wflow/test/land_process.jl
+++ b/Wflow/test/land_process.jl
@@ -447,11 +447,12 @@ end
 end
 
 @testitem "unit: correct_infiltration" begin
+    # Test with infiltration from surface water and precipitation
     potential_infiltration = 10.0
-    infilt_available_surfacewater = 2.0
-    water_flux_surface = 10.0
-    actual_infiltration = 6.0
-    infiltexcess = 1.0
+    potential_infiltration_surfacewater = 2.0
+    water_flux_surface_input = 10.0
+    actual_infiltration_input = 6.0 # this includes the infiltration from surface water
+    infiltexcess_input = 1.0
 
     infilt_surfacewater,
     actual_infiltration,
@@ -459,23 +460,25 @@ end
     excesswater,
     water_flux_surface = Wflow.correct_infiltration(
         potential_infiltration,
-        infilt_available_surfacewater,
-        water_flux_surface,
-        actual_infiltration,
-        infiltexcess,
+        potential_infiltration_surfacewater,
+        water_flux_surface_input,
+        actual_infiltration_input,
+        infiltexcess_input,
     )
 
-    @test infilt_surfacewater ≈ 1.2
-    @test actual_infiltration ≈ 4.8
-    @test infiltexcess ≈ 0.8
-    @test water_flux_surface ≈ 8.0
+    @test infilt_surfacewater == 1.2 # this is the infiltration from surface water
+    @test actual_infiltration ≈ 4.8 # this excludes the infiltration from surface water
+    @test infilt_surfacewater + actual_infiltration ≈ actual_infiltration_input
+    @test infiltexcess == 0.8
+    @test water_flux_surface == 8.0
     @test excesswater ≈ 2.4
 
-    potential_infiltration = 0.0
-    infilt_available_surfacewater = 1.0
-    water_flux_surface = 3.0
-    actual_infiltration = 0.0
-    infiltexcess = 0.0
+    # Test with infiltration from only precipitation, no infiltration from surface water
+    potential_infiltration = 10.0
+    potential_infiltration_surfacewater = 0.0
+    water_flux_surface_input = 10.0
+    actual_infiltration_input = 6.0
+    infiltexcess_input = 1.0
 
     infilt_surfacewater,
     actual_infiltration,
@@ -483,17 +486,44 @@ end
     excesswater,
     water_flux_surface = Wflow.correct_infiltration(
         potential_infiltration,
-        infilt_available_surfacewater,
-        water_flux_surface,
-        actual_infiltration,
-        infiltexcess,
+        potential_infiltration_surfacewater,
+        water_flux_surface_input,
+        actual_infiltration_input,
+        infiltexcess_input,
     )
 
     @test infilt_surfacewater == 0.0
+    @test actual_infiltration == actual_infiltration_input
+    @test infilt_surfacewater + actual_infiltration == actual_infiltration_input
+    @test infiltexcess == infiltexcess_input
+    @test water_flux_surface == water_flux_surface_input
+    @test excesswater == 3.0
+
+    # Test with infiltration from only surface water, no infiltration from precipitation
+    potential_infiltration = 10.0
+    potential_infiltration_surfacewater = 10.0
+    water_flux_surface_input = 10.0
+    actual_infiltration_input = 6.0
+    infiltexcess_input = 1.0
+
+    infilt_surfacewater,
+    actual_infiltration,
+    infiltexcess,
+    excesswater,
+    water_flux_surface = Wflow.correct_infiltration(
+        potential_infiltration,
+        potential_infiltration_surfacewater,
+        water_flux_surface_input,
+        actual_infiltration_input,
+        infiltexcess_input,
+    )
+
+    @test infilt_surfacewater == 6.0
     @test actual_infiltration == 0.0
+    @test infilt_surfacewater + actual_infiltration == 6.0
     @test infiltexcess == 0.0
-    @test water_flux_surface ≈ 2.0
-    @test excesswater ≈ 2.0
+    @test water_flux_surface == 0.0
+    @test excesswater == 0.0
 end
 
 @testitem "unit: correct_overland_flow_level" begin
diff --git a/Wflow/test/run_sbm.jl b/Wflow/test/run_sbm.jl
index 58cbe2ce2..2184a1faa 100644
--- a/Wflow/test/run_sbm.jl
+++ b/Wflow/test/run_sbm.jl
@@ -823,3 +823,75 @@ end
     end
     Wflow.close_files(model; delete_output = false)
 end
+
+@testitem "reinfiltration" begin
+    @testset "Kinematic wave overland flow" begin
+        tomlpath = joinpath(@__DIR__, "sbm_config.toml")
+        config = Wflow.Config(tomlpath)
+        config.dir_output = mktempdir()
+        config.model.reinfiltration_surfacewater__flag = true
+
+        idxs = [53, 54, 55]
+
+        model = Wflow.Model(config)
+        Wflow.run_timestep!(model)
+        Wflow.run_timestep!(model)
+
+        (; soil) = model.land
+        # get total available infiltration
+        @test soil.variables.potential_infiltration[idxs] ≈
+              [0.07163384298140926, 12.55529617902951, 1.094994363241204]
+        @test soil.variables.potential_infiltration_surfacewater[idxs] ≈
+              [0.0, 12.485540999074937, 1.0244892972592872]
+        # get actual infiltration
+        @test soil.variables.actinfilt[idxs] ≈
+              [0.07163384298140926, 0.06975517995457337, 0.07050506598191691]
+        # there is a lot of infiltration coming from surface water in cell 54
+        @test soil.variables.infilt_surfacewater[idxs] ≈
+              [0.0, 12.485540999074937, 1.0244892972592872]
+
+        (; h) = model.routing.overland_flow.variables
+        @test h[idxs] ≈ [0.0, 0.039396885875137864, 0.02395111033358323]
+
+        Wflow.close_files(model; delete_output = false)
+    end
+
+    @testset "Local inertial overland flow" begin
+        tomlpath = joinpath(@__DIR__, "sbm_river-land-local-inertial_config.toml")
+        config = Wflow.Config(tomlpath)
+        config.dir_output = mktempdir()
+        config.model.reinfiltration_surfacewater__flag = true
+
+        idxs = [13929, 13930, 13931]
+
+        model = Wflow.Model(config)
+        Wflow.run_timestep!(model)
+        old_h = copy(model.routing.overland_flow.variables.h)
+        Wflow.run_timestep!(model)
+
+        (; soil) = model.land
+        # get total available infiltration
+        @test soil.variables.potential_infiltration[idxs] ≈
+              [0.16281059646400808, 1.9744062829658247, 0.16363083612044332]
+        @test soil.variables.potential_infiltration_surfacewater[idxs] ≈
+              [0.007480075393492462, 1.9744062829658247, 0.0]
+        # get actual infiltration
+        @test soil.variables.actinfilt[idxs] ≈
+              [0.15533052107051562, 0.0, 0.16363083612044332]
+        # there is a lot of infiltration coming from surface water in cell 54
+        @test soil.variables.infilt_surfacewater[idxs] ≈
+              [0.007480075393492462, 1.9744062829658247, 0.0]
+
+        (; h) = model.routing.overland_flow.variables
+        # all available surface water was infiltrated
+        @test h[idxs] ≈ [0.0, 0.0, 0.0]
+        decreased_h = h[idxs] .< old_h[idxs]
+        equal_h = h[idxs] .== old_h[idxs]
+        # test that h decreased in the first two cells, but not in the third cell (had no
+        # surface water infiltration)
+        @test decreased_h == [true, true, false]
+        @test equal_h == [false, false, true]
+
+        Wflow.close_files(model; delete_output = false)
+    end
+end

From 71964b5138bbe42627b9ebb6f235f57cfee8e68a Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 23 Jul 2026 16:11:10 +0200
Subject: [PATCH 17/24] process comments part1

---
 Wflow/src/config_structure.jl                 |  2 +-
 Wflow/src/routing/surface/surface_kinwave.jl  | 74 +++++++++++-----
 .../routing/surface/surface_local_inertial.jl | 29 +++---
 Wflow/src/routing/surface/surface_routing.jl  |  5 ++
 Wflow/src/soil/soil.jl                        | 53 +++++------
 Wflow/src/soil/soil_process.jl                | 32 +------
 Wflow/test/land_process.jl                    | 88 ++++++++++++-------
 Wflow/test/run_sbm.jl                         | 40 ++++++---
 Wflow/test/soil.jl                            |  2 +-
 docs/model_docs/land/landhydrology_sbm.qmd    |  8 +-
 10 files changed, 183 insertions(+), 150 deletions(-)

diff --git a/Wflow/src/config_structure.jl b/Wflow/src/config_structure.jl
index 050ec666e..740720a4d 100644
--- a/Wflow/src/config_structure.jl
+++ b/Wflow/src/config_structure.jl
@@ -70,7 +70,7 @@ end
     snow_gravitational_transport__flag::Bool = false
     glacier__flag::Bool = false
     soil_infiltration_reduction__flag::Bool = false
-    reinfiltration_surfacewater__flag::Bool = false
+    land_surface_water_reinfiltration__flag::Bool = false
     soil_layer__thickness::Vector{Int} = [100, 300, 800]
     saturated_hydraulic_conductivity_profile::VerticalConductivityProfile.T =
         VerticalConductivityProfile.exponential
diff --git a/Wflow/src/routing/surface/surface_kinwave.jl b/Wflow/src/routing/surface/surface_kinwave.jl
index a8afaff0a..71a910664 100644
--- a/Wflow/src/routing/surface/surface_kinwave.jl
+++ b/Wflow/src/routing/surface/surface_kinwave.jl
@@ -654,30 +654,58 @@ get_inflow_reservoir(::AbstractRiverFlowModel, ::GroundwaterFlowModel, inds::Vec
 Update overland flow water level and discharge for KinWaveOverlandFlow model based on
 surface water infiltration.
 """
-function correct_overland_flow_level!(
-    model::SbmSoilModel,
-    overland_flow::KinWaveOverlandFlowModel,
+function update_overland_flow_and_depth!(
+    overland_flow_model::KinWaveOverlandFlowModel,
+    soil_model::SbmSoilModel,
     domain::Domain,
-    config::Config,
 )
-    v = model.variables
-
-    if config.model.reinfiltration_surfacewater__flag
-        (; surface_flow_width) = domain.land.parameters
-        n = length(surface_flow_width)
-        threaded_foreach(1:n; basesize = 1000) do i
-            q, h = correct_overland_flow_level(
-                overland_flow.variables.h[i],
-                v.infilt_surfacewater[i],
-                domain.land.parameters.river_fraction[i],
-                surface_flow_width[i],
-                overland_flow.parameters.alpha[i],
-                overland_flow.parameters.beta,
-            )
-            if !isnothing(q)
-                overland_flow.variables.flow.q[i] = q
-                overland_flow.variables.h[i] = h
-            end
-        end
+    (; infilt_surfacewater) = soil_model.variables
+    n = length(infilt_surfacewater)
+    threaded_foreach(1:n; basesize = 1000) do i
+        update_overland_flow_and_depth!(
+            overland_flow_model,
+            infilt_surfacewater[i],
+            domain.land.parameters,
+            i,
+        )
     end
 end
+
+"""
+Update overland flow water level and discharge in-place for a single cell based on
+surface water infiltration.
+"""
+function update_overland_flow_and_depth!(
+    overland_flow_model::KinWaveOverlandFlowModel,
+    infilt_surfacewater,
+    land_parameters,
+    i,
+)
+    if infilt_surfacewater > 0.0
+        # Get original h_land in mm
+        original_h_land = overland_flow_model.variables.h[i] * 1000.0
+
+        # Correct values for river fraction to ensure correct water accounting
+        infiltrated_surfacewater =
+            (infilt_surfacewater / (1.0 - land_parameters.river_fraction[i]))
+        # Calculate new h_land in m
+        h = (original_h_land - infiltrated_surfacewater) / 1000.0
+
+        q = ifelse(
+            land_parameters.surface_flow_width[i] > 0.0 && h > 0.0,
+            # Compute cross-sectional area from h
+            pow(
+                (h * land_parameters.surface_flow_width[i]) /
+                overland_flow_model.parameters.alpha[i],
+                1.0 / overland_flow_model.parameters.beta,
+            ),
+            0.0,
+        )
+        # set q to 0.0 if it is below the minimum flow threshold
+        q = ifelse(q < KIN_WAVE_MIN_FLOW, 0.0, q)
+
+        overland_flow_model.variables.flow.q[i] = q
+        overland_flow_model.variables.h[i] = h
+    end
+    return nothing
+end
diff --git a/Wflow/src/routing/surface/surface_local_inertial.jl b/Wflow/src/routing/surface/surface_local_inertial.jl
index 18209fe3e..6124c1325 100644
--- a/Wflow/src/routing/surface/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface/surface_local_inertial.jl
@@ -1023,8 +1023,7 @@ function local_inertial_update_water_depth!(
                 ) * dt
             # Apply surface water infiltration correction for river cells
             if land_bc.infiltration_volume[i] > 0.0
-                land_v.storage[i] =
-                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
+                land_v.storage[i] -= land_bc.infiltration_volume[i]
             end
             if land_v.storage[i] < 0.0
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
@@ -1072,15 +1071,14 @@ function local_inertial_update_water_depth!(
                     land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] +
                     land_bc.runoff[i]
                 ) * dt
+            if land_bc.infiltration_volume[i] > 0.0
+                land_v.storage[i] -= land_bc.infiltration_volume[i]
+            end
             if land_v.storage[i] < 0.0
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
             end
             # Apply surface water infiltration correction if available
-            if land_bc.infiltration_volume[i] > 0.0
-                land_v.storage[i] =
-                    max(0.0, land_v.storage[i] - land_bc.infiltration_volume[i])
-            end
             land_v.h[i] = land_v.storage[i] / (x_length[i] * y_length[i])
         end
     end
@@ -1370,19 +1368,16 @@ Update overland flow water level and storage for LocalInertialOverlandFlow model
 surface water infiltration. For local inertial flow, the infiltration is applied by updating
 the boundary conditions which are then used in the water depth update function.
 """
-function correct_overland_flow_level!(
-    model::SbmSoilModel,
-    overland_flow::LocalInertialOverlandFlowModel,
+function update_overland_flow_and_depth!(
+    overland_flow_model::LocalInertialOverlandFlowModel,
+    soil_model::SbmSoilModel,
     domain::Domain,
-    config::Config,
 )
-    (; infilt_surfacewater) = model.variables
+    (; infilt_surfacewater) = soil_model.variables
     (; area) = domain.land.parameters
 
-    if config.model.reinfiltration_surfacewater__flag
-        # Update the boundary condition for surface water infiltration
-        # This will be used in local_inertial_update_water_depth!
-        overland_flow.boundary_conditions.infiltration_volume .=
-            infilt_surfacewater .* area .* 0.001
-    end
+    # Update the boundary condition for surface water infiltration
+    # This will be used in local_inertial_update_water_depth!
+    overland_flow_model.boundary_conditions.infiltration_volume .=
+        infilt_surfacewater .* area .* 0.001
 end
diff --git a/Wflow/src/routing/surface/surface_routing.jl b/Wflow/src/routing/surface/surface_routing.jl
index 5ea734167..fd3e510b4 100644
--- a/Wflow/src/routing/surface/surface_routing.jl
+++ b/Wflow/src/routing/surface/surface_routing.jl
@@ -10,6 +10,11 @@ function surface_routing!(model)
     (; overland_flow, river_flow, subsurface_flow) = routing
     (; reservoir) = river_flow.boundary_conditions
 
+    # correct overland flow water levels in case of reinfiltration
+    if config.model.land_surface_water_reinfiltration__flag
+        update_overland_flow_and_depth!(overland_flow, soil, domain)
+    end
+
     dt = tosecond(clock.dt)
     # update lateral inflow for kinematic wave overland flow
     update_lateral_inflow!(
diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 02a45e437..6b432eb51 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -693,13 +693,12 @@ function update_bc_soil_model!(
     evaporation!(demand.paddy, potential_soilevaporation)
     potential_soilevaporation .= potential_soilevaporation .- get_evaporation(demand.paddy)
 
-    water_flux_surface .=
-        max.(
-            runoff.boundary_conditions.water_flux_surface .+
-            get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
-            runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
-            0.0,
-        )
+    water_flux_surface .= max.(
+        runoff.boundary_conditions.water_flux_surface .+
+        get_irrigation_allocated(allocation) .- runoff.variables.runoff_river .-
+        runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
+        0.0,
+    )
     return nothing
 end
 
@@ -773,22 +772,29 @@ function update_available_for_infiltration!(
     return nothing
 end
 
-function correct_infiltration!(model::SbmSoilModel)
-    v = model.variables
-    (; water_flux_surface) = model.boundary_conditions
+function update_infiltration_fluxes!(soil_model::SbmSoilModel)
+    (;
+        infilt_surfacewater,
+        actinfilt,
+        infiltexcess,
+        excesswater,
+        potential_infiltration,
+        potential_infiltration_surfacewater,
+    ) = soil_model.variables
+    (; water_flux_surface) = soil_model.boundary_conditions
 
-    n = length(v.actinfilt)
+    n = length(actinfilt)
     threaded_foreach(1:n; basesize = 1000) do i
-        v.infilt_surfacewater[i],
-        v.actinfilt[i],
-        v.infiltexcess[i],
-        v.excesswater[i],
-        water_flux_surface[i] = correct_infiltration(
-            v.potential_infiltration[i],
-            v.potential_infiltration_surfacewater[i],
+        infilt_surfacewater[i],
+        actinfilt[i],
+        infiltexcess[i],
+        excesswater[i],
+        water_flux_surface[i] = update_infiltration_fluxes(
+            potential_infiltration[i],
+            potential_infiltration_surfacewater[i],
             water_flux_surface[i],
-            v.actinfilt[i],
-            v.infiltexcess[i],
+            actinfilt[i],
+            infiltexcess[i],
         )
     end
 end
@@ -1196,7 +1202,7 @@ function update_soil_water_flow!(
         soil_model,
         domain,
         runoff,
-        config.model.reinfiltration_surfacewater__flag,
+        config.model.land_surface_water_reinfiltration__flag,
     )
 
     infiltration!(soil_model)
@@ -1210,7 +1216,7 @@ function update_soil_water_flow!(
 
     # Correct fluxes in case of reinfiltration, also to ensure correct soil and path
     # infiltration, and excesswater
-    correct_infiltration!(soil_model)
+    update_infiltration_fluxes!(soil_model)
 
     actual_infiltration_soil_path!(soil_model)
 
@@ -1392,9 +1398,6 @@ function update_soil_water_storage!(
     update_runoff!(demand.paddy, v.runoff)
     @. v.net_runoff = v.runoff - ae_openw_l
 
-    # correct overland flow water levels in case of reinfiltration
-    correct_overland_flow_level!(soil_model, overland_flow, domain, config)
-
     return nothing
 end
 
diff --git a/Wflow/src/soil/soil_process.jl b/Wflow/src/soil/soil_process.jl
index 5b132ae24..13395fd1b 100644
--- a/Wflow/src/soil/soil_process.jl
+++ b/Wflow/src/soil/soil_process.jl
@@ -278,7 +278,7 @@ infiltration. The correction factor is based on the ratio of
 Returns `infilt_surfacewater`, corrected `actinfilt`, corrected `infiltexcess`, `excesswater`,
 and corrected `water_flux_surface`.
 """
-function correct_infiltration(
+function update_infiltration_fluxes(
     potential_infiltration,
     potential_infiltration_surfacewater,
     water_flux_surface,
@@ -304,33 +304,3 @@ function correct_infiltration(
 
     return infilt_surfacewater, actinfilt, infiltexcess, excesswater, water_flux_surface
 end
-
-function correct_overland_flow_level(
-    overlandflow_depth,
-    infilt_surfacewater,
-    river_fraction,
-    surface_flow_width,
-    alpha,
-    beta,
-)
-    if infilt_surfacewater > 0.0
-        # Get original h_land in mm
-        original_h_land = overlandflow_depth * 1000.0
-
-        # Correct values for river fraction to ensure correct water accounting
-        infiltrated_surfacewater = (infilt_surfacewater / (1.0 - river_fraction))
-        # Calculate new h_land in m
-        h = (original_h_land - infiltrated_surfacewater) / 1000.0
-
-        q = ifelse(
-            surface_flow_width > 0.0 && h > 0.0 && alpha > 0.0 && beta != 0.0,
-            # Compute cross-sectional area from h
-            pow.(max.((h * surface_flow_width) / alpha, 1e-10), 1.0 / beta),
-            0.0,  # Set q to 0.0 if conditions are not met
-        )
-    else
-        q = nothing
-        h = overlandflow_depth
-    end
-    return q, h
-end
diff --git a/Wflow/test/land_process.jl b/Wflow/test/land_process.jl
index 9534f39e0..a245ad38e 100644
--- a/Wflow/test/land_process.jl
+++ b/Wflow/test/land_process.jl
@@ -446,7 +446,7 @@ end
     @test actinfiltpath == 0.0
 end
 
-@testitem "unit: correct_infiltration" begin
+@testitem "unit: update_infiltration_fluxes" begin
     # Test with infiltration from surface water and precipitation
     potential_infiltration = 10.0
     potential_infiltration_surfacewater = 2.0
@@ -458,7 +458,7 @@ end
     actual_infiltration,
     infiltexcess,
     excesswater,
-    water_flux_surface = Wflow.correct_infiltration(
+    water_flux_surface = Wflow.update_infiltration_fluxes(
         potential_infiltration,
         potential_infiltration_surfacewater,
         water_flux_surface_input,
@@ -484,7 +484,7 @@ end
     actual_infiltration,
     infiltexcess,
     excesswater,
-    water_flux_surface = Wflow.correct_infiltration(
+    water_flux_surface = Wflow.update_infiltration_fluxes(
         potential_infiltration,
         potential_infiltration_surfacewater,
         water_flux_surface_input,
@@ -510,7 +510,7 @@ end
     actual_infiltration,
     infiltexcess,
     excesswater,
-    water_flux_surface = Wflow.correct_infiltration(
+    water_flux_surface = Wflow.update_infiltration_fluxes(
         potential_infiltration,
         potential_infiltration_surfacewater,
         water_flux_surface_input,
@@ -526,36 +526,56 @@ end
     @test excesswater == 0.0
 end
 
-@testitem "unit: correct_overland_flow_level" begin
-    overlandflow_depth = 0.02
-    infilt_surfacewater = 5.0
+@testitem "unit: update_overland_flow_and_depth!" begin
+    using Wflow:
+        KinWaveOverlandFlowModel,
+        ManningFlowParameters,
+        OverLandFlowVariables,
+        FlowVariables,
+        LandFlowBC,
+        TimeStepping
+
+    n = 1
+    infiltration_amount = 5.0 # mm
+    original_depth = 0.02 # m
     river_fraction = 0.2
-    surface_flow_width = 10.0
-    alpha = 2.0
-    beta = 0.5
-
-    q, h = Wflow.correct_overland_flow_level(
-        overlandflow_depth,
-        infilt_surfacewater,
-        river_fraction,
-        surface_flow_width,
-        alpha,
-        beta,
+    expected_water_depth =
+        original_depth - ((infiltration_amount * 1e-3) / (1 - river_fraction))
+
+    flow_vars = FlowVariables(; n)
+    flow_vars.q[1] = 0.0
+    variables = OverLandFlowVariables(; n, flow = flow_vars)
+    variables.h[1] = original_depth
+
+    mannings_n = [0.072]
+    slope = [0.01]
+    parameters = ManningFlowParameters(mannings_n, slope)
+    parameters.alpha[1] = 2.0
+
+    boundary_conditions = LandFlowBC(; n)
+    timestepping =
+        TimeStepping(; adaptive = false, dt_fixed = 900.0, stable_timesteps = zeros(n))
+
+    overland_flow_model =
+        KinWaveOverlandFlowModel(; timestepping, boundary_conditions, parameters, variables)
+
+    land_parameters = (; river_fraction = [river_fraction], surface_flow_width = [10.0])
+
+    # Test with positive infiltration
+    Wflow.update_overland_flow_and_depth!(
+        overland_flow_model,
+        infiltration_amount,
+        land_parameters,
+        1,
     )
-
-    @test h ≈ 0.01375
-    @test q ≈ 0.0047265625
-
-    infilt_surfacewater = 0.0
-    q, h = Wflow.correct_overland_flow_level(
-        overlandflow_depth,
-        infilt_surfacewater,
-        river_fraction,
-        surface_flow_width,
-        alpha,
-        beta,
-    )
-
-    @test q === nothing
-    @test h == overlandflow_depth
+    @test overland_flow_model.variables.h[1] ≈ expected_water_depth
+    @test overland_flow_model.variables.h[1] ≈ 0.01375
+    @test overland_flow_model.variables.flow.q[1] ≈ 0.011537751232883156
+
+    # Test with zero infiltration (no update should occur)
+    overland_flow_model.variables.h[1] = original_depth
+    overland_flow_model.variables.flow.q[1] = 0.1
+    Wflow.update_overland_flow_and_depth!(overland_flow_model, 0.0, land_parameters, 1)
+    @test overland_flow_model.variables.flow.q[1] == 0.1
+    @test overland_flow_model.variables.h[1] == original_depth
 end
diff --git a/Wflow/test/run_sbm.jl b/Wflow/test/run_sbm.jl
index 2184a1faa..d0a7f203b 100644
--- a/Wflow/test/run_sbm.jl
+++ b/Wflow/test/run_sbm.jl
@@ -829,12 +829,13 @@ end
         tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         config = Wflow.Config(tomlpath)
         config.dir_output = mktempdir()
-        config.model.reinfiltration_surfacewater__flag = true
+        config.model.land_surface_water_reinfiltration__flag = true
 
         idxs = [53, 54, 55]
 
         model = Wflow.Model(config)
         Wflow.run_timestep!(model)
+        # old_h = copy(model.routing.overland_flow.variables.h)
         Wflow.run_timestep!(model)
 
         (; soil) = model.land
@@ -850,8 +851,17 @@ end
         @test soil.variables.infilt_surfacewater[idxs] ≈
               [0.0, 12.485540999074937, 1.0244892972592872]
 
-        (; h) = model.routing.overland_flow.variables
-        @test h[idxs] ≈ [0.0, 0.039396885875137864, 0.02395111033358323]
+        # (; h) = model.routing.overland_flow.variables
+        # @test h[idxs] ≈ [0.0, 0.039396885875137864, 0.02395111033358323]
+        # @test h[idxs] - old_h[idxs] ≈ [0.0, -0.039396885875137864, -0.02395111033358323]
+        # # test that h decreased in the first two cells, but not in the third cell (had no
+        # # surface water infiltration)
+        # decreased_h = h[idxs] .< old_h[idxs]
+        # equal_h = h[idxs] .== old_h[idxs]
+        # increased_h = h[idxs] .> old_h[idxs]
+        # @test decreased_h == [false, true, true]
+        # @test equal_h == [true, false, false]
+        # @test increased_h == [false, false, false]
 
         Wflow.close_files(model; delete_output = false)
     end
@@ -860,13 +870,13 @@ end
         tomlpath = joinpath(@__DIR__, "sbm_river-land-local-inertial_config.toml")
         config = Wflow.Config(tomlpath)
         config.dir_output = mktempdir()
-        config.model.reinfiltration_surfacewater__flag = true
+        config.model.land_surface_water_reinfiltration__flag = true
 
         idxs = [13929, 13930, 13931]
 
         model = Wflow.Model(config)
         Wflow.run_timestep!(model)
-        old_h = copy(model.routing.overland_flow.variables.h)
+        # old_h = copy(model.routing.overland_flow.variables.h)
         Wflow.run_timestep!(model)
 
         (; soil) = model.land
@@ -882,15 +892,17 @@ end
         @test soil.variables.infilt_surfacewater[idxs] ≈
               [0.007480075393492462, 1.9744062829658247, 0.0]
 
-        (; h) = model.routing.overland_flow.variables
-        # all available surface water was infiltrated
-        @test h[idxs] ≈ [0.0, 0.0, 0.0]
-        decreased_h = h[idxs] .< old_h[idxs]
-        equal_h = h[idxs] .== old_h[idxs]
-        # test that h decreased in the first two cells, but not in the third cell (had no
-        # surface water infiltration)
-        @test decreased_h == [true, true, false]
-        @test equal_h == [false, false, true]
+        # (; h) = model.routing.overland_flow.variables
+        # # all available surface water was infiltrated
+        # @test h[idxs] ≈ [0.0, 0.0, 0.0]
+        # decreased_h = h[idxs] .< old_h[idxs]
+        # equal_h = h[idxs] .== old_h[idxs]
+        # increased_h = h[idxs] .> old_h[idxs]
+        # # test that h decreased in the first two cells, but not in the third cell (had no
+        # # surface water infiltration)
+        # @test decreased_h == [true, true, false]
+        # @test equal_h == [false, false, true]
+        # @test increased_h == [false, false, false]
 
         Wflow.close_files(model; delete_output = false)
     end
diff --git a/Wflow/test/soil.jl b/Wflow/test/soil.jl
index 44558ba6f..3c2d6f650 100644
--- a/Wflow/test/soil.jl
+++ b/Wflow/test/soil.jl
@@ -299,7 +299,7 @@ end
     config = Wflow.Config(;
         model = Wflow.ModelSection(;
             type = Wflow.ModelType.sbm,
-            reinfiltration_surfacewater__flag = false,
+            land_surface_water_reinfiltration__flag = false,
         ),
         input = Wflow.InputSection(;
             path_forcing = "",
diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index 6eb308550..af943de14 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -960,14 +960,14 @@ let
 end
 ```
 
-#### Infiltration from surface water
+#### Land surface water infiltration
 
-Infiltration from surface water (overland flow) is supported when enabling the following
-setting in the TOML file:
+Infiltration from land surface water (overland flow) is supported when enabling the
+following setting in the TOML file:
 
 ```toml
 [model]
-reinfiltration_surfacewater__flag = true
+land_surface_water_reinfiltration__flag = true
 ```
 
 This can be done with both `kinematic_wave` and `local_inertial` land routing settings. When

From 4a081fe19b9e9d4225b5ea5053c3c4dec7ae0bbe Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 23 Jul 2026 17:13:28 +0200
Subject: [PATCH 18/24] Update soil.jl

---
 Wflow/src/soil/soil.jl | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 6b432eb51..1b50702a5 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -1214,9 +1214,13 @@ function update_soil_water_flow!(
     # actual infiltration and excess water
     actual_infiltration!(soil_model)
 
-    # Correct fluxes in case of reinfiltration, also to ensure correct soil and path
-    # infiltration, and excesswater
-    update_infiltration_fluxes!(soil_model)
+    # Correct fluxes in case of reinfiltration, and only compute excess water if reinfiltration
+    # is not enabled
+    if config.model.land_surface_water_reinfiltration__flag
+        update_infiltration_fluxes!(soil_model)
+    else
+        @. v.excesswater = water_flux_surface - v.actinfilt - v.infiltexcess
+    end
 
     actual_infiltration_soil_path!(soil_model)
 

From c3243d1f15ad4524b366eec1ffbb7b1d11a24b8c Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 17:51:06 +0200
Subject: [PATCH 19/24] move infilt to boundary

---
 Wflow/src/soil/soil.jl | 26 +++++++++++++-------------
 Wflow/test/run_sbm.jl  |  8 ++++----
 2 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 1b50702a5..145f34aa6 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -45,10 +45,6 @@ abstract type AbstractSoilModel end
     infiltexcess::Vector{Float64} = fill(MISSING_VALUE, n)
     # Infiltration from surface water [mm Δt⁻¹]
     infilt_surfacewater::Vector{Float64} = fill(0.0, n)
-    # Potential infiltration originating from surface water [mm Δt⁻¹]
-    potential_infiltration_surfacewater::Vector{Float64} = fill(0.0, n)
-    # Total water available for infiltration [mm Δt⁻¹]
-    potential_infiltration::Vector{Float64} = fill(0.0, n)
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{Float64} = fill(MISSING_VALUE, n)
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -204,6 +200,10 @@ end
     potential_transpiration::Vector{Float64} = fill(MISSING_VALUE, n)
     # Potential soil evaporation rate [mm Δt⁻¹]
     potential_soilevaporation::Vector{Float64} = fill(MISSING_VALUE, n)
+    # Potential infiltration originating from surface water [mm Δt⁻¹]
+    potential_infiltration_surfacewater::Vector{Float64} = fill(0.0, n)
+    # Total water available for infiltration [mm Δt⁻¹]
+    potential_infiltration::Vector{Float64} = fill(0.0, n)
 end
 
 "Exponential depth profile of vertical hydraulic conductivity at the soil surface"
@@ -754,19 +754,20 @@ function update_available_for_infiltration!(
     do_surface_water_infiltration::Bool,
 )
     v = model.variables
-    (; water_flux_surface) = model.boundary_conditions
+    (; water_flux_surface, potential_infiltration_surfacewater, potential_infiltration) =
+        model.boundary_conditions
     (; waterdepth_land) = runoff.boundary_conditions
     (; river_fraction) = domain.land.parameters
 
-    n = length(v.potential_infiltration)
+    n = length(potential_infiltration)
     threaded_foreach(1:n; basesize = 1000) do i
-        v.potential_infiltration_surfacewater[i] = 0.0
+        potential_infiltration_surfacewater[i] = 0.0
         if do_surface_water_infiltration
-            v.potential_infiltration_surfacewater[i] =
+            potential_infiltration_surfacewater[i] =
                 waterdepth_land[i] * (1.0 - river_fraction[i]) * 0.95
-            water_flux_surface[i] += v.potential_infiltration_surfacewater[i]
+            water_flux_surface[i] += potential_infiltration_surfacewater[i]
         end
-        v.potential_infiltration[i] = water_flux_surface[i]
+        potential_infiltration[i] = water_flux_surface[i]
     end
 
     return nothing
@@ -778,10 +779,9 @@ function update_infiltration_fluxes!(soil_model::SbmSoilModel)
         actinfilt,
         infiltexcess,
         excesswater,
-        potential_infiltration,
-        potential_infiltration_surfacewater,
     ) = soil_model.variables
-    (; water_flux_surface) = soil_model.boundary_conditions
+    (; water_flux_surface, potential_infiltration, potential_infiltration_surfacewater) =
+        soil_model.boundary_conditions
 
     n = length(actinfilt)
     threaded_foreach(1:n; basesize = 1000) do i
diff --git a/Wflow/test/run_sbm.jl b/Wflow/test/run_sbm.jl
index d0a7f203b..07b7a4498 100644
--- a/Wflow/test/run_sbm.jl
+++ b/Wflow/test/run_sbm.jl
@@ -840,9 +840,9 @@ end
 
         (; soil) = model.land
         # get total available infiltration
-        @test soil.variables.potential_infiltration[idxs] ≈
+        @test soil.boundary_conditions.potential_infiltration[idxs] ≈
               [0.07163384298140926, 12.55529617902951, 1.094994363241204]
-        @test soil.variables.potential_infiltration_surfacewater[idxs] ≈
+        @test soil.boundary_conditions.potential_infiltration_surfacewater[idxs] ≈
               [0.0, 12.485540999074937, 1.0244892972592872]
         # get actual infiltration
         @test soil.variables.actinfilt[idxs] ≈
@@ -881,9 +881,9 @@ end
 
         (; soil) = model.land
         # get total available infiltration
-        @test soil.variables.potential_infiltration[idxs] ≈
+        @test soil.boundary_conditions.potential_infiltration[idxs] ≈
               [0.16281059646400808, 1.9744062829658247, 0.16363083612044332]
-        @test soil.variables.potential_infiltration_surfacewater[idxs] ≈
+        @test soil.boundary_conditions.potential_infiltration_surfacewater[idxs] ≈
               [0.007480075393492462, 1.9744062829658247, 0.0]
         # get actual infiltration
         @test soil.variables.actinfilt[idxs] ≈

From 02542daffddd0d5a387d1a818eda601061d79177 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 18:03:08 +0200
Subject: [PATCH 20/24] move to boundary conditions

---
 Wflow/src/sbm.jl       |  2 ++
 Wflow/src/soil/soil.jl | 24 ++++++++++--------------
 Wflow/test/soil.jl     | 33 ++++++++++++++++++++++++++++++++-
 3 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/Wflow/src/sbm.jl b/Wflow/src/sbm.jl
index 5ee580aaf..cbf767294 100644
--- a/Wflow/src/sbm.jl
+++ b/Wflow/src/sbm.jl
@@ -123,6 +123,8 @@ function update_land_hydrology_model!(
         soil,
         atmospheric_forcing,
         (; interception, runoff, demand, allocation),
+        domain,
+        config,
     )
 
     update_soil_water_flow!(
diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index 145f34aa6..c807e2882 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -681,6 +681,8 @@ function update_bc_soil_model!(
     soil_model::SbmSoilModel,
     atmospheric_forcing::AtmosphericForcing,
     external_models::NamedTuple,
+    domain::Domain,
+    config::Config,
 )
     (; interception, runoff, demand, allocation) = external_models
     (; potential_transpiration, water_flux_surface, potential_soilevaporation) =
@@ -699,6 +701,13 @@ function update_bc_soil_model!(
         runoff.variables.runoff_land .+ get_water_depth(demand.paddy),
         0.0,
     )
+    # update available for infiltration in case surface water infiltration is enabled
+    update_available_for_infiltration!(
+        soil_model,
+        domain,
+        runoff,
+        config.model.land_surface_water_reinfiltration__flag,
+    )
     return nothing
 end
 
@@ -774,12 +783,7 @@ function update_available_for_infiltration!(
 end
 
 function update_infiltration_fluxes!(soil_model::SbmSoilModel)
-    (;
-        infilt_surfacewater,
-        actinfilt,
-        infiltexcess,
-        excesswater,
-    ) = soil_model.variables
+    (; infilt_surfacewater, actinfilt, infiltexcess, excesswater) = soil_model.variables
     (; water_flux_surface, potential_infiltration, potential_infiltration_surfacewater) =
         soil_model.boundary_conditions
 
@@ -1197,14 +1201,6 @@ function update_soil_water_flow!(
         modelsnow = config.model.snow__flag,
         soil_infiltration_reduction = config.model.soil_infiltration_reduction__flag,
     )
-    # update available for infiltration in case surface water infiltration is enabled
-    update_available_for_infiltration!(
-        soil_model,
-        domain,
-        runoff,
-        config.model.land_surface_water_reinfiltration__flag,
-    )
-
     infiltration!(soil_model)
     # unsaturated zone flow
     unsaturated_zone_flow!(soil_model)
diff --git a/Wflow/test/soil.jl b/Wflow/test/soil.jl
index 3c2d6f650..429f861eb 100644
--- a/Wflow/test/soil.jl
+++ b/Wflow/test/soil.jl
@@ -63,11 +63,42 @@
 
     external_models = (; interception, runoff, demand, allocation)
 
-    Wflow.update_bc_soil_model!(soil_model, atmospheric_forcing, external_models)
+    domain = Wflow.Domain(;
+        land = Wflow.DomainLand(;
+            parameters = Wflow.LandParameters(; river_fraction = [0.1]),
+        ),
+    )
+    config = Wflow.Config(;
+        model = Wflow.ModelSection(;
+            type = Wflow.ModelType.sbm,
+            land_surface_water_reinfiltration__flag = false,
+        ),
+        input = Wflow.InputSection(;
+            path_forcing = "",
+            path_static = "",
+            basin__local_drain_direction = "",
+            river_location__mask = "",
+            subbasin_location__count = "",
+            forcing = Wflow.InputEntries(),
+            static = Wflow.InputEntries(),
+            location_maps = Wflow.PropertyDict(Dict{String, Any}()),
+        ),
+        path = "",
+    )
+
+    Wflow.update_bc_soil_model!(
+        soil_model,
+        atmospheric_forcing,
+        external_models,
+        domain,
+        config,
+    )
 
     @test soil_model.boundary_conditions.potential_transpiration[1] ≈ 3.456247174877943
     @test soil_model.boundary_conditions.potential_soilevaporation[1] ≈ 1.472182114066203
     @test soil_model.boundary_conditions.water_flux_surface[1] ≈ 0.02411574274509466
+    @test soil_model.boundary_conditions.potential_infiltration[1] ≈ 0.02411574274509466
+    @test soil_model.boundary_conditions.potential_infiltration_surfacewater[1] == 0.0
 end
 
 @testitem "unit: unsaturated_zone_flow!" begin

From cc9250eecb2f4c220bfda91391c962b33a786c18 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 18:10:02 +0200
Subject: [PATCH 21/24] max_reinfiltration_factor

---
 Wflow/src/soil/soil.jl                     | 5 ++++-
 docs/model_docs/land/landhydrology_sbm.qmd | 8 +++++---
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/Wflow/src/soil/soil.jl b/Wflow/src/soil/soil.jl
index c807e2882..effeec168 100644
--- a/Wflow/src/soil/soil.jl
+++ b/Wflow/src/soil/soil.jl
@@ -768,12 +768,15 @@ function update_available_for_infiltration!(
     (; waterdepth_land) = runoff.boundary_conditions
     (; river_fraction) = domain.land.parameters
 
+    # Fraction of overland flow that is made available for reinfiltration [-]
+    max_reinfiltration_fraction = 0.95
+
     n = length(potential_infiltration)
     threaded_foreach(1:n; basesize = 1000) do i
         potential_infiltration_surfacewater[i] = 0.0
         if do_surface_water_infiltration
             potential_infiltration_surfacewater[i] =
-                waterdepth_land[i] * (1.0 - river_fraction[i]) * 0.95
+                waterdepth_land[i] * (1.0 - river_fraction[i]) * max_reinfiltration_fraction
             water_flux_surface[i] += potential_infiltration_surfacewater[i]
         end
         potential_infiltration[i] = water_flux_surface[i]
diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index af943de14..209ce0022 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -972,9 +972,11 @@ land_surface_water_reinfiltration__flag = true
 
 This can be done with both `kinematic_wave` and `local_inertial` land routing settings. When
 enabling this process, the amount of overland flow is added as water that can potentially
-infiltrate. During infiltration, the amount of infiltrated water originating from overland
-flow is tracked using a simple ratio
-(`available_from_overland_flow/total_potential_infiltration`). After the infiltration
+infiltrate. To keep some water in the overland flow (and avoid draining a cell completely in
+one timestep), only 95% of the available overland flow depth (scaled by the non-river cell
+fraction `1 - river_fraction`) is made available for reinfiltration. During infiltration,
+the amount of infiltrated water originating from overland flow is tracked using a simple
+ratio (`available_from_overland_flow/total_potential_infiltration`). After the infiltration
 process, the volume of water in the overland flow is corrected to remove this water from the
 overland flow.
 

From df1770f5b198ae3b2dca956eb0f8767260269419 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 18:15:47 +0200
Subject: [PATCH 22/24] update docs

---
 docs/model_docs/land/landhydrology_sbm.qmd | 42 +++++++++++++++++-----
 1 file changed, 33 insertions(+), 9 deletions(-)

diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index 209ce0022..fc69e870f 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -970,15 +970,39 @@ following setting in the TOML file:
 land_surface_water_reinfiltration__flag = true
 ```
 
-This can be done with both `kinematic_wave` and `local_inertial` land routing settings. When
-enabling this process, the amount of overland flow is added as water that can potentially
-infiltrate. To keep some water in the overland flow (and avoid draining a cell completely in
-one timestep), only 95% of the available overland flow depth (scaled by the non-river cell
-fraction `1 - river_fraction`) is made available for reinfiltration. During infiltration,
-the amount of infiltrated water originating from overland flow is tracked using a simple
-ratio (`available_from_overland_flow/total_potential_infiltration`). After the infiltration
-process, the volume of water in the overland flow is corrected to remove this water from the
-overland flow.
+This can be done with both `kinematic_wave` and `local_inertial` land routing settings.
+When enabling this process, part of the overland flow depth is added to the water available
+for infiltration. To keep some water in the overland flow (and avoid draining a cell
+completely in one time step), the potential infiltration from surface water
+`potential_infiltration_surfacewater` $\SIb{}{mm\ t^{-1}}$ is computed from the overland
+flow depth `waterdepth_land` $\SIb{}{mm}$ as:
+
+$$
+    \subtext{i}{pot, sw} = \subtext{h}{land} (1 - \subtext{f}{river}) \subtext{f}{reinfilt, max},
+$$
+
+where $\SIb{\subtext{f}{river}}{-}$ is the river fraction of the cell and
+$\SIb{\subtext{f}{reinfilt, max}}{-}$ is the maximum reinfiltration fraction (fixed at
+$0.95$). This amount is added to the water flux at the soil surface `water_flux_surface`
+$\SIb{}{mm\ t^{-1}}$, so that the total water available for infiltration
+`potential_infiltration` $\SIb{}{mm\ t^{-1}}$ becomes:
+
+$$
+    \subtext{i}{pot} = \subtext{q}{surface} + \subtext{i}{pot, sw}.
+$$
+
+The regular infiltration routine then computes the actual infiltration `actinfilt`
+$\SIb{}{mm\ t^{-1}}$ and infiltration excess `infiltexcess` $\SIb{}{mm\ t^{-1}}$. After
+infiltration, the share originating from overland flow is tracked with the ratio
+`infilt_ratio` = `actinfilt` / `potential_infiltration`, giving the infiltration from
+surface water:
+
+$$
+    \subtext{i}{sw} = \subtext{i}{pot, sw} \cdot \frac{\subtext{i}{act}}{\subtext{i}{pot}}.
+$$
+
+Finally, `actinfilt`, `infiltexcess` and `water_flux_surface` are corrected by
+$\subtext{i}{pot, sw}$ to avoid double counting.
 
 ### Capillary rise
 

From 81c90885b6af211a798ac36918401473f2180342 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 18:18:51 +0200
Subject: [PATCH 23/24] clarify docs

---
 docs/model_docs/land/landhydrology_sbm.qmd | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/docs/model_docs/land/landhydrology_sbm.qmd b/docs/model_docs/land/landhydrology_sbm.qmd
index fc69e870f..b6e74ae80 100644
--- a/docs/model_docs/land/landhydrology_sbm.qmd
+++ b/docs/model_docs/land/landhydrology_sbm.qmd
@@ -1002,7 +1002,8 @@ $$
 $$
 
 Finally, `actinfilt`, `infiltexcess` and `water_flux_surface` are corrected by
-$\subtext{i}{pot, sw}$ to avoid double counting.
+$\subtext{i}{pot, sw}$ to avoid double counting of water already originating from surface
+water.
 
 ### Capillary rise
 

From 3ad8b1b5681df6342a92e12040d0d2ef9e8044c0 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 26 Aug 2026 18:29:40 +0200
Subject: [PATCH 24/24] infiltration_volume to boundary conditions

---
 .../routing/surface/surface_local_inertial.jl | 40 +++++++------------
 1 file changed, 14 insertions(+), 26 deletions(-)

diff --git a/Wflow/src/routing/surface/surface_local_inertial.jl b/Wflow/src/routing/surface/surface_local_inertial.jl
index 6124c1325..bc8be3f8f 100644
--- a/Wflow/src/routing/surface/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface/surface_local_inertial.jl
@@ -753,6 +753,10 @@ function update_bc_overland_flow_model!(
         net_runoff / 1000.0 * area / dt + net_runoff_river * area * 0.001 / dt
     overland_flow_model.boundary_conditions.runoff[river_indices] .+=
         get_flux_to_river(subsurface_flow, river_indices)
+
+    # infiltration volume [m³] from surface water, applied in `local_inertial_update_water_depth!`
+    @. overland_flow_model.boundary_conditions.infiltration_volume =
+        soil.variables.infilt_surfacewater * area * 0.001
     return nothing
 end
 
@@ -1020,11 +1024,7 @@ function local_inertial_update_water_depth!(
                     sum_at(river_v.q, edges_at_node.dst[inds_river[i]]) + land_v.qx[xd] -
                     land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] + land_bc.runoff[i] -
                     river_bc.abstraction[inds_river[i]]
-                ) * dt
-            # Apply surface water infiltration correction for river cells
-            if land_bc.infiltration_volume[i] > 0.0
-                land_v.storage[i] -= land_bc.infiltration_volume[i]
-            end
+                ) * dt - land_bc.infiltration_volume[i]
             if land_v.storage[i] < 0.0
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
@@ -1070,15 +1070,11 @@ function local_inertial_update_water_depth!(
                 (
                     land_v.qx[xd] - land_v.qx[i] + land_v.qy[yd] - land_v.qy[i] +
                     land_bc.runoff[i]
-                ) * dt
-            if land_bc.infiltration_volume[i] > 0.0
-                land_v.storage[i] -= land_bc.infiltration_volume[i]
-            end
+                ) * dt - land_bc.infiltration_volume[i]
             if land_v.storage[i] < 0.0
                 land_v.error[i] = land_v.error[i] + abs(land_v.storage[i])
                 land_v.storage[i] = 0.0 # set storage to zero
             end
-            # Apply surface water infiltration correction if available
             land_v.h[i] = land_v.storage[i] / (x_length[i] * y_length[i])
         end
     end
@@ -1364,20 +1360,12 @@ function FloodPlainModel(
 end
 
 """
-Update overland flow water level and storage for LocalInertialOverlandFlow model based on
-surface water infiltration. For local inertial flow, the infiltration is applied by updating
-the boundary conditions which are then used in the water depth update function.
+No-op for the `LocalInertialOverlandFlowModel`. Surface water infiltration is applied via
+the `infiltration_volume` boundary condition set in `update_bc_overland_flow_model!` and
+subtracted from land storage in `local_inertial_update_water_depth!`.
 """
-function update_overland_flow_and_depth!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
-    soil_model::SbmSoilModel,
-    domain::Domain,
-)
-    (; infilt_surfacewater) = soil_model.variables
-    (; area) = domain.land.parameters
-
-    # Update the boundary condition for surface water infiltration
-    # This will be used in local_inertial_update_water_depth!
-    overland_flow_model.boundary_conditions.infiltration_volume .=
-        infilt_surfacewater .* area .* 0.001
-end
+update_overland_flow_and_depth!(
+    ::LocalInertialOverlandFlowModel,
+    ::SbmSoilModel,
+    ::Domain,
+) = nothing