From df51bea1fda0719111d5acd8428af81e672b2710 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 3 Sep 2024 09:39:59 +0200
Subject: [PATCH 01/25] Update sbm.jl

---
 src/sbm.jl | 34 +++++++++++++++++++++++++++++-----
 1 file changed, 29 insertions(+), 5 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 60bfc6f24..e1654da66 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -144,6 +144,8 @@
     infiltsoilpath::Vector{T}
     # Infiltration excess water [mm Δt⁻¹]
     infiltexcess::Vector{T}
+    # Infiltration from surface water [mm Δt⁻¹]
+    infilt_surfacewater::Vector{T}
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{T}
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -705,6 +707,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         ae_openw_l = fill(mv, n),
         ae_openw_r = fill(mv, n),
         avail_forinfilt = fill(mv, n),
+        infilt_surfacewater = fill(mv, n),
         actinfilt = fill(mv, n),
         actinfiltsoil = fill(mv, n),
         actinfiltpath = fill(mv, n),
@@ -859,6 +862,7 @@ function update_until_recharge(sbm::SBM, config)
     transfermethod = get(config.model, "transfermethod", false)::Bool
     ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code?
     ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
+    do_surface_water_infiltration = get(config.model, "surface_water_infiltration", false)::Bool
 
     threaded_foreach(1:sbm.n, basesize = 250) do i
         if modelsnow
@@ -906,6 +910,14 @@ function update_until_recharge(sbm::SBM, config)
             sbm.waterfrac[i] * sbm.potential_evaporation[i],
         )
 
+        # Update land waterlevel
+        waterlevel_land = sbm.waterlevel_land[i] - ae_openw_l
+
+        if do_surface_water_infiltration
+            # Add land waterlevel to infiltration
+            avail_forinfilt += waterlevel_land
+        end
+
         # evap available for soil evaporation
         soilevap_fraction = max(
             sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.waterfrac[i] -
@@ -1111,7 +1123,18 @@ function update_until_recharge(sbm::SBM, config)
         end
 
         actinfilt = infiltsoilpath - du
-        excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du
+
+        # Scale infiltration from surface water based on the ratio between actinfil and
+        # infiltsoilpath (and prevent division by zero)
+        if do_surface_water_infiltration
+            infilt_ratio = iszero(infiltsoilpath) ? 0.0 : actinfilt / infiltsoilpath
+            infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
+            # Subtract waterlevel_land from this, as this water is already excess water
+            excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
+        else
+            infilt_surfacewater = 0.0
+            excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du
+        end
 
         # Separation between compacted and non compacted areas (correction with the satflow du)
         # This is required for D-Emission/Delwaq
@@ -1180,6 +1203,7 @@ function update_until_recharge(sbm::SBM, config)
         sbm.avail_forinfilt[i] = avail_forinfilt
         sbm.actinfilt[i] = actinfilt
         sbm.infiltexcess[i] = infiltexcess
+        sbm.infilt_surfacewater[i] = infilt_surfacewater
         sbm.recharge[i] = recharge
         sbm.transpiration[i] = transpiration
         sbm.soilevap[i] = soilevap
@@ -1295,7 +1319,7 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
         sbm.exfiltsatwater[i] = exfiltsatwater[i]
         sbm.exfiltustore[i] = exfiltustore
         sbm.runoff[i] = runoff
-        sbm.net_runoff[i] = runoff - sbm.ae_openw_l[i]
+        sbm.net_runoff[i] = runoff - sbm.ae_openw_l[i] - sbm.infilt_surfacewater[i]
         sbm.vwc[i] = vwc
         sbm.vwc_perc[i] = vwc_perc
         sbm.rootstore[i] = rootstore
@@ -1422,13 +1446,13 @@ function update_water_demand(sbm::SBM)
                         sbm.nonpaddy.maximum_irrigation_rate[i]
                     if depletion >= raw # start irrigation
                         irri_dem_gross += depletion
-                        # add depletion to irrigation gross demand when the maximum irrigation rate has been 
+                        # add depletion to irrigation gross demand when the maximum irrigation rate has been
                         # applied at the previous time step (to get volumetric water content at field capacity)
                     elseif depletion > 0.0 && max_irri_rate_applied # continue irrigation
                         irri_dem_gross += depletion
                     end
                 end
-                # limit irrigation demand to infiltration capacity    
+                # limit irrigation demand to infiltration capacity
                 infiltration_capacity =
                     sbm.soilinfredu[i] * (1.0 - sbm.pathfrac[i]) * sbm.infiltcapsoil[i]
                 irri_dem_gross = min(irri_dem_gross, infiltration_capacity)
@@ -1459,7 +1483,7 @@ function update_water_demand(sbm::SBM)
             end
             sbm.paddy.demand_gross[i] = irri_dem_gross
         end
-        # update gross water demands 
+        # update gross water demands
         sbm.allocation.irri_demand_gross[i] = irri_dem_gross
         sbm.allocation.nonirri_demand_gross[i] = industry_dem + domestic_dem + livestock_dem
         sbm.allocation.total_gross_demand[i] =

From ebe92bdd5c9dcce814970212b029df70ea6ed405 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Mon, 9 Sep 2024 09:59:20 +0200
Subject: [PATCH 02/25] add support for h_thres for overland flow

---
 src/flow.jl | 67 ++++++++++++++++++++++++++++++++++++++++++-----------
 src/sbm.jl  |  6 +++--
 2 files changed, 58 insertions(+), 15 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index c778b9f6c..5e5a24adc 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -63,6 +63,8 @@ end
     cel::Vector{T} | "m s-1"                        # Celerity of the kinematic wave
     to_river::Vector{T} | "m3 s-1"                  # Part of overland flow [m³ s⁻¹] that flows to the river
     kinwave_it::Bool | "-" | 0 | "none" | "none"    # Boolean for iterations kinematic wave
+    h_thresh::Vector{T} | "m"                       # Threshold for water level before flow occurs
+    pond_height::Vector{T} | "m"                    # Waterlevel of the pond (water stored on land surface)
 end
 
 function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, tstep, dt)
@@ -75,6 +77,9 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
         ncread(nc, config, "lateral.land.n"; sel = inds, defaults = 0.072, type = Float)
     n = length(inds)
 
+    h_thresh =
+        ncread(nc, config, "lateral.land.h_thresh"; sel = inds, defaults = 0, type = Float)
+
     sf_land = SurfaceFlowLand(
         beta = Float(0.6),
         sl = sl,
@@ -97,6 +102,8 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
         cel = zeros(Float, n),
         to_river = zeros(Float, n),
         kinwave_it = iterate,
+        h_thresh = h_thresh,
+        pond_height = zeros(Float, n),
     )
 
     return sf_land
@@ -199,6 +206,12 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     @. sf.alpha = sf.alpha_term * pow(sf.width, sf.alpha_pow)
     @. sf.qlat = sf.inwater / sf.dl
 
+    # Convert threshold to volume
+    threshold_volume = sf.h_thresh .* sf.width .* sf.dl
+
+    # Get pond volume at the start of the time step
+    pond_volume_current = sf.pond_height .* sf.width .* sf.dl
+
     sf.q_av .= 0.0
     sf.h_av .= 0.0
     sf.to_river .= 0.0
@@ -227,23 +240,51 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         )
                     end
 
-                    sf.q[v] = kinematic_wave(
-                        sf.qin[v],
-                        sf.q[v],
-                        sf.qlat[v],
-                        sf.alpha[v],
-                        sf.beta,
-                        dt,
-                        sf.dl[v],
+                    # Check potential pond volume, to see if flow is allowed to occur
+                    pond_volume_pot = max(
+                        pond_volume_current[v] +
+                        (sf.qlat[v] * sf.dl[v] * dt) +
+                        (sf.qin[v] * dt),
+                        0.0,
                     )
 
-                    # update h, only if surface width > 0.0
-                    if sf.width[v] > 0.0
-                        crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
-                        sf.h[v] = crossarea / sf.width[v]
+                    # Start kinematic wave if pond volume exceeds threshold
+                    if pond_volume_pot >= threshold_volume[v]
+                        q_ = kinematic_wave(
+                            sf.qin[v],
+                            sf.q[v],
+                            sf.qlat[v],
+                            sf.alpha[v],
+                            sf.beta,
+                            dt,
+                            sf.dl[v],
+                        )
+
+                        # update h, only if surface width > 0.0
+                        if sf.width[v] > 0.0
+                            crossarea = sf.alpha[v] * pow(q_, sf.beta)
+                            h_ = crossarea / sf.width[v]
+                        else
+                            h_ = 0.0
+                        end
+
+                    else
+                        # No flow if pond volume is below threshold
+                        q_ = 0.0
+                        h_ = 0.0
+                        # Update pond volume
+                        pond_volume_current[v] = pond_volume_pot
                     end
+
+                    # Update values
+                    sf.q[v] = q_
+                    sf.h[v] = h_
+                    sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v])
+
+                    # Update average values
                     sf.q_av[v] += sf.q[v]
-                    sf.h_av[v] += sf.h[v]
+                    sf.h_av[v] += (sf.h[v] + sf.pond_height[v])
+
                 end
             end
         end
diff --git a/src/sbm.jl b/src/sbm.jl
index eea77ac0f..b278d7ea3 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -872,7 +872,8 @@ function update_until_recharge(sbm::SBM, config)
     transfermethod = get(config.model, "transfermethod", false)::Bool
     ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code?
     ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
-    do_surface_water_infiltration = get(config.model, "surface_water_infiltration", false)::Bool
+    do_surface_water_infiltration =
+        get(config.model, "surface_water_infiltration", false)::Bool
 
     threaded_foreach(1:sbm.n, basesize = 250) do i
         if modelsnow
@@ -1140,7 +1141,8 @@ function update_until_recharge(sbm::SBM, config)
             infilt_ratio = iszero(infiltsoilpath) ? 0.0 : actinfilt / infiltsoilpath
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
             # Subtract waterlevel_land from this, as this water is already excess water
-            excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
+            excesswater =
+                avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
         else
             infilt_surfacewater = 0.0
             excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du

From 6bd036c20dddf0ff3d28d0916912e04bb1cb6cab Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Mon, 9 Sep 2024 10:22:06 +0200
Subject: [PATCH 03/25] update test and changelog

---
 docs/src/changelog.md | 12 +++++++++++-
 test/bmi.jl           |  4 ++--
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index de9f4cffb..d6c53462e 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -5,6 +5,16 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## Unreleased
+
+### Added
+- Support for reinfiltration of surface overland flow water. Can be set through the
+  `surface_water_infiltration` key in the `[model]` section. This way, the water stored on the
+  land surface is added to the water that is allowed to infiltrate into the soil.
+- Support for adding a water level flow-threshold for kinematic wave overland flow. This
+  setting can be changed by changing the `lateral.land.h_thresh` key in the settings file. This
+  can be added as a map or as a fixed value (through the use of `h_thresh.value`).
+
 ## v0.8.1 - 2024-08-27
 
 ### Fixed
@@ -35,7 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   is fixed by using `divrem` for the computation of the number of `steps` in this function.
   An error is thrown when the absolute remainder of `divrem` is larger than `eps()`, or when
   the number of `steps` is negative.
- - Fixed internal and external broken links in docs. 
+ - Fixed internal and external broken links in docs.
  - The internal time step of the local inertial model (`stable_timestep` function) can get
   zero when `LoopVectorization` is applied (`@tturbo`) to the for loop of these functions.
   This issue occured on a virtual machine, Windows 10 Enterprise, with Intel(R) Xeon(R) Gold
diff --git a/test/bmi.jl b/test/bmi.jl
index ba0f3751b..dc78246a1 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -20,8 +20,8 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "model information functions" begin
             @test BMI.get_component_name(model) == "sbm"
-            @test BMI.get_input_item_count(model) == 207
-            @test BMI.get_output_item_count(model) == 207
+            @test BMI.get_input_item_count(model) == 210
+            @test BMI.get_output_item_count(model) == 210
             to_check = [
                 "vertical.nlayers",
                 "vertical.theta_r",

From de4c2293a8246acf1fc15cc3468c79f14304ad86 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 17 Sep 2024 13:26:27 +0200
Subject: [PATCH 04/25] clean up overland flow threshold

---
 src/flow.jl | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 5e5a24adc..32e3588a9 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -250,7 +250,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
 
                     # Start kinematic wave if pond volume exceeds threshold
                     if pond_volume_pot >= threshold_volume[v]
-                        q_ = kinematic_wave(
+                        sf.q[v] = kinematic_wave(
                             sf.qin[v],
                             sf.q[v],
                             sf.qlat[v],
@@ -262,23 +262,19 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
 
                         # update h, only if surface width > 0.0
                         if sf.width[v] > 0.0
-                            crossarea = sf.alpha[v] * pow(q_, sf.beta)
-                            h_ = crossarea / sf.width[v]
-                        else
-                            h_ = 0.0
+                            crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
+                            sf.h[v] = crossarea / sf.width[v]
                         end
 
                     else
                         # No flow if pond volume is below threshold
-                        q_ = 0.0
-                        h_ = 0.0
+                        sf.q[v] = 0.0
+                        sf.h[v] = 0.0
                         # Update pond volume
                         pond_volume_current[v] = pond_volume_pot
                     end
 
                     # Update values
-                    sf.q[v] = q_
-                    sf.h[v] = h_
                     sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v])
 
                     # Update average values

From 6b70f1b3524f6ef9c194ea9f955ab312d56ffd9a Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 17 Sep 2024 13:26:55 +0200
Subject: [PATCH 05/25] clean code related to `excesswater`

---
 src/sbm.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index b278d7ea3..a4f7348ea 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1140,13 +1140,11 @@ function update_until_recharge(sbm::SBM, config)
         if do_surface_water_infiltration
             infilt_ratio = iszero(infiltsoilpath) ? 0.0 : actinfilt / infiltsoilpath
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
-            # Subtract waterlevel_land from this, as this water is already excess water
-            excesswater =
-                avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
         else
             infilt_surfacewater = 0.0
-            excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du
         end
+        # Subtract waterlevel_land from this, as this water is already excess water
+        excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
 
         # Separation between compacted and non compacted areas (correction with the satflow du)
         # This is required for D-Emission/Delwaq

From 68aed28baccc6fa2ca84de47c45b0d4ba62a52a4 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 17 Sep 2024 13:29:40 +0200
Subject: [PATCH 06/25] improve infilt_ratio

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

diff --git a/src/sbm.jl b/src/sbm.jl
index a4f7348ea..8cd2c72c8 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1138,7 +1138,7 @@ function update_until_recharge(sbm::SBM, config)
         # Scale infiltration from surface water based on the ratio between actinfil and
         # infiltsoilpath (and prevent division by zero)
         if do_surface_water_infiltration
-            infilt_ratio = iszero(infiltsoilpath) ? 0.0 : actinfilt / infiltsoilpath
+            infilt_ratio = iszero(actinfilt) ? 0.0 : actinfilt / avail_forinfilt
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
         else
             infilt_surfacewater = 0.0

From 127b68b443b2df38139a4d84c70d841a3689af60 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 17 Sep 2024 15:31:55 +0200
Subject: [PATCH 07/25] fix `excesswater`

---
 src/sbm.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 8cd2c72c8..86f12e7ce 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1140,11 +1140,13 @@ function update_until_recharge(sbm::SBM, config)
         if do_surface_water_infiltration
             infilt_ratio = iszero(actinfilt) ? 0.0 : actinfilt / avail_forinfilt
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
+            # Subtract waterlevel_land from this, as this water is already excess water
+            excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
         else
             infilt_surfacewater = 0.0
+            # Calculate excess water
+            excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du
         end
-        # Subtract waterlevel_land from this, as this water is already excess water
-        excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
 
         # Separation between compacted and non compacted areas (correction with the satflow du)
         # This is required for D-Emission/Delwaq

From 7f2d86c5b29ac9c8e78dfe7eb3d8412458eabbd1 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 25 Sep 2024 10:12:27 +0200
Subject: [PATCH 08/25] fix water balance issue; improve readability of
 excesswater computation

---
 src/sbm.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 86f12e7ce..4ced00e82 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1138,14 +1138,16 @@ function update_until_recharge(sbm::SBM, config)
         # Scale infiltration from surface water based on the ratio between actinfil and
         # infiltsoilpath (and prevent division by zero)
         if do_surface_water_infiltration
-            infilt_ratio = iszero(actinfilt) ? 0.0 : actinfilt / avail_forinfilt
+            infilt_ratio = iszero(avail_forinfilt) ? 0.0 : actinfilt / avail_forinfilt
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
-            # Subtract waterlevel_land from this, as this water is already excess water
-            excesswater = avail_forinfilt - waterlevel_land - infiltsoilpath - infiltexcess + du
+            # correct avail_forinfilt and actinfilt for the parts coming from surface water
+            # as this is already considered to be excess water
+            excesswater = (avail_forinfilt - waterlevel_land) - (actinfilt - infilt_surfacewater) - infiltexcess
+
         else
             infilt_surfacewater = 0.0
             # Calculate excess water
-            excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du
+            excesswater = avail_forinfilt - actinfilt - infiltexcess
         end
 
         # Separation between compacted and non compacted areas (correction with the satflow du)

From 6e1fd92204ebe8a19b34392fd58f0c95d1849e9f Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Sep 2024 14:28:29 +0200
Subject: [PATCH 09/25] move correction to runoff calculation

---
 src/sbm.jl       | 33 ++++++++++++++++++++++-----------
 src/sbm_model.jl |  1 +
 2 files changed, 23 insertions(+), 11 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 4ced00e82..8309d01f4 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1140,16 +1140,12 @@ function update_until_recharge(sbm::SBM, config)
         if do_surface_water_infiltration
             infilt_ratio = iszero(avail_forinfilt) ? 0.0 : actinfilt / avail_forinfilt
             infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
-            # correct avail_forinfilt and actinfilt for the parts coming from surface water
-            # as this is already considered to be excess water
-            excesswater = (avail_forinfilt - waterlevel_land) - (actinfilt - infilt_surfacewater) - infiltexcess
-
         else
             infilt_surfacewater = 0.0
-            # Calculate excess water
-            excesswater = avail_forinfilt - actinfilt - infiltexcess
         end
 
+        excesswater = avail_forinfilt - actinfilt - infiltexcess
+
         # Separation between compacted and non compacted areas (correction with the satflow du)
         # This is required for D-Emission/Delwaq
         if infiltsoil + infiltpath > 0.0
@@ -1251,7 +1247,7 @@ function update_until_recharge(sbm::SBM, config)
     end
 end
 
-function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
+function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater, config)
 
     threaded_foreach(1:sbm.n, basesize = 1000) do i
         usl, n_usl = set_layerthickness(zi[i], sbm.sumlayers[i], sbm.act_thickl[i])
@@ -1272,24 +1268,39 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
 
         ustoredepth = sum(@view usld[1:n_usl])
 
+        # Correct the excess water to only include the water that is not (yet) surface water,
+        # otherwise this water would be accounted for twice. This is only relevant if
+        # do_surface_water_infiltration is true
+        do_surface_water_infiltration =
+            get(config.model, "surface_water_infiltration", false)::Bool
+        if do_surface_water_infiltration
+            contribution_surfacewater = sbm.waterlevel_land[i] - sbm.ae_openw_l[i]
+            correction_surfacewater =
+                iszero(sbm.avail_forinfilt[i]) ? 1.0 :
+                1.0 - (contribution_surfacewater / sbm.avail_forinfilt[i])
+        else
+            correction_surfacewater = 1.0
+        end
+
         if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i]
             paddy_h_add =
                 exfiltustore +
                 exfiltsatwater[i] +
-                sbm.excesswater[i] +
+                sbm.excesswater[i] * correction_surfacewater +
                 sbm.runoff_land[i] +
-                sbm.infiltexcess[i]
+                sbm.infiltexcess[i] * correction_surfacewater
             runoff = max(paddy_h_add - sbm.paddy.h_max[i], 0.0)
             sbm.paddy.h[i] = paddy_h_add - runoff
         else
             runoff =
                 exfiltustore +
                 exfiltsatwater[i] +
-                sbm.excesswater[i] +
+                sbm.excesswater[i] * correction_surfacewater +
                 sbm.runoff_land[i] +
-                sbm.infiltexcess[i]
+                sbm.infiltexcess[i] * correction_surfacewater
         end
 
+
         # volumetric water content per soil layer and root zone
         vwc = sbm.vwc[i]
         vwc_perc = sbm.vwc_perc[i]
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index ec1032894..c706fa03c 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -518,6 +518,7 @@ function update_after_subsurfaceflow(
         vertical,
         lateral.subsurface.zi * 1000.0,
         lateral.subsurface.exfiltwater * 1000.0,
+        config,
     )
 
     ssf_toriver = lateral.subsurface.to_river ./ tosecond(basetimestep)

From 79ed50be65a6bb5dcd3afb241ec83cce21dd430e Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Sep 2024 14:28:54 +0200
Subject: [PATCH 10/25] add flow threshold for overland flow

---
 src/flow.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 32e3588a9..0fa2028c3 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -1091,9 +1091,17 @@ function initialize_shallowwater_land(
     froude_limit = get(config.model, "froude_limit", true)::Bool # limit flow to subcritical according to Froude number
     alpha = get(config.model, "inertial_flow_alpha", 0.7)::Float64 # stability coefficient for model time step (0.2-0.7)
     theta = get(config.model, "inertial_flow_theta", 0.8)::Float64 # weighting factor
-    h_thresh = get(config.model, "h_thresh", 1.0e-03)::Float64 # depth threshold for flow at link
+    # depth threshold for flow at link
+    h_thresh = ncread(
+        nc,
+        config,
+        "lateral.land.h_thresh";
+        sel = inds,
+        defaults = 1.0e-03,
+        type = Float,
+    )
 
-    @info "Local inertial approach is used for overlandflow." alpha theta h_thresh froude_limit
+    @info "Local inertial approach is used for overlandflow." alpha theta froude_limit
 
     n_land =
         ncread(nc, config, "lateral.land.n"; sel = inds, defaults = 0.072, type = Float)
@@ -1294,7 +1302,7 @@ function shallowwater_update(
             zs_max = max(zs_x, zs_xu)
             hf = (zs_max - sw.zx_max[i])
 
-            if hf > sw.h_thresh
+            if hf > sw.h_thresh[i]
                 length = T(0.5) * (sw.xl[i] + sw.xl[xu]) # can be precalculated
                 sw.qx[i] = local_inertial_flow(
                     sw.theta,
@@ -1334,7 +1342,7 @@ function shallowwater_update(
             zs_max = max(zs_y, zs_yu)
             hf = (zs_max - sw.zy_max[i])
 
-            if hf > sw.h_thresh
+            if hf > sw.h_thresh[i]
                 length = T(0.5) * (sw.yl[i] + sw.yl[yu]) # can be precalculated
                 sw.qy[i] = local_inertial_flow(
                     sw.theta,

From de29798d98f5688ff78989233f969ed8089ba32b Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Sep 2024 16:12:58 +0200
Subject: [PATCH 11/25] fix `h_thresh` type in struct definition

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

diff --git a/src/flow.jl b/src/flow.jl
index 0fa2028c3..d9d8f5764 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -1052,7 +1052,7 @@ const dirs = (:yd, :xd, :xu, :yu)
     g::T | "m2 s-1" | 0 | "scalar"                          # acceleration due to gravity
     theta::T | "-" | 0 | "scalar"                           # weighting factor (de Almeida et al., 2012)
     alpha::T | "-" | 0 | "scalar"                           # stability coefficient (de Almeida et al., 2012)
-    h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
+    h_thresh::Vector{T} | "m" | 0 | "scalar"                # depth threshold for calculating flow
     dt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     qy0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in y direction at previous time step
     qx0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in x direction at previous time step

From 2cd95774aac28b1a5ab27686493d0abe020f88b6 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 26 Sep 2024 16:16:26 +0200
Subject: [PATCH 12/25] update changelog

---
 docs/src/changelog.md | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index d6c53462e..e9560baab 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -11,9 +11,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Support for reinfiltration of surface overland flow water. Can be set through the
   `surface_water_infiltration` key in the `[model]` section. This way, the water stored on the
   land surface is added to the water that is allowed to infiltrate into the soil.
-- Support for adding a water level flow-threshold for kinematic wave overland flow. This
-  setting can be changed by changing the `lateral.land.h_thresh` key in the settings file. This
-  can be added as a map or as a fixed value (through the use of `h_thresh.value`).
+- Support for adding a water level flow-threshold for both kinematic wave and local inertial
+  overland flow. This setting can be changed by changing the `lateral.land.h_thresh` key in the
+  settings file. This can be added as a map or as a fixed value (through the use of
+  `h_thresh.value`).
 
 ## v0.8.1 - 2024-08-27
 

From f6bfe8a6dd51798c01779b260ef14a4c470867c2 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 09:23:57 +0200
Subject: [PATCH 13/25] update docs, fix gwf error

---
 docs/src/model_docs/params_lateral.md     |  2 ++
 docs/src/model_docs/params_vertical.md    |  1 +
 docs/src/user_guide/additional_options.md | 28 +++++++++++++++++++++++
 src/sbm_gwf_model.jl                      |  1 +
 4 files changed, 32 insertions(+)

diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 6456500e8..4c58258f0 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -74,6 +74,8 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `cel`           | celerity of kinematic wave | m s``^{-1}`` | - |
 | `to_river`      | part of overland flow that flows to the river | m``^3`` s``^{-1}`` | - |
 | `kinwave_it`    | boolean for kinematic wave iterations | - | false |
+| **`h_thresh`**  | threshold for water level before flow occurs | m | 0 |
+| `pond_height`   | waterlevel of the pond (water stored on land surface) | m | - |
 
 ### [Reservoirs](@id reservoir_params)
 The Table below shows the parameters (fields) of struct `SimpleReservoir`, including a
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 12e52bcd6..e52431555 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -117,6 +117,7 @@ profile `kv` is used and `z_layered` is required as input.
 | `actinfiltpath` | actual infiltration into compacted fraction  | mm Δt``^{-1}`` | - |
 | `infiltsoilpath` | infiltration into the unsaturated zone | mm Δt``^{-1}`` | - |
 | `infiltexcess` | infiltration excess water | mm Δt``^{-1}`` | - |
+| `infilt_surfacewater` | part of infiltration that originates from surface water | mm Δt``^{-1}`` | - |
 | `excesswater` | water that cannot infiltrate due to saturated soil (saturation excess) | mm Δt``^{-1}`` | - |
 | `exfiltsatwater` | water exfiltrating during saturation excess conditions | mm Δt``^{-1}`` | - |
 | `exfiltustore` | water exfiltrating from unsaturated store because of change in water table | mm Δt``^{-1}`` | - |
diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md
index 6b5baa03b..2d01eef05 100644
--- a/docs/src/user_guide/additional_options.md
+++ b/docs/src/user_guide/additional_options.md
@@ -183,6 +183,34 @@ irrigation_trigger = "irrigation_trigger"
 h = "h_paddy"
 ```
 
+## Enabling re-infiltration of overland flow
+To allow for re-infiltration of overland flow, the following model flag needs to be set:
+
+```toml
+[model]
+surface_water_infiltration = true
+```
+
+When using this option, the average surface water level (of the previous time step) is added to
+the volume of water that is allowed to infiltrate (`avail_forinfilt`). To keep track of the
+infiltrated water that originates from the surface water, the variable `infilt_surfacewater` is
+used. This variable is then subtracted from the net runoff, such that the overland flow is
+corrected for the loss of this water. Lastly, some corrections are applied on the
+`excesswater`, to prevent the surface water from being counted twice as excess water. This
+option works for both kinematic wave and local inertial overland flow.
+
+## Enabling a water level threshold for overland flow
+A water level threshold can be set for the overland flow. As long as this threshold is not
+exceeded, flow is not allowed to occur, and the water will be considered a "pond" (the variable
+`lateral.land.pond_height` can be used to keep track of its water level for kinematic wave
+overland flow, and the variable `lateral.land.h` can be used). When this threshold is exceeded,
+flow is allowed to occur. This flow threshold can be set as a map from the staticmaps:
+
+```toml
+[input.lateral.land]
+h_thresh = "overland_flow_threshold"
+```
+
 ## [Using multithreading] (@id multi_threading)
 
 ### Using wflow in Julia
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 2dfe4a2ae..f313fc745 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -600,6 +600,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
         vertical,
         (network.land.altitude .- aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
         exfiltwater .* 1000.0,
+        config,
     )
 
     ssf_toriver = zeros(vertical.n)

From ae37c4a6a851b2b1de9754cdbe1808c3b0c3f649 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 11:56:37 +0200
Subject: [PATCH 14/25] add tests

---
 test/run_sbm.jl | 78 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)

diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index d50ccb752..e77e56a8a 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -540,3 +540,81 @@ Wflow.close_files(model, delete_output = false)
 
     Wflow.close_files(model, delete_output = false)
 end
+
+
+# test infiltration of surface water
+@testset "surface water infiltration" begin
+    i_surface = 498 # pixel with surface water, infiltration and excesswater
+    i_nosurface = 2444 # pixel without surface water
+
+    tomlpath = joinpath(@__DIR__, "sbm_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.model.surface_water_infiltration = true
+
+    model = Wflow.initialize_sbm_model(config)
+
+    model = Wflow.run_timestep(model)
+    model = Wflow.run_timestep(model)
+
+    @unpack vertical = model
+
+    # Test is surface water is indeed added to avail_forinfilt
+    @test vertical.waterlevel_land[i_surface] ≈ 16.191219691562683
+    @test vertical.waterlevel_land[i_nosurface] == 0.0
+
+    @test vertical.avail_forinfilt[i_surface] ≈ 16.24809470532854
+    @test vertical.avail_forinfilt[i_nosurface] ≈ 0.08956446991756856
+
+    # Test if the paritioning of infiltrated water is correct
+    @test vertical.infilt_surfacewater[i_surface] ≈ 9.1388544728681
+    @test vertical.infilt_surfacewater[i_nosurface] == 0.0
+
+    # Test excesswater and runoff for two locations
+    @test vertical.excesswater[i_surface] ≈ 7.077138112678853
+    @test vertical.runoff[i_surface] ≈ 88.75818647205926
+
+    # Test is net_runoff is indeed lower than runoff
+    @test vertical.net_runoff[i_surface] ≈ 79.61933199919116
+    @test vertical.net_runoff[i_nosurface] == 0.0
+
+    Wflow.close_files(model, delete_output = false)
+end
+
+
+# test flow threshold for kinematic overland flow
+@testset "surface water infiltration" begin
+    idx = 498 # pixel with surface water, infiltration and excesswater
+    threshold =  0.02
+
+    tomlpath = joinpath(@__DIR__, "sbm_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.input.lateral.land.h_thresh = Dict("value" => threshold)
+
+    mod_threshold = Wflow.initialize_sbm_model(config)
+
+
+    # run two timesteps for both models
+    mod_threshold = Wflow.run_timestep(mod_threshold)
+    mod_threshold = Wflow.run_timestep(mod_threshold)
+    Wflow.close_files(mod_threshold, delete_output = false)
+
+
+    config = Wflow.Config(tomlpath)
+    mod_nothreshold = Wflow.initialize_sbm_model(config)
+
+    mod_nothreshold = Wflow.run_timestep(mod_nothreshold)
+    mod_nothreshold = Wflow.run_timestep(mod_nothreshold)
+    Wflow.close_files(mod_nothreshold, delete_output = false)
+
+
+    @unpack lateral = mod_threshold
+    lateral_threshold = lateral
+    @unpack lateral = mod_nothreshold
+    lateral_nothreshold = lateral
+
+    @test lateral_threshold.land.q_av[idx] ≈ 1.3364541870336923
+    @test lateral_threshold.land.q_av[idx] < lateral_nothreshold.land.q_av[idx]
+    @test lateral_threshold.land.pond_height[idx] ≈ 0.01591595571103741
+    @test lateral_nothreshold.land.pond_height[idx] == 0.0
+
+end
\ No newline at end of file

From 8560072f376735ef8210eb3ebc3a7b4ec65ad23b Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 12:00:26 +0200
Subject: [PATCH 15/25] update docs

---
 docs/src/user_guide/additional_options.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md
index 2d01eef05..3bfa29c6e 100644
--- a/docs/src/user_guide/additional_options.md
+++ b/docs/src/user_guide/additional_options.md
@@ -211,6 +211,8 @@ flow is allowed to occur. This flow threshold can be set as a map from the stati
 h_thresh = "overland_flow_threshold"
 ```
 
+It should be noted that (for kinematic wave overland flow) flow is allowed to occur when the `pond_height` plus the incoming fluxes would be greater than the `h_thresh`. As a result, flow can occur at pond heights that are slightly lower than the set threshold.
+
 ## [Using multithreading] (@id multi_threading)
 
 ### Using wflow in Julia

From e5067f756cd45ae104c4e5c5636fdb5fdc1011a6 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 12:01:49 +0200
Subject: [PATCH 16/25] update WflowServer tests

---
 server/test/client.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/server/test/client.jl b/server/test/client.jl
index 36e6b5dc5..f99c73dbf 100644
--- a/server/test/client.jl
+++ b/server/test/client.jl
@@ -32,8 +32,8 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 207)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 207)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 210)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 210)
     to_check = [
         "vertical.nlayers",
         "vertical.theta_r",

From db41548f577f284b99b8ca6353c1eda59bc2c9cc Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 13:00:17 +0200
Subject: [PATCH 17/25] wrap lines

---
 docs/src/user_guide/additional_options.md | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md
index 3bfa29c6e..60e59cace 100644
--- a/docs/src/user_guide/additional_options.md
+++ b/docs/src/user_guide/additional_options.md
@@ -211,7 +211,9 @@ flow is allowed to occur. This flow threshold can be set as a map from the stati
 h_thresh = "overland_flow_threshold"
 ```
 
-It should be noted that (for kinematic wave overland flow) flow is allowed to occur when the `pond_height` plus the incoming fluxes would be greater than the `h_thresh`. As a result, flow can occur at pond heights that are slightly lower than the set threshold.
+It should be noted that (for kinematic wave overland flow) flow is allowed to occur when the
+`pond_height` plus the incoming fluxes would be greater than the `h_thresh`. As a result, flow
+can occur at pond heights that are slightly lower than the set threshold.
 
 ## [Using multithreading] (@id multi_threading)
 

From b5d4889cdf47e6ece77d11f34994fd61d0771043 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 27 Sep 2024 13:01:03 +0200
Subject: [PATCH 18/25] Update run_sbm.jl

---
 test/run_sbm.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index e77e56a8a..f97cf9df2 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -617,4 +617,4 @@ end
     @test lateral_threshold.land.pond_height[idx] ≈ 0.01591595571103741
     @test lateral_nothreshold.land.pond_height[idx] == 0.0
 
-end
\ No newline at end of file
+end

From 56e7d326ada885f254ccad595bca4fb06c57f2a5 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 11 Oct 2024 17:10:12 +0200
Subject: [PATCH 19/25] fix computation of pond volume to cell area

---
 src/flow.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index d9d8f5764..64e13aece 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -196,7 +196,7 @@ end
 
 
 function update(sf::SurfaceFlowLand, network, frac_toriver)
-    @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
+    @unpack area, graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
         network
 
     ns = length(subdomain_order)
@@ -207,10 +207,10 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     @. sf.qlat = sf.inwater / sf.dl
 
     # Convert threshold to volume
-    threshold_volume = sf.h_thresh .* sf.width .* sf.dl
+    threshold_volume = sf.h_thresh .* area
 
     # Get pond volume at the start of the time step
-    pond_volume_current = sf.pond_height .* sf.width .* sf.dl
+    pond_volume_current = sf.pond_height .* area
 
     sf.q_av .= 0.0
     sf.h_av .= 0.0
@@ -275,7 +275,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                     end
 
                     # Update values
-                    sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v])
+                    sf.pond_height[v] = pond_volume_current[v] / area[v]
 
                     # Update average values
                     sf.q_av[v] += sf.q[v]

From 4a39e4ad45771ac5392cb09022e87f717d58210b Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 11 Oct 2024 17:10:41 +0200
Subject: [PATCH 20/25] fix such that only water above the threshold is allowed
 to flow

---
 src/flow.jl | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 64e13aece..80fe97495 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -240,20 +240,25 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         )
                     end
 
+                    # Calculate potential inflow
+                    pot_inflow = (sf.qlat[v] * sf.dl[v] * dt) + (sf.qin[v] * dt)
+
                     # Check potential pond volume, to see if flow is allowed to occur
                     pond_volume_pot = max(
-                        pond_volume_current[v] +
-                        (sf.qlat[v] * sf.dl[v] * dt) +
-                        (sf.qin[v] * dt),
+                        pond_volume_current[v] + pot_inflow,
                         0.0,
                     )
 
                     # Start kinematic wave if pond volume exceeds threshold
                     if pond_volume_pot >= threshold_volume[v]
+                        # Calculate fraction of pond volume that is above threshold
+                        flowing_fraction = iszero(pot_inflow) ? 1.0 :
+                            (pond_volume_pot - threshold_volume[v]) / pot_inflow
+
                         sf.q[v] = kinematic_wave(
-                            sf.qin[v],
+                            sf.qin[v] * flowing_fraction,
                             sf.q[v],
-                            sf.qlat[v],
+                            sf.qlat[v] * flowing_fraction,
                             sf.alpha[v],
                             sf.beta,
                             dt,
@@ -266,6 +271,8 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                             sf.h[v] = crossarea / sf.width[v]
                         end
 
+                        # Update pond volume
+                        pond_volume_current[v] = threshold_volume[v]
                     else
                         # No flow if pond volume is below threshold
                         sf.q[v] = 0.0

From 01bd103b4e71672b63a5416a271171b2e80ac702 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Wed, 16 Oct 2024 13:10:03 +0200
Subject: [PATCH 21/25] revert area calculation back to width*dl

---
 src/flow.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 80fe97495..0cdec5c79 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -196,7 +196,7 @@ end
 
 
 function update(sf::SurfaceFlowLand, network, frac_toriver)
-    @unpack area, graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
+    @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
         network
 
     ns = length(subdomain_order)
@@ -207,10 +207,10 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     @. sf.qlat = sf.inwater / sf.dl
 
     # Convert threshold to volume
-    threshold_volume = sf.h_thresh .* area
+    threshold_volume = sf.h_thresh .* sf.width .* sf.dl
 
     # Get pond volume at the start of the time step
-    pond_volume_current = sf.pond_height .* area
+    pond_volume_current = sf.pond_height .* sf.width .* sf.dl
 
     sf.q_av .= 0.0
     sf.h_av .= 0.0
@@ -282,7 +282,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                     end
 
                     # Update values
-                    sf.pond_height[v] = pond_volume_current[v] / area[v]
+                    sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v])
 
                     # Update average values
                     sf.q_av[v] += sf.q[v]

From 26f5004e9479baab156623d6c19a2dc42e88a206 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Fri, 18 Oct 2024 13:41:39 +0200
Subject: [PATCH 22/25] remove dependency of waterfrac; fix flowing_fraction

---
 src/flow.jl |  7 +++++--
 src/sbm.jl  | 38 +++++++++++++++++++++++---------------
 2 files changed, 28 insertions(+), 17 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 0cdec5c79..650bd8df3 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -252,8 +252,11 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                     # Start kinematic wave if pond volume exceeds threshold
                     if pond_volume_pot >= threshold_volume[v]
                         # Calculate fraction of pond volume that is above threshold
-                        flowing_fraction = iszero(pot_inflow) ? 1.0 :
-                            (pond_volume_pot - threshold_volume[v]) / pot_inflow
+                        if pot_inflow <= 0.0
+                            flowing_fraction = 1.0
+                        else
+                            flowing_fraction = (pond_volume_pot - threshold_volume[v]) / pot_inflow
+                        end
 
                         sf.q[v] = kinematic_wave(
                             sf.qin[v] * flowing_fraction,
diff --git a/src/sbm.jl b/src/sbm.jl
index 8309d01f4..432408e40 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -156,6 +156,8 @@
     infiltexcess::Vector{T}
     # Infiltration from surface water [mm Δt⁻¹]
     infilt_surfacewater::Vector{T}
+    # Contribution from surface water [mm Δt⁻¹]
+    contribution_surfacewater::Vector{T}
     # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹]
     excesswater::Vector{T}
     # Water exfiltrating during saturation excess conditions [mm Δt⁻¹]
@@ -718,6 +720,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         ae_openw_r = fill(mv, n),
         avail_forinfilt = fill(mv, n),
         infilt_surfacewater = fill(mv, n),
+        contribution_surfacewater = fill(mv, n),
         actinfilt = fill(mv, n),
         actinfiltsoil = fill(mv, n),
         actinfiltpath = fill(mv, n),
@@ -916,26 +919,31 @@ function update_until_recharge(sbm::SBM, config)
             sbm.waterlevel_river[i] * sbm.riverfrac[i],
             sbm.riverfrac[i] * sbm.potential_evaporation[i],
         )
-        ae_openw_l = min(
-            sbm.waterlevel_land[i] * sbm.waterfrac[i],
-            sbm.waterfrac[i] * sbm.potential_evaporation[i],
+        # Determine remaining fraction of available potential_evaporation
+        soilevap_fraction = max(
+            sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.glacierfrac[i],
+            0.0,
         )
+        # Determine amount of surface water
+        surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i])
 
-        # Update land waterlevel
-        waterlevel_land = sbm.waterlevel_land[i] - ae_openw_l
+        # Calculate open water evaporation
+        ae_openw_l = min(
+            surface_water,
+            sbm.potential_evaporation[i] * soilevap_fraction,
+        )
 
+        # Update land waterlevel and scale to part of cell not covered by rivers
         if do_surface_water_infiltration
+            contribution_surfacewater = surface_water - ae_openw_l
             # Add land waterlevel to infiltration
-            avail_forinfilt += waterlevel_land
+            avail_forinfilt += contribution_surfacewater
+        else
+            contribution_surfacewater = 0.0
         end
 
         # evap available for soil evaporation
-        soilevap_fraction = max(
-            sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.waterfrac[i] -
-            sbm.glacierfrac[i],
-            0.0,
-        )
-        potsoilevap = soilevap_fraction * sbm.potential_evaporation[i]
+        potsoilevap = (sbm.potential_evaporation[i] * soilevap_fraction) - ae_openw_l
 
         if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i]
             evap_paddy_water = min(sbm.paddy.h[i], potsoilevap)
@@ -1139,7 +1147,7 @@ function update_until_recharge(sbm::SBM, config)
         # infiltsoilpath (and prevent division by zero)
         if do_surface_water_infiltration
             infilt_ratio = iszero(avail_forinfilt) ? 0.0 : actinfilt / avail_forinfilt
-            infilt_surfacewater = max(0.0, waterlevel_land * infilt_ratio)
+            infilt_surfacewater = max(0.0, contribution_surfacewater * infilt_ratio)
         else
             infilt_surfacewater = 0.0
         end
@@ -1214,6 +1222,7 @@ function update_until_recharge(sbm::SBM, config)
         sbm.actinfilt[i] = actinfilt
         sbm.infiltexcess[i] = infiltexcess
         sbm.infilt_surfacewater[i] = infilt_surfacewater
+        sbm.contribution_surfacewater[i] = contribution_surfacewater
         sbm.recharge[i] = recharge
         sbm.transpiration[i] = transpiration
         sbm.soilevap[i] = soilevap
@@ -1274,10 +1283,9 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater, config)
         do_surface_water_infiltration =
             get(config.model, "surface_water_infiltration", false)::Bool
         if do_surface_water_infiltration
-            contribution_surfacewater = sbm.waterlevel_land[i] - sbm.ae_openw_l[i]
             correction_surfacewater =
                 iszero(sbm.avail_forinfilt[i]) ? 1.0 :
-                1.0 - (contribution_surfacewater / sbm.avail_forinfilt[i])
+                1.0 - (sbm.contribution_surfacewater[i] / sbm.avail_forinfilt[i])
         else
             correction_surfacewater = 1.0
         end

From dd921b77cbabfb2c985212b3b83b780c264d2ac1 Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 24 Oct 2024 17:19:19 +0200
Subject: [PATCH 23/25] fix open water evaporation

---
 src/sbm.jl | 31 +++++++++++++++++++++++--------
 1 file changed, 23 insertions(+), 8 deletions(-)

diff --git a/src/sbm.jl b/src/sbm.jl
index 432408e40..3c32b9500 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -909,7 +909,18 @@ function update_until_recharge(sbm::SBM, config)
         ustoredepth = sum(@view sbm.ustorelayerdepth[i][1:sbm.nlayers[i]])
 
         runoff_river = min(1.0, sbm.riverfrac[i]) * avail_forinfilt
-        runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
+
+        # Calculate the initial capacity of the unsaturated store
+        ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
+        # Add a flag if the soil is saturated
+        check_saturated = ustorecapacity < 0.1
+
+        if check_saturated
+            runoff_land = min(1.0, 1.0 - sbm.riverfrac[i]) * avail_forinfilt
+        else
+            runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt
+        end
+
         if !isnothing(sbm.paddy) || !isnothing(sbm.nonpaddy)
             avail_forinfilt = avail_forinfilt + sbm.allocation.irri_alloc[i]
         end
@@ -928,10 +939,17 @@ function update_until_recharge(sbm::SBM, config)
         surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i])
 
         # Calculate open water evaporation
-        ae_openw_l = min(
-            surface_water,
-            sbm.potential_evaporation[i] * soilevap_fraction,
-        )
+        if check_saturated
+            ae_openw_l = min(
+                surface_water,
+                sbm.potential_evaporation[i] * soilevap_fraction,
+            )
+        else
+            ae_openw_l = min(
+                surface_water,
+                sbm.potential_evaporation[i] * sbm.waterfrac[i],
+            )
+        end
 
         # Update land waterlevel and scale to part of cell not covered by rivers
         if do_surface_water_infiltration
@@ -954,9 +972,6 @@ function update_until_recharge(sbm::SBM, config)
             evap_paddy_water = 0.0
         end
 
-        # Calculate the initial capacity of the unsaturated store
-        ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
-
         # Calculate the infiltration flux into the soil column
         infiltsoilpath,
         infiltsoil,

From d05f673db91e86e055925e6d1750837ebef85eca Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Thu, 24 Oct 2024 17:19:35 +0200
Subject: [PATCH 24/25] move pond calculations to fully inside for loop

---
 src/flow.jl | 23 +++++++++++------------
 1 file changed, 11 insertions(+), 12 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 650bd8df3..594c32c5a 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -206,12 +206,6 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     @. sf.alpha = sf.alpha_term * pow(sf.width, sf.alpha_pow)
     @. sf.qlat = sf.inwater / sf.dl
 
-    # Convert threshold to volume
-    threshold_volume = sf.h_thresh .* sf.width .* sf.dl
-
-    # Get pond volume at the start of the time step
-    pond_volume_current = sf.pond_height .* sf.width .* sf.dl
-
     sf.q_av .= 0.0
     sf.h_av .= 0.0
     sf.to_river .= 0.0
@@ -240,22 +234,26 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         )
                     end
 
+                    # Convert threshold to volume
+                    threshold_volume = sf.h_thresh[v] * sf.width[v] * sf.dl[v]
+                    pond_volume_current = sf.pond_height[v] * sf.width[v] * sf.dl[v]
+
                     # Calculate potential inflow
                     pot_inflow = (sf.qlat[v] * sf.dl[v] * dt) + (sf.qin[v] * dt)
 
                     # Check potential pond volume, to see if flow is allowed to occur
                     pond_volume_pot = max(
-                        pond_volume_current[v] + pot_inflow,
+                        pond_volume_current + pot_inflow,
                         0.0,
                     )
 
                     # Start kinematic wave if pond volume exceeds threshold
-                    if pond_volume_pot >= threshold_volume[v]
+                    if pond_volume_pot >= threshold_volume
                         # Calculate fraction of pond volume that is above threshold
                         if pot_inflow <= 0.0
                             flowing_fraction = 1.0
                         else
-                            flowing_fraction = (pond_volume_pot - threshold_volume[v]) / pot_inflow
+                            flowing_fraction = (pond_volume_pot - threshold_volume) / pot_inflow
                         end
 
                         sf.q[v] = kinematic_wave(
@@ -275,17 +273,18 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         end
 
                         # Update pond volume
-                        pond_volume_current[v] = threshold_volume[v]
+                        pond_volume_current = threshold_volume
+
                     else
                         # No flow if pond volume is below threshold
                         sf.q[v] = 0.0
                         sf.h[v] = 0.0
                         # Update pond volume
-                        pond_volume_current[v] = pond_volume_pot
+                        pond_volume_current = pond_volume_pot
                     end
 
                     # Update values
-                    sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v])
+                    sf.pond_height[v] = pond_volume_current / (sf.width[v] * sf.dl[v])
 
                     # Update average values
                     sf.q_av[v] += sf.q[v]

From 2dfa845b968c853746a1ab8569eccacf5ebf181a Mon Sep 17 00:00:00 2001
From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>
Date: Tue, 29 Oct 2024 15:45:48 +0100
Subject: [PATCH 25/25] fix flowing_fraction and volume calculation

---
 src/flow.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/flow.jl b/src/flow.jl
index 594c32c5a..c30b1697b 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -254,6 +254,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                             flowing_fraction = 1.0
                         else
                             flowing_fraction = (pond_volume_pot - threshold_volume) / pot_inflow
+                            flowing_fraction = min(flowing_fraction, 1.0)
                         end
 
                         sf.q[v] = kinematic_wave(
@@ -297,7 +298,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     sf.q_av ./= its
     sf.h_av ./= its
     sf.to_river ./= its
-    sf.volume .= sf.dl .* sf.width .* sf.h
+    sf.volume .= sf.dl .* sf.width .* (sf.h .+ sf.pond_height)
 end
 
 function update(sf::SurfaceFlowRiver, network, doy)