diff --git a/server/test/client.jl b/server/test/client.jl
index c89ed699f..c5436931c 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" => 206)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 206)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 204)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 204)
     to_check = [
         "vertical.soil.parameters.nlayers",
         "vertical.soil.parameters.theta_r",
@@ -81,7 +81,7 @@ vwc_1_size = 0
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.198747900215207f0, 2.6880427720508515f0, 3.4848783702629564f0]
+          [2.1007361866518766, 2.5702292750107687, 3.2904803551115727]
     msg =
         (fn = "set_value", name = "vertical.soil.variables.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")
diff --git a/src/Wflow.jl b/src/Wflow.jl
index 84b5af460..91555a124 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -49,7 +49,7 @@ using Parameters: @with_kw
 using Polyester: @batch
 using ProgressLogging: @progress
 using StaticArrays: SVector, pushfirst, setindex
-using Statistics: mean, median, quantile!
+using Statistics: mean, median, quantile!, quantile
 using TerminalLoggers
 using TOML: TOML
 
diff --git a/src/flow.jl b/src/flow.jl
index f05fe2b7e..81df9c049 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -6,7 +6,12 @@
     volume::Vector{T} | "m3"                # Kinematic wave volume [m³] (based on water depth h)
     h::Vector{T} | "m"                      # Water depth [m]
     h_av::Vector{T} | "m"                   # Average water depth [m]
-    celerity::Vector{T} | "m s-1"           # Celerity of the kinematic wave
+end
+
+@with_kw struct TimeStepping{T}
+    dt::Union{Vector{T}, T}
+    adaptive::Bool = true
+    cfl::T = 0.70
 end
 
 function FlowVariables(n)
@@ -18,7 +23,6 @@ function FlowVariables(n)
         volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
-        celerity = zeros(Float, n),
     )
     return variables
 end
@@ -32,23 +36,19 @@ end
     alpha_pow::T                            # Used in the power part of alpha [-]
     alpha_term::Vector{T} | "-"             # Term used in computation of alpha [-]
     alpha::Vector{T} | "s3/5 m1/5"          # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation
-    its::Int                                # Number of fixed iterations [-]
-    kinwave_it::Bool                        # Boolean for iterations kinematic wave
 end
 
-function ManningFlowParameters(slope, mannings_n, dl, width, iterate, tstep, dt)
+function ManningFlowParameters(slope, mannings_n, dl, width)
     n = length(slope)
     parameters = ManningFlowParameters(;
         beta = Float(0.6),
         slope,
         mannings_n,
         dl = dl,
-        its = tstep > 0 ? Int(cld(dt, tstep)) : tstep,
         width,
         alpha_pow = Float((2.0 / 3.0) * 0.6),
         alpha_term = fill(mv, n),
         alpha = fill(mv, n),
-        kinwave_it = iterate,
     )
     return parameters
 end
@@ -68,7 +68,7 @@ function Base.getproperty(v::RiverFlowParameters, s::Symbol)
     end
 end
 
-function RiverFlowParameters(nc, config, inds, dl, width, iterate, tstep)
+function RiverFlowParameters(nc, config, inds, dl, width)
     n_river = ncread(
         nc,
         config,
@@ -102,9 +102,7 @@ function RiverFlowParameters(nc, config, inds, dl, width, iterate, tstep)
     )
     clamp!(slope, 0.00001, Inf)
 
-    dt = Float64(config.timestepsecs)
-    flow_parameter_set =
-        ManningFlowParameters(slope, n_river, dl, width, iterate, tstep, dt)
+    flow_parameter_set = ManningFlowParameters(slope, n_river, dl, width)
     parameters =
         RiverFlowParameters(; flow = flow_parameter_set, bankfull_depth = bankfull_depth)
     return parameters
@@ -136,6 +134,7 @@ function RiverFlowBC(n, reservoir, reservoir_index, lake, lake_index)
 end
 
 @with_kw struct SurfaceFlowRiver{T, R, L, A}
+    timestepping::TimeStepping{T}
     boundary_conditions::RiverFlowBC{T, R, L}
     parameters::RiverFlowParameters{T}
     variables::FlowVariables{T}
@@ -152,26 +151,35 @@ function SurfaceFlowRiver(
     reservoir,
     lake_index,
     lake,
-    iterate,
-    tstep,
 )
-    @info "Kinematic wave approach is used for river flow." iterate
-    if tstep > 0
-        @info "Using a fixed sub-timestep (seconds) $tstep for kinematic wave river flow."
-    end
     n = length(inds)
+
+    adaptive = get(config.model, "kin_wave_iteration", false)::Bool
+    @info "Kinematic wave approach is used for river flow, adaptive timestepping = $adaptive."
+
+    if adaptive
+        dt = zeros(n)
+    else
+        dt_fixed = get(config.model, "kw_river_tstep", 900.0)
+        @info "Using a fixed sub-timestep (seconds) $dt_fixed for kinematic wave river flow."
+        dt = dt_fixed
+    end
+
+    timestepping = TimeStepping(; dt, adaptive)
+
     do_water_demand = haskey(config.model, "water_demand")
     allocation = do_water_demand ? AllocationRiver(n) : NoAllocationRiver{Float}()
 
     variables = FlowVariables(n)
-    parameters = RiverFlowParameters(nc, config, inds, dl, width, iterate, tstep)
-    bc = RiverFlowBC(n, reservoir, reservoir_index, lake, lake_index)
+    parameters = RiverFlowParameters(nc, config, inds, dl, width)
+    boundary_conditions = RiverFlowBC(n, reservoir, reservoir_index, lake, lake_index)
 
     sf_river = SurfaceFlowRiver(;
-        boundary_conditions = bc,
-        parameters = parameters,
-        variables = variables,
-        allocation = allocation,
+        timestepping,
+        boundary_conditions,
+        parameters,
+        variables,
+        allocation,
     )
 
     return sf_river
@@ -197,17 +205,13 @@ end
 end
 
 @with_kw struct SurfaceFlowLand{T}
+    timestepping::TimeStepping{T}
     boundary_conditions::LandFlowBC{T}
     parameters::ManningFlowParameters{T}
     variables::LandFlowVariables{T}
 end
 
-function SurfaceFlowLand(nc, config, inds; sl, dl, width, iterate, tstep)
-    @info "Kinematic wave approach is used for overland flow." iterate
-    if tstep > 0
-        @info "Using a fixed sub-timestep (seconds) $tstep for kinematic wave overland flow."
-    end
-
+function SurfaceFlowLand(nc, config, inds; sl, dl, width)
     n_land = ncread(
         nc,
         config,
@@ -216,29 +220,78 @@ function SurfaceFlowLand(nc, config, inds; sl, dl, width, iterate, tstep)
         defaults = 0.072,
         type = Float,
     )
+
     n = length(inds)
+    adaptive = get(config.model, "kin_wave_iteration", false)::Bool
+    @info "Kinematic wave approach is used for river flow, adaptive timestepping = $adaptive."
 
+    if adaptive
+        dt = zeros(n)
+    else
+        dt_fixed = get(config.model, "kw_river_tstep", 3600.0)
+        @info "Using a fixed sub-timestep (seconds) $dt_fixed for kinematic wave river flow."
+        dt = dt_fixed
+    end
+
+    timestepping = TimeStepping(; dt, adaptive)
     variables = LandFlowVariables(; flow = FlowVariables(n), to_river = zeros(Float, n))
     dt = Float64(config.timestepsecs)
-    parameters = ManningFlowParameters(sl, n_land, dl, width, iterate, tstep, dt)
-    bc = LandFlowBC(; inwater = zeros(Float, n))
-    sf_land = SurfaceFlowLand(;
-        boundary_conditions = bc,
-        variables = variables,
-        parameters = parameters,
-    )
+    parameters = ManningFlowParameters(sl, n_land, dl, width)
+    boundary_conditions = LandFlowBC(; inwater = zeros(Float, n))
+    sf_land = SurfaceFlowLand(; timestepping, boundary_conditions, variables, parameters)
 
     return sf_land
 end
 
-function update!(sf::SurfaceFlowLand, network, frac_toriver, dt)
+function surfaceflow_land_update!(sf::SurfaceFlowLand, network, frac_toriver, dt)
     (; subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes) = network
 
-    (; inwater) = sf.boundary_conditions
-    (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, width, dl) = sf.parameters
-    (; h, h_av, q, q_av, qin, qlat, volume, to_river) = sf.variables
+    (; beta, alpha, width, dl) = sf.parameters
+    (; h, h_av, q, q_av, qin, qlat, to_river) = sf.variables
 
     ns = length(subdomain_order)
+    qin .= 0.0
+    for k in 1:ns
+        threaded_foreach(eachindex(subdomain_order[k]); basesize = 1) do i
+            m = subdomain_order[k][i]
+            for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
+                # for a river cell without a reservoir or lake part of the upstream
+                # surface flow goes to the river (frac_toriver) and part goes to the
+                # surface flow reservoir (1.0 - frac_toriver), upstream nodes with a
+                # reservoir or lake are excluded
+                to_river[v] +=
+                    sum_at(
+                        i -> q[i] * frac_toriver[i],
+                        upstream_nodes[n],
+                        eltype(to_river),
+                    ) * dt
+                if width[v] > 0.0
+                    qin[v] = sum_at(
+                        i -> q[i] * (1.0 - frac_toriver[i]),
+                        upstream_nodes[n],
+                        eltype(q),
+                    )
+                end
+
+                q[v] = kinematic_wave(qin[v], q[v], qlat[v], alpha[v], beta, dt, dl[v])
+
+                # update h, only if surface width > 0.0
+                if width[v] > 0.0
+                    crossarea = alpha[v] * pow(q[v], beta)
+                    h[v] = crossarea / width[v]
+                end
+                q_av[v] += q[v] * dt
+                h_av[v] += h[v] * dt
+            end
+        end
+    end
+end
+
+function update!(sf::SurfaceFlowLand, network, frac_toriver, dt)
+    (; inwater) = sf.boundary_conditions
+    (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, width, dl) = sf.parameters
+    (; h, h_av, q_av, qlat, volume, to_river) = sf.variables
+    (; adaptive) = sf.timestepping
 
     @. alpha_term = pow(mannings_n / sqrt(slope), beta)
     # use fixed alpha value based flow width
@@ -249,52 +302,23 @@ function update!(sf::SurfaceFlowLand, network, frac_toriver, dt)
     h_av .= 0.0
     to_river .= 0.0
 
-    dt_s, its = stable_timestep(sf, dt)
-    for _ in 1:its
-        qin .= 0.0
-        for k in 1:ns
-            threaded_foreach(eachindex(subdomain_order[k]); basesize = 1) do i
-                m = subdomain_order[k][i]
-                for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
-                    # for a river cell without a reservoir or lake part of the upstream
-                    # surface flow goes to the river (frac_toriver) and part goes to the
-                    # surface flow reservoir (1.0 - frac_toriver), upstream nodes with a
-                    # reservoir or lake are excluded
-                    to_river[v] += sum_at(
-                        i -> q[i] * frac_toriver[i],
-                        upstream_nodes[n],
-                        eltype(to_river),
-                    )
-                    if width[v] > 0.0
-                        qin[v] = sum_at(
-                            i -> q[i] * (1.0 - frac_toriver[i]),
-                            upstream_nodes[n],
-                            eltype(q),
-                        )
-                    end
-
-                    q[v] =
-                        kinematic_wave(qin[v], q[v], qlat[v], alpha[v], beta, dt_s, dl[v])
-
-                    # update h, only if surface width > 0.0
-                    if width[v] > 0.0
-                        crossarea = alpha[v] * pow(q[v], beta)
-                        h[v] = crossarea / width[v]
-                    end
-                    q_av[v] += q[v]
-                    h_av[v] += h[v]
-                end
-            end
+    t = 0.0
+    while t < dt
+        dt_s = adaptive ? stable_timestep(sf, 0.02) : sf.timestepping.dt
+        if t + dt_s > dt
+            dt_s = dt - t
         end
+        surfaceflow_land_update!(sf, network, frac_toriver, dt_s)
+        t = t + dt_s
     end
-    q_av ./= its
-    h_av ./= its
-    to_river ./= its
+    q_av ./= dt
+    h_av ./= dt
+    to_river ./= dt
     volume .= dl .* width .* h
     return nothing
 end
 
-function update!(sf::SurfaceFlowRiver, network, doy, dt)
+function surfaceflow_river_update!(sf::SurfaceFlowRiver, network, doy, dt)
     (; graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes) = network
 
     (;
@@ -308,17 +332,102 @@ function update!(sf::SurfaceFlowRiver, network, doy, dt)
         inflow_wb,
     ) = sf.boundary_conditions
 
-    (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, width, dl, bankfull_depth) =
-        sf.parameters
+    (; beta, alpha, width, dl) = sf.parameters
     (; h, h_av, q, q_av, qin, qlat, volume) = sf.variables
 
     ns = length(subdomain_order)
+    qin .= 0.0
+    for k in 1:ns
+        threaded_foreach(eachindex(subdomain_order[k]); basesize = 1) do i
+            m = subdomain_order[k][i]
+            for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
+                # sf.qin by outflow from upstream reservoir or lake location is added
+                qin[v] += sum_at(q, upstream_nodes[n])
+                # Inflow supply/abstraction is added to qlat (divide by flow length)
+                # If inflow < 0, abstraction is limited
+                if inflow[v] < 0.0
+                    max_abstract =
+                        min((inwater[v] + qin[v] + volume[v] / dt) * 0.80, -inflow[v])
+                    _inflow = -max_abstract / sf.dl[v]
+                else
+                    _inflow = inflow[v] / dl[v]
+                end
+                _inflow -= abstraction[v] / dl[v]
+
+                q[v] = kinematic_wave(
+                    qin[v],
+                    q[v],
+                    qlat[v] + _inflow,
+                    alpha[v],
+                    beta,
+                    dt,
+                    dl[v],
+                )
+
+                if !isnothing(reservoir) && reservoir_index[v] != 0
+                    # run reservoir model and copy reservoir outflow to inflow (qin) of
+                    # downstream river cell
+                    i = reservoir_index[v]
+                    update!(reservoir, i, q[v] + inflow_wb[v], dt)
+
+                    downstream_nodes = outneighbors(graph, v)
+                    n_downstream = length(downstream_nodes)
+                    if n_downstream == 1
+                        j = only(downstream_nodes)
+                        qin[j] = reservoir.variables.outflow[i]
+                    elseif n_downstream == 0
+                        error(
+                            """A reservoir without a downstream river node is not supported.
+                            Add a downstream river node or move the reservoir to an upstream node (model schematization).
+                            """,
+                        )
+                    else
+                        error("bifurcations not supported")
+                    end
+
+                elseif !isnothing(lake) && lake_index[v] != 0
+                    # run lake model and copy lake outflow to inflow (qin) of downstream river
+                    # cell
+                    i = lake_index[v]
+                    update!(lake, i, q[v] + inflow_wb[v], doy, dt)
+
+                    downstream_nodes = outneighbors(graph, v)
+                    n_downstream = length(downstream_nodes)
+                    if n_downstream == 1
+                        j = only(downstream_nodes)
+                        qin[j] = lake.variables.outflow[i]
+                    elseif n_downstream == 0
+                        error(
+                            """A lake without a downstream river node is not supported.
+                            Add a downstream river node or move the lake to an upstream node (model schematization).
+                            """,
+                        )
+                    else
+                        error("bifurcations not supported")
+                    end
+                end
+                # update h
+                crossarea = alpha[v] * pow(q[v], beta)
+                h[v] = crossarea / width[v]
+                volume[v] = dl[v] * width[v] * h[v]
+                q_av[v] += q[v] * dt
+                h_av[v] += h[v] * dt
+            end
+        end
+    end
+end
+
+function update!(sf::SurfaceFlowRiver, network, doy, dt)
+    (; reservoir, lake, inwater) = sf.boundary_conditions
+
+    (; alpha_term, mannings_n, slope, beta, alpha_pow, alpha, width, dl, bankfull_depth) =
+        sf.parameters
+    (; h, h_av, q_av, qlat, volume) = sf.variables
+    (; adaptive) = sf.timestepping
 
     @. alpha_term = pow(mannings_n / sqrt(slope), beta)
     # use fixed alpha value based on 0.5 * bankfull_depth
     @. alpha = alpha_term * pow(width + bankfull_depth, alpha_pow)
-
-    # move qlat to boundary conditions function?
     @. qlat = inwater / dl
 
     q_av .= 0.0
@@ -336,124 +445,43 @@ function update!(sf::SurfaceFlowRiver, network, doy, dt)
         lake.variables.actevap .= 0.0
     end
 
-    dt_s, its = stable_timestep(sf, dt)
-    for _ in 1:its
-        qin .= 0.0
-        for k in 1:ns
-            threaded_foreach(eachindex(subdomain_order[k]); basesize = 1) do i
-                m = subdomain_order[k][i]
-                for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
-                    # sf.qin by outflow from upstream reservoir or lake location is added
-                    qin[v] += sum_at(q, upstream_nodes[n])
-                    # Inflow supply/abstraction is added to qlat (divide by flow length)
-                    # If inflow < 0, abstraction is limited
-                    if inflow[v] < 0.0
-                        max_abstract =
-                            min((inwater[v] + qin[v] + volume[v] / dt_s) * 0.80, -inflow[v])
-                        _inflow = -max_abstract / sf.dl[v]
-                    else
-                        _inflow = inflow[v] / dl[v]
-                    end
-                    _inflow -= abstraction[v] / dl[v]
-
-                    q[v] = kinematic_wave(
-                        qin[v],
-                        q[v],
-                        qlat[v] + _inflow,
-                        alpha[v],
-                        beta,
-                        dt_s,
-                        dl[v],
-                    )
-
-                    if !isnothing(reservoir) && reservoir_index[v] != 0
-                        # run reservoir model and copy reservoir outflow to inflow (qin) of
-                        # downstream river cell
-                        i = reservoir_index[v]
-                        update!(reservoir, i, q[v] + inflow_wb[v], dt_s)
-
-                        downstream_nodes = outneighbors(graph, v)
-                        n_downstream = length(downstream_nodes)
-                        if n_downstream == 1
-                            j = only(downstream_nodes)
-                            qin[j] = reservoir.variables.outflow[i]
-                        elseif n_downstream == 0
-                            error(
-                                """A reservoir without a downstream river node is not supported.
-                                Add a downstream river node or move the reservoir to an upstream node (model schematization).
-                                """,
-                            )
-                        else
-                            error("bifurcations not supported")
-                        end
-
-                    elseif !isnothing(lake) && lake_index[v] != 0
-                        # run lake model and copy lake outflow to inflow (qin) of downstream river
-                        # cell
-                        i = lake_index[v]
-                        update!(lake, i, q[v] + inflow_wb[v], doy, dt_s)
-
-                        downstream_nodes = outneighbors(graph, v)
-                        n_downstream = length(downstream_nodes)
-                        if n_downstream == 1
-                            j = only(downstream_nodes)
-                            qin[j] = lake.variables.outflow[i]
-                        elseif n_downstream == 0
-                            error(
-                                """A lake without a downstream river node is not supported.
-                                Add a downstream river node or move the lake to an upstream node (model schematization).
-                                """,
-                            )
-                        else
-                            error("bifurcations not supported")
-                        end
-                    end
-
-                    # update h
-                    crossarea = alpha[v] * pow(q[v], beta)
-                    h[v] = crossarea / width[v]
-                    volume[v] = dl[v] * width[v] * h[v]
-                    q_av[v] += q[v]
-                    h_av[v] += h[v]
-                end
-            end
+    t = 0.0
+    while t < dt
+        dt_s = adaptive ? stable_timestep(sf, 0.05) : sf.timestepping.dt
+        if t + dt_s > dt
+            dt_s = dt - t
         end
+        surfaceflow_river_update!(sf, network, doy, dt_s)
+        t = t + dt_s
     end
-    q_av ./= its
-    h_av ./= its
+    q_av ./= dt
+    h_av ./= dt
     volume .= dl .* width .* h
     return nothing
 end
 
-function stable_timestep(sf::S, dt) where {S <: Union{SurfaceFlowLand, SurfaceFlowRiver}}
-    (; q, celerity) = sf.variables
-    (; alpha, beta, dl, its, kinwave_it) = sf.parameters
+function stable_timestep(sf::S, p) where {S <: Union{SurfaceFlowLand, SurfaceFlowRiver}}
+    (; q) = sf.variables
+    (; alpha, beta, dl) = sf.parameters
+    (; dt) = sf.timestepping
 
     n = length(q)
-    # two options for iteration, fixed or based on courant number.
-    if kinwave_it
-        if its > 0
-            n_iterations = its
-        else
-            # calculate celerity
-            courant = zeros(n)
-            for v in 1:n
-                if q[v] > 0.0
-                    celerity[v] = 1.0 / (alpha[v] * beta * pow(q[v], (beta - 1.0)))
-                    courant[v] = (dt / dl[v]) * celerity[v]
-                end
-            end
-            filter!(x -> x ≠ 0.0, courant)
-            n_iterations =
-                isempty(courant) ? 1 : ceil(Int, (1.25 * quantile!(courant, 0.95)))
+    dt .= Inf
+    for i in 1:n
+        if q[i] > 0.0
+            c = 1.0 / (alpha[i] * beta * pow(q[i], (beta - 1.0)))
+            dt[i] = (dl[i] / c)
         end
+    end
+    dt_filtered = filter(x -> !isinf(x), dt)
+
+    if !isempty(dt_filtered)
+        dt_s = quantile!(dt_filtered, p)
     else
-        n_iterations = 1
+        dt_s = 600.0
     end
 
-    # sub time step
-    dt_sub = dt / n_iterations
-    return dt_sub, n_iterations
+    return dt_s
 end
 
 @get_units @grid_loc struct KhExponential{T}
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 2d7c23906..008061e4f 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -138,8 +138,6 @@ function initialize_sbm_gwf_model(config::Config)
             sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
-            iterate = kinwave_it,
-            tstep = kw_land_tstep,
         )
     elseif land_routing == "local-inertial"
         index_river_nf = rev_inds_riv[inds] # not filtered (with zeros)
@@ -177,8 +175,6 @@ function initialize_sbm_gwf_model(config::Config)
             reservoir = reservoirs,
             lake_index = lakeindex,
             lake = lakes,
-            iterate = kinwave_it,
-            tstep = kw_river_tstep,
         )
     elseif river_routing == "local-inertial"
         rf, nodes_at_link = ShallowWaterRiver(
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 604bb1052..154adca9e 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -159,8 +159,6 @@ function initialize_sbm_model(config::Config)
             sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
-            iterate = kinwave_it,
-            tstep = kw_land_tstep,
         )
     elseif land_routing == "local-inertial"
         index_river_nf = rev_inds_riv[inds] # not filtered (with zeros)
@@ -196,8 +194,6 @@ function initialize_sbm_model(config::Config)
             reservoir = reservoirs,
             lake_index = lakeindex,
             lake = lakes,
-            iterate = kinwave_it,
-            tstep = kw_river_tstep,
         )
     elseif river_routing == "local-inertial"
         rf, nodes_at_link = ShallowWaterRiver(
diff --git a/test/bmi.jl b/test/bmi.jl
index aa7719c77..d11ad473f 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -17,8 +17,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) == 206
-            @test BMI.get_output_item_count(model) == 206
+            @test BMI.get_input_item_count(model) == 204
+            @test BMI.get_output_item_count(model) == 204
             to_check = [
                 "vertical.soil.parameters.nlayers",
                 "vertical.soil.parameters.theta_r",
@@ -82,7 +82,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
                 "lateral.river.variables.q",
                 zeros(Float, 3),
                 [1, 100, 5617],
-            ) ≈ [0.623325399343309, 5.227139951657074, 0.027942874327781947]
+            ) ≈ [0.6525631197206111, 7.493760826794606, 0.02319714614721354]
             BMI.set_value(
                 model,
                 "vertical.soil.variables.zi",
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 8a1df28b2..2d138ed06 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -14,26 +14,26 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-02T00:00:00")
-    @test row.Q ≈ 8.15299947254324f0
+    @test row.Q ≈ 14.283603959245331f0
     @test row.volume ≈ 2.7535003939625636f7
     @test row.temp_bycoord ≈ 2.390000104904175f0
     @test row.vwc_layer2_bycoord ≈ 0.25938809638672006f0
     @test row.temp_byindex ≈ 2.390000104904175f0
-    @test row.Q_6336050 ≈ 0.006583064321841488f0
-    @test row.Q_6336510 ≈ 0.029864230092642642f0
-    @test row.Q_6836100 ≈ 0.19995488963854305f0
-    @test row.Q_6336500 ≈ 0.006277726622788425f0
-    @test row.Q_6836190 ≈ 0.0031262850749354237f0
-    @test row.Q_6336800 ≈ 0.008278375560053742f0
-    @test row.Q_6336900 ≈ 0.0066141980189014385f0
-    @test row.Q_6336930 ≈ 0.09141703511009937f0
-    @test row.Q_6336910 ≈ 0.007475453481320056f0
-    @test row.Q_6136500 ≈ 0.001834989281902289f0
-    @test row.Q_6136520 ≈ 0.0022266031120691397f0
-    @test row.Q_6136150 ≈ 0.006310361139139334f0
-    @test row.Q_6136151 ≈ 0.007946301730645885f0
-    @test row.Q_6136160 ≈ 3.927719795530719f0
-    @test row.Q_6136202 ≈ 1.4162246003743886f0
+    @test row.Q_6336050 ≈ 0.0077488610359143706f0
+    @test row.Q_6336510 ≈ 0.029475452452069627f0
+    @test row.Q_6836100 ≈ 0.018908801563313697f0
+    @test row.Q_6336500 ≈ 0.006568230884299706f0
+    @test row.Q_6836190 ≈ 0.004495680971248625f0
+    @test row.Q_6336800 ≈ 0.00916756751611816f0
+    @test row.Q_6336900 ≈ 0.008200812406926841f0
+    @test row.Q_6336930 ≈ 0.026477954532813916f0
+    @test row.Q_6336910 ≈ 0.009127079213306247f0
+    @test row.Q_6136500 ≈ 0.001154334630885042f0
+    @test row.Q_6136520 ≈ 0.0014432161225527867f0
+    @test row.Q_6136150 ≈ 0.0073146444422822225f0
+    @test row.Q_6136151 ≈ 0.006614628762328246f0
+    @test row.Q_6136160 ≈ 6.186486795085275f0
+    @test row.Q_6136202 ≈ 1.8472398588556238f0
     @test row.recharge_1 ≈ -0.0020800523945940217f0
 end
 
@@ -41,26 +41,26 @@ end
     ds = model.writer.dataset_scalar
     @test ds["time"][1] == DateTime("2000-01-02T00:00:00")
     @test ds["Q"][:][1:20] ≈ [
-        0.7425387f0,
-        1.4162246f0,
-        1.4425076f0,
-        1.4044669f0,
-        5.738109f0,
-        2.7616737f0,
-        2.1128905f0,
-        4.105428f0,
-        0.008651769f0,
-        3.9277198f0,
-        4.069447f0,
-        0.006356805f0,
-        0.007946302f0,
-        0.008135906f0,
-        0.0037393502f0,
-        0.70888275f0,
-        0.0024000728f0,
-        1.3347782f0,
-        3.8374817f0,
-        1.676597f0,
+        0.79544294f0,
+        1.8472399f0,
+        1.8459954f0,
+        1.5499605f0,
+        13.334657f0,
+        4.865496f0,
+        0.09143141f0,
+        4.437221f0,
+        0.010515818f0,
+        6.1864867f0,
+        6.1764274f0,
+        0.0060072034f0,
+        0.0066146287f0,
+        0.006407934f0,
+        0.00416659f0,
+        1.2095966f0,
+        0.002213421f0,
+        1.832764f0,
+        8.601765f0,
+        2.7730224f0,
     ]
     @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
     @test ds["temp_index"][:] ≈ [2.39f0]
@@ -84,7 +84,7 @@ end
     @test snow.variables.snow_storage[5] ≈ 3.768513390588815f0
     @test mean(snow.variables.snow_storage) ≈ 0.038019723676094325f0
     @test sbm.variables.total_storage[50063] ≈ 559.9035608052374f0
-    @test sbm.variables.total_storage[429] ≈ 597.4578475404879f0 # river cell
+    @test sbm.variables.total_storage[429] ≈ 597.1603591432968f0 # river cell
 end
 
 # run the second timestep
@@ -103,7 +103,7 @@ Wflow.run_timestep!(model)
     @test snow.variables.snow_storage[5] ≈ 3.843412524052313f0
     @test mean(snow.variables.snow_storage) ≈ 0.03461317061870949f0
     @test sbm.variables.total_storage[50063] ≈ 560.0152135062889f0
-    @test sbm.variables.total_storage[429] ≈ 617.2238533241972f0 # river cell
+    @test sbm.variables.total_storage[429] ≈ 616.8798940139915f0 # river cell
 end
 
 @testset "subsurface flow" begin
@@ -116,7 +116,7 @@ end
 
 @testset "overland flow" begin
     q = model.lateral.land.variables.q_av
-    @test sum(q) ≈ 291.4923871784623f0
+    @test sum(q) ≈ 285.59889671823953f0
     @test q[26625] ≈ 0.0
     @test q[39308] ≈ 0.0
     @test q[network.land.order[end]] ≈ 1.0f-30
@@ -124,16 +124,16 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.variables.q_av
-    @test sum(q) ≈ 3625.0013368279815f0
-    @test q[1622] ≈ 0.0006503254947860838f0
-    @test q[43] ≈ 12.06416878694095f0
-    @test q[network.river.order[end]] ≈ 0.039200124520463835f0
+    @test sum(q) ≈ 3848.832553197247f0
+    @test q[1622] ≈ 0.0007520477205381629f0
+    @test q[43] ≈ 11.936893577401598f0
+    @test q[network.river.order[end]] ≈ 0.04418878169426842f0
 end
 
 @testset "reservoir simple" begin
     res = model.lateral.river.boundary_conditions.reservoir
     @test res.variables.outflow[1] ≈ 0.21750000119148086f0
-    @test res.boundary_conditions.inflow[1] ≈ 43.18479982574888f0
+    @test res.boundary_conditions.inflow[1] ≈ 44.312772458823474f0
     @test res.variables.volume[1] ≈ 2.751299001489657f7
     @test res.boundary_conditions.precipitation[1] ≈ 0.17999997735023499f0
     @test res.boundary_conditions.evaporation[1] ≈ 0.5400000810623169f0
@@ -168,10 +168,10 @@ end
 
 @testset "river flow at basin outlets and downstream of one pit" begin
     q = model.lateral.river.variables.q_av
-    @test q[4009] ≈ 8.543145028037452f0 # pit/ outlet, CartesianIndex(141, 228)
+    @test q[4009] ≈ 8.556311783589425f0 # pit/ outlet, CartesianIndex(141, 228)
     @test q[4020] ≈ 0.006779014715290862f0 # downstream of pit 4009, CartesianIndex(141, 229)
-    @test q[2508] ≈ 150.5640617045796f0 # pit/ outlet
-    @test q[5808] ≈ 0.12438899196040518f0 # pit/ outlet
+    @test q[2508] ≈ 150.5479054707659f0 # pit/ outlet
+    @test q[5808] ≈ 0.12337347858092176f0 # pit/ outlet
 end
 
 # test changing forcing and cyclic LAI parameter
@@ -219,7 +219,7 @@ Wflow.run_timestep!(model)
 
 @testset "river inflow (cyclic)" begin
     @test model.lateral.river.boundary_conditions.inflow[44] ≈ 0.75
-    @test model.lateral.river.variables.q_av[44] ≈ 10.723729440690567f0
+    @test model.lateral.river.variables.q_av[44] ≈ 10.554267976107754f0
 end
 
 # test fixed forcing (precipitation = 2.5)
@@ -265,14 +265,14 @@ Wflow.run_timestep!(model)
 
 @testset "river flow and depth (local inertial)" begin
     q = model.lateral.river.variables.q_av
-    @test sum(q) ≈ 3922.0644366362544f0
-    @test q[1622] ≈ 7.315676375562105f-5
-    @test q[43] ≈ 11.92787156357907f0
-    @test q[501] ≈ 3.57855182713785f0
+    @test sum(q) ≈ 3854.717278465037f0
+    @test q[1622] ≈ 7.311274274855274f-5
+    @test q[43] ≈ 11.716766734364437f0
+    @test q[501] ≈ 3.4819773071884716f0
     h = model.lateral.river.variables.h_av
     @test h[1622] ≈ 0.001987887644883981f0
-    @test h[43] ≈ 0.4366415244811759f0
-    @test h[501] ≈ 0.057317706869865745f0
+    @test h[43] ≈ 0.43311924038778815f0
+    @test h[501] ≈ 0.05635210581824346f0
     q_channel = model.lateral.river.variables.q_channel_av
     @test q ≈ q_channel
 end
@@ -419,16 +419,16 @@ Wflow.run_timestep!(model)
 
 @testset "river flow (local inertial) with floodplain schematization simulation" begin
     q = model.lateral.river.variables.q_av
-    @test sum(q) ≈ 3910.4728717811836f0
-    @test q[1622] ≈ 7.315676384849305f-5
-    @test q[43] ≈ 11.92787156357908f0
-    @test q[501] ≈ 3.510668846752431f0
-    @test q[5808] ≈ 0.002223993845806248f0
+    @test sum(q) ≈ 3843.944494991296f0
+    @test q[1622] ≈ 7.311274278896417f-5
+    @test q[43] ≈ 11.716766734364418f0
+    @test q[501] ≈ 3.424364314225289f0
+    @test q[5808] ≈ 0.0022274679501657693f0
     h = model.lateral.river.variables.h_av
     @test h[1622] ≈ 0.001987887580593841f0
-    @test h[43] ≈ 0.436641524481545f0
-    @test h[501] ≈ 0.05670770509802258f0
-    @test h[5808] ≈ 0.005929945681367346f0
+    @test h[43] ≈ 0.433119240388070f0
+    @test h[501] ≈ 0.055832770820860404f0
+    @test h[5808] ≈ 0.005935591961908253f0
 end
 
 # set boundary condition local inertial routing from netCDF file
@@ -440,15 +440,15 @@ Wflow.run_timestep!(model)
 
 @testset "change boundary condition for local inertial routing (including floodplain)" begin
     q = model.lateral.river.variables.q_av
-    @test sum(q) ≈ 3910.683449719468f0
-    @test q[1622] ≈ 7.315757521099307f-5
-    @test q[43] ≈ 11.927871563591228f0
-    @test q[501] ≈ 3.5106678593721496f0
+    @test sum(q) ≈ 3844.1544889903134f0
+    @test q[1622] ≈ 7.311151138631112f-5
+    @test q[43] ≈ 11.716766734416717f0
+    @test q[501] ≈ 3.424329413571391f0
     @test q[5808] ≈ 0.060518234525259465f0
     h = model.lateral.river.variables.h_av
     @test h[1622] ≈ 0.0019878952928530183f0
-    @test h[43] ≈ 0.4366415249636809f0
-    @test h[501] ≈ 0.056707564314724804f0
+    @test h[43] ≈ 0.4331192403230577f0
+    @test h[501] ≈ 0.0558281185092927f0
     @test h[5808] ≈ 2.0000006940603936f0
 end
 Wflow.close_files(model; delete_output = false)
@@ -535,9 +535,9 @@ Wflow.close_files(model; delete_output = false)
         Wflow.run_timestep!(model)
         Wflow.run_timestep!(model)
         q = model.lateral.river.variables.q_av
-        @test sum(q) ≈ 3159.38300016008f0
-        @test q[1622] ≈ 0.0005972577112819149f0
-        @test q[43] ≈ 10.017642376280731f0
+        @test sum(q) ≈ 3302.1390089922525f0
+        @test q[1622] ≈ 0.0006990391393246408f0
+        @test q[43] ≈ 9.673592182882691f0
     end
 
     Wflow.close_files(model; delete_output = false)
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 9dd58e73d..7fd9eb4d3 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -13,7 +13,7 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-06-01T00:00:00")
-    @test row.Q_av ≈ 0.01620324716944374f0
+    @test row.Q_av ≈ 0.01619703129434486f0
     @test row.head ≈ 1.6471323360175287f0
 end
 
@@ -40,18 +40,18 @@ end
 
 @testset "overland flow (kinematic wave)" begin
     q = model.lateral.land.variables.q_av
-    @test sum(q) ≈ 2.229860508650628f-7
+    @test sum(q) ≈ 2.2321111203610908f-7
 end
 
 @testset "river domain (kinematic wave)" begin
     q = model.lateral.river.variables.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.035443370536496675f0
-    @test q[6] ≈ 0.008031554512314907f0
+    @test sum(q) ≈ 0.035468154534622556f0
+    @test q[6] ≈ 0.00803825101724232f0
     @test river.variables.volume[6] ≈ 4.532124903256408f0
-    @test river.boundary_conditions.inwater[6] ≈ 0.0004073892212290558f0
+    @test river.boundary_conditions.inwater[6] ≈ 0.00040826084140456616f0
     @test q[13] ≈ 0.0006017024138583771f0
-    @test q[network.river.order[end]] ≈ 0.008559590281509943f0
+    @test q[network.river.order[end]] ≈ 0.00856866488665273f0
 end
 
 @testset "groundwater" begin
@@ -158,7 +158,7 @@ end
 
 @testset "overland flow warm start (kinematic wave)" begin
     q = model.lateral.land.variables.q_av
-    @test sum(q) ≈ 1.4589771292158736f-5
+    @test sum(q) ≈ 1.4224503548471601f-5
 end
 
 @testset "river domain warm start (kinematic wave)" begin
@@ -167,7 +167,7 @@ end
     @test sum(q) ≈ 0.01191742350356312f0
     @test q[6] ≈ 0.0024353072305122064f0
     @test river.variables.volume[6] ≈ 2.2277585577366357f0
-    @test river.boundary_conditions.inwater[6] ≈ -1.3019072795599315f-5
+    @test river.boundary_conditions.inwater[6] ≈ -1.3042629584651168f-5
     @test q[13] ≈ 7.332742814063803f-5
     @test q[network.river.order[end]] ≈ 0.002472526149620472f0
 end
diff --git a/test/run_sbm_piave.jl b/test/run_sbm_piave.jl
index 450f31767..ac666bd0a 100644
--- a/test/run_sbm_piave.jl
+++ b/test/run_sbm_piave.jl
@@ -27,52 +27,52 @@ Wflow.close_files(model; delete_output = false)
 @testset "piave with and without water demand" begin
     idx = 1:3:28
     @test q_demand[idx] ≈ [
-        218.52770747627918f0,
-        193.02890386240665f0,
-        271.68778616960685f0,
-        161.80734173386602f0,
-        164.81279958487485f0,
-        150.4149716788231f0,
-        121.02659706677429f0,
-        108.15426854851842f0,
-        174.63022487569546f0,
-        108.20883498755789f0,
+        218.52013823809472f0,
+        193.0134951603773f0,
+        272.4111837647947f0,
+        161.88264628787172f0,
+        164.8199089671644f0,
+        150.2681168314876f0,
+        121.20070337007452f0,
+        108.10106381132412f0,
+        175.13799714754256f0,
+        108.26190463186364f0,
     ]
     @test q_[idx] ≈ [
-        219.9042687883228f0,
-        197.73238933696658f0,
-        276.99163278840734f0,
-        165.74244080863346f0,
-        172.99856395134805f0,
-        153.91963517651158f0,
-        128.52653903788593f0,
-        112.02923578000491f0,
-        178.19599851038708f0,
-        109.99414262238396f0,
+        219.87655632704903f0,
+        197.7038754807009f0,
+        277.7110869134211f0,
+        165.79913971520423f0,
+        173.04466296857905f0,
+        153.84187794694486f0,
+        128.71609293239374f0,
+        112.02394903669563f0,
+        178.8207608992179f0,
+        110.0540286256144f0,
     ]
     @test riv_vol_demand[idx] ≈ [
-        60125.790043761706f0,
-        55476.580177102805f0,
-        62195.90597660078f0,
-        48963.00368970478f0,
-        51871.46897847274f0,
-        43794.737044826295f0,
-        41710.964839545784f0,
-        39623.67826731185f0,
-        44719.74807290461f0,
-        39494.768128540185f0,
+        60108.85346812383f0,
+        55483.137044195726f0,
+        62098.68355844664f0,
+        48989.58100040463f0,
+        51879.53087453071f0,
+        43784.59321641879f0,
+        41700.35166194379f0,
+        39630.55231098758f0,
+        44664.51932077705f0,
+        39490.08957716613f0,
     ]
     @test riv_vol[idx] ≈ [
-        60620.42943610538f0,
-        55927.61042067993f0,
-        62448.21099253601f0,
-        49283.04780863544f0,
-        52293.53301649927f0,
-        44029.56815119591f0,
-        41996.885591977836f0,
-        39751.39389784692f0,
-        44619.74874948371f0,
-        39285.722489743566f0,
+        60612.8152223446f0,
+        55934.74021061718f0,
+        62385.538865370385f0,
+        49305.35863834837f0,
+        52307.51338053202f0,
+        44027.55259408914f0,
+        41998.69770760942f0,
+        39768.94289967264f0,
+        44586.63915460806f0,
+        39296.07358095679f0,
     ]
     @test ssf_vol_demand[idx] ≈ [
         244149.81771426898f0,