diff --git a/Wflow/src/Wflow.jl b/Wflow/src/Wflow.jl
index d8a95d442..1c8989014 100644
--- a/Wflow/src/Wflow.jl
+++ b/Wflow/src/Wflow.jl
@@ -75,6 +75,8 @@ const VERSION =
     VersionNumber(TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"])
 
 const GRAVITATIONAL_ACCELERATION = 9.80665 # m s⁻²
+# local drain direction pit [-]
+const LDD_PIT = 5
 
 mutable struct Clock{T}
     time::T
@@ -139,7 +141,7 @@ include("routing/routing.jl")
 include("domain.jl")
 
 """
-    Model{R <: Routing, L <: AbstractLandModel, M <: AbstractMassBalance, T <: AbstractModelType} <: AbstractModel{T}
+    Model{R <: Routing, L <: AbstractLandModel, M <: AbstractMassBalance, W <: Writer, T <: AbstractModelType} <: AbstractModel{T}
 
 Composite type that represents all different aspects of a Wflow Model, such as the network,
 parameters, clock, configuration and input and output.
@@ -209,8 +211,10 @@ include("routing/subsurface/lateral_subsurface_flow.jl")
 include("routing/subsurface/subsurface_process.jl")
 include("routing/subsurface/boundary_conditions.jl")
 include("routing/surface/reservoir.jl")
+include("routing/surface/floodplain.jl")
+include("routing/surface/surface_flow.jl")
 include("routing/surface/surface_kinwave.jl")
-include("routing/surface/surface_local_inertial.jl")
+include("routing/surface/surface_staggered_scheme.jl")
 include("routing/surface/surface_routing.jl")
 include("routing/surface/surface_process.jl")
 include("demand/water_demand.jl")
diff --git a/Wflow/src/bmi.jl b/Wflow/src/bmi.jl
index 6315678c1..b8f723730 100644
--- a/Wflow/src/bmi.jl
+++ b/Wflow/src/bmi.jl
@@ -312,9 +312,9 @@ function BMI.get_grid_edge_count(model::Model, grid::Int)
     if grid == 3
         return ne(domain.river.network.graph)
     elseif grid == 4
-        return length(domain.land.network.edge_indices.idx_right)
+        return length(domain.land.network.edge_indices.ind_x_up)
     elseif grid == 5
-        return length(domain.land.network.edge_indices.idx_up)
+        return length(domain.land.network.edge_indices.ind_y_up)
     elseif grid in 0:2 || grid == 6
         @warn("edges are not provided for grid type $grid (variables are located at nodes)")
     else
@@ -334,16 +334,16 @@ function BMI.get_grid_edge_nodes(model::Model, grid::Int, edge_nodes::Vector{Int
         edge_nodes[range(2, n; step = 2)] = nodes_at_edge.dst
         return edge_nodes
     elseif grid == 4
-        idx_right = domain.land.network.edge_indices.idx_right
+        ind_x_up = domain.land.network.edge_indices.ind_x_up
         edge_nodes[range(1, n; step = 2)] = 1:m
-        idx_right[idx_right .== m + 1] .= -999
-        edge_nodes[range(2, n; step = 2)] = idx_right
+        ind_x_up[ind_x_up .== m + 1] .= -999
+        edge_nodes[range(2, n; step = 2)] = ind_x_up
         return edge_nodes
     elseif grid == 5
-        idx_up = domain.land.network.edge_indices.idx_up
+        ind_y_up = domain.land.network.edge_indices.ind_y_up
         edge_nodes[range(1, n; step = 2)] = 1:m
-        idx_up[idx_up .== m + 1] .= -999
-        edge_nodes[range(2, n; step = 2)] = idx_up
+        ind_y_up[ind_y_up .== m + 1] .= -999
+        edge_nodes[range(2, n; step = 2)] = ind_y_up
         return edge_nodes
     elseif grid in 0:2 || grid == 6
         @warn("edges are not provided for grid type $grid (variables are located at nodes)")
@@ -408,7 +408,7 @@ Return the grid element type of a model variable (PropertyLens `var`) based on a
 function grid_element_type(
     ::T,
     var::PropertyLens,
-) where {T <: Union{LocalInertialRiverFlowModel, LocalInertialOverlandFlowModel}}
+) where {T <: Union{RiverFlowModel{<:LocalInertial}, OverlandFlowModel{<:LocalInertial}}}
     vars = (PropertyLens(x) for x in (:q, :q_average, :qx, :qy))
     element_type = if var in vars
         "edge"
diff --git a/Wflow/src/config_structure.jl b/Wflow/src/config_structure.jl
index 873f18fbc..b4081d168 100644
--- a/Wflow/src/config_structure.jl
+++ b/Wflow/src/config_structure.jl
@@ -6,7 +6,7 @@ For configuration files we use TOML.
 =#
 
 # Option enumerators
-@enumx RoutingType kinematic_wave local_inertial
+@enumx RoutingType kinematic_wave manning_staggered local_inertial
 @enumx ModelType sbm sbm_gwf sediment
 @enumx CalendarType standard gregorian proleptic_gregorian julian noleap _365_day all_leap _366_day _360_day
 @enumx GwfConductivityProfileType uniform exponential
@@ -91,11 +91,14 @@ end
     river_local_inertial_flow__alpha_coefficient::Float64 = 0.7
     land_local_inertial_flow__alpha_coefficient::Float64 = 0.7
     land_local_inertial_flow__theta_coefficient::Float64 = 1.0
-    river_water_flow_threshold__depth = 1e-3
     land_surface_water_flow_threshold__depth = 1e-3
     river_water_flow__froude_limit_flag = true
     land_surface_water_flow__froude_limit_flag = true
+    # Local inertial and staggered manning river flow routing
+    river_water_flow_threshold__depth = 1e-3
     floodplain_1d__flag::Bool = false
+    # Staggered manning river flow routing
+    river_staggered_manning_flow__alpha_coefficient::Float64 = 0.7
     # Groundwater flow
     conductivity_profile::GwfConductivityProfileType.T = GwfConductivityProfileType.uniform
     drain__flag::Bool = false
diff --git a/Wflow/src/domain.jl b/Wflow/src/domain.jl
index 0089798f5..419389b0f 100644
--- a/Wflow/src/domain.jl
+++ b/Wflow/src/domain.jl
@@ -107,7 +107,8 @@ function Domain(dataset::NCDataset, config::Config, ::Union{SbmModel, SbmGwfMode
     if river_routing == RoutingType.kinematic_wave
         @reset network_river.upstream_nodes =
             filter_upstream_nodes(network_river.graph, pits[network_river.indices])
-    elseif river_routing == RoutingType.local_inertial
+    elseif river_routing == RoutingType.local_inertial ||
+           river_routing == RoutingType.manning_staggered
         nodes_at_edge, index_pit = NodesAtEdge(network_river)
         @reset network_river.nodes_at_edge = nodes_at_edge
         @reset network_river.pit_indices = network_river.indices[index_pit]
diff --git a/Wflow/src/mass_balance.jl b/Wflow/src/mass_balance.jl
index ade36dd27..0e1ce064f 100644
--- a/Wflow/src/mass_balance.jl
+++ b/Wflow/src/mass_balance.jl
@@ -125,13 +125,12 @@ function compute_total_storage!(
 end
 
 """
-    get_storage(river_flow_model::LocalInertialRiverFlowModel, i)
-    get_storage(river_flow_model::KinWaveRiverFlowModel, i)
+    get_storage(river_flow_model::RiverFlowModel, i)
 
-Return storage of a river flow model at index `i`. For `LocalInertialRiverFlowModel` floodplain
-storage is added to river storage if an optional floodplain is included.
+Return storage of a river flow model at index `i`, floodplain storage is added to river
+storage if an optional floodplain is included.
 """
-function get_storage(river_flow_model::LocalInertialRiverFlowModel, i)
+function get_storage(river_flow_model::AbstractRiverFlowModel, i)
     (; storage) = river_flow_model.variables
     if isnothing(river_flow_model.floodplain)
         return storage[i]
@@ -140,8 +139,6 @@ function get_storage(river_flow_model::LocalInertialRiverFlowModel, i)
         return total_storage
     end
 end
-get_storage(river_flow_model::KinWaveRiverFlowModel, i) =
-    river_flow_model.variables.storage[i]
 
 """
 Save river (+ floodplain) storage at previous time step as `storage_prev` of river
@@ -291,7 +288,7 @@ compute_flow_balance!(reservoir_model::Nothing, water_balance::NoMassBalance, dt
 
 "Compute water mass balance error and relative error for river kinematic wave routing."
 function compute_flow_balance!(
-    river_flow_model::KinWaveRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:KinematicWave},
     water_balance::MassBalance,
     network::NetworkRiver,
     dt::Float64,
@@ -299,12 +296,16 @@ function compute_flow_balance!(
     (; storage_prev, error, relative_error) = water_balance
     (; inwater, external_inflow, actual_external_abstraction_average, abstraction) =
         river_flow_model.boundary_conditions
-    (; qin_average, q_average, storage) = river_flow_model.variables
+    (; qin_average, q_average) = river_flow_model.variables
 
     for i in eachindex(storage_prev)
         total_in = inwater[i] + qin_average[i] + max(0.0, external_inflow[i])
         total_out = q_average[i] + actual_external_abstraction_average[i] + abstraction[i]
-        storage_rate = (storage[i] - storage_prev[i]) / dt
+        storage = river_flow_model.variables.storage[i]
+        if !isnothing(river_flow_model.floodplain)
+            storage += river_flow_model.floodplain.variables.storage[i]
+        end
+        storage_rate = (storage - storage_prev[i]) / dt
         error[i], relative_error[i] =
             compute_mass_balance_error(total_in, total_out, storage_rate)
     end
@@ -316,7 +317,7 @@ Compute water mass balance error and relative error for river (and floodplain) l
 inertial routing.
 """
 function compute_flow_balance!(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     water_balance::MassBalance,
     network::NetworkRiver,
     dt::Float64,
@@ -326,18 +327,17 @@ function compute_flow_balance!(
         river_flow_model.boundary_conditions
     (; edges_at_node) = network
 
-    q_av_average = river_flow_model.variables.q_average
-    actual_external_abstraction_av_average = actual_external_abstraction_average
+    (; q_average) = river_flow_model.variables
 
     for i in river_flow_model.parameters.active_n
         total_in = 0.0
         total_out = 0.0
-        q_src = sum_at(q_av_average, edges_at_node.src[i])
+        q_src = sum_at(q_average, edges_at_node.src[i])
         total_in, total_out = add_inflow(total_in, total_out, [q_src, inwater[i]])
         total_in += max(0.0, external_inflow[i])
-        q_dst = sum_at(q_av_average, edges_at_node.dst[i])
+        q_dst = sum_at(q_average, edges_at_node.dst[i])
         total_in, total_out = add_outflow(total_in, total_out, q_dst)
-        total_out += actual_external_abstraction_av_average[i] + abstraction[i]
+        total_out += actual_external_abstraction_average[i] + abstraction[i]
         storage = river_flow_model.variables.storage[i]
 
         if !isnothing(river_flow_model.floodplain)
@@ -380,7 +380,7 @@ Compute water mass balance error and relative error for overland flow kinematic
 routing.
 """
 function compute_flow_balance!(
-    overland_flow_model::KinWaveOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     water_balance::MassBalance,
     dt::Float64,
 )
@@ -404,8 +404,8 @@ Compute water mass balance error and relative error for 1D river local inertial
 computed for each land cell (total storage) considering both river and overland flow.
 """
 function compute_flow_balance!(
-    river_flow_model::LocalInertialRiverFlowModel,
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    river_flow_model::RiverFlowModel{<:LocalInertial},
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     water_balance::MassBalance,
     domain::Domain,
     dt::Float64,
@@ -426,14 +426,14 @@ function compute_flow_balance!(
     qy_av_average = qy_average
 
     for i in 1:(overland_flow_model.parameters.n)
-        idx_down = indices.idx_down[i]
-        idx_left = indices.idx_left[i]
+        ind_y_down = indices.ind_y_down[i]
+        ind_x_down = indices.ind_x_down[i]
         total_in = 0.0
         total_out = 0.0
         total_in, total_out = add_inflow(
             total_in,
             total_out,
-            [qx_av_average[idx_left], qy_av_average[idx_down], runoff[i]],
+            [qx_av_average[ind_x_down], qy_av_average[ind_y_down], runoff[i]],
         )
         total_in, total_out =
             add_outflow(total_in, total_out, [qx_av_average[i], qy_av_average[i]])
@@ -512,7 +512,9 @@ end
 
 function compute_flow_routing_balance!(
     model::Model{R},
-) where {R <: Routing{<:LocalInertialOverlandFlowModel, <:LocalInertialRiverFlowModel}}
+) where {
+    R <: Routing{<:OverlandFlowModel{<:LocalInertial}, <:RiverFlowModel{<:LocalInertial}},
+}
     (; river_flow, overland_flow, subsurface_flow) = model.routing
     (; reservoir) = river_flow.boundary_conditions
     (; overland_water_balance, reservoir_water_balance, subsurface_water_balance) =
diff --git a/Wflow/src/network.jl b/Wflow/src/network.jl
index b8cdac242..be8260b03 100644
--- a/Wflow/src/network.jl
+++ b/Wflow/src/network.jl
@@ -1,6 +1,6 @@
 # maps the fields of struct `EdgeConnectivity` to the defined Wflow cartesian indices of
 # const `neighbors`.
-const DIRS = (:idx_down, :idx_left, :idx_right, :idx_up)
+const DIRS = (:ind_y_down, :ind_x_down, :ind_x_up, :ind_y_up)
 
 """
 Struct for storing 2D staggered grid edge connectivity in `x` and `y` directions. For
@@ -15,19 +15,21 @@ Edges without neighbors are handled by an extra index (at `n + 1`, with `n` edge
 linear index `i` of the `EdgeConnectivity` fields represents the edge between node index `i`
 and the neighboring nodes in the CartesianIndex(-1,0) and CartesianIndex(0,-1) directions.
 The edges are defined as follows:
-- `idx_right` is the edge between node `i` and node `idx_right` in the `CartesianIndex(1,0)` direction.
-- `idx_left` is the edge between node `idx_left` in the `CartesianIndex(-1,0)` direction and the
-  neighboring node (CartesianIndex(-2,0) direction).
-- `idx_up` is the edge between node `i` and node `idx_up` in the `CartesianIndex(0,1)` direction.
-- `idx_down` is the edge between node `idx_down` in the `CartesianIndex(0,-1)` direction and the
-  neighboring node (`CartesianIndex(0,-2)` direction).
+- `ind_x_up` is the edge between node `i` and node `ind_x_up` in the `CartesianIndex(1,0)`
+  direction.
+- `ind_x_down` is the edge between node `ind_x_down` in the `CartesianIndex(-1,0)` direction
+  and the neighboring node (CartesianIndex(-2,0) direction).
+- `ind_y_up` is the edge between node `i` and node `ind_y_up` in the `CartesianIndex(0,1)`
+  direction.
+- `ind_y_down` is the edge between node `ind_y_down` in the `CartesianIndex(0,-1)` direction
+  and the neighboring node (`CartesianIndex(0,-2)` direction).
 """
 @with_kw struct EdgeConnectivity
     n::Int
-    idx_right::Vector{Int} = zeros(Int, n)
-    idx_left::Vector{Int} = zeros(Int, n)
-    idx_up::Vector{Int} = zeros(Int, n)
-    idx_down::Vector{Int} = zeros(Int, n)
+    ind_x_up::Vector{Int} = zeros(Int, n)
+    ind_x_down::Vector{Int} = zeros(Int, n)
+    ind_y_up::Vector{Int} = zeros(Int, n)
+    ind_y_down::Vector{Int} = zeros(Int, n)
 end
 
 "Struct for storing source `src` node and destination `dst` node of an edge."
@@ -115,7 +117,7 @@ nthreads > 1 to run the kinematic wave parallel, otherwise it is equal to the co
 domain.
 """
 function network_subdomains(config::Config, network::NetworkLand)
-    pit_inds = findall(x -> x == 5, network.local_drain_direction)
+    pit_inds = findall(x -> x == LDD_PIT, network.local_drain_direction)
     order_of_subdomains, subdomain_inds, toposort_subdomain = kinwave_set_subdomains(
         network.graph,
         network.order,
@@ -261,7 +263,7 @@ nthreads > 1 to run the kinematic wave parallel, otherwise it is equal to the co
 domain.
 """
 function network_subdomains(config::Config, network::NetworkRiver)
-    pit_inds = findall(x -> x == 5, network.local_drain_direction)
+    pit_inds = findall(x -> x == LDD_PIT, network.local_drain_direction)
     order_of_subdomains, subdomain_inds, toposort_subdomain = kinwave_set_subdomains(
         network.graph,
         network.order,
@@ -277,7 +279,7 @@ end
 
 "Initialize `NodesAtEdge`"
 function NodesAtEdge(network::NetworkRiver)
-    index_pit = findall(x -> x == 5, network.local_drain_direction)
+    index_pit = findall(x -> x == LDD_PIT, network.local_drain_direction)
     add_vertex_edge_graph!(network.graph, index_pit)
     nodes_at_edge = NodesAtEdge(; adjacent_nodes_at_edge(network.graph)...)
     return nodes_at_edge, index_pit
diff --git a/Wflow/src/routing/initialize_routing.jl b/Wflow/src/routing/initialize_routing.jl
index 135d5ca8f..98f515d71 100644
--- a/Wflow/src/routing/initialize_routing.jl
+++ b/Wflow/src/routing/initialize_routing.jl
@@ -52,9 +52,9 @@ function initialize_overland_flow(dataset::NCDataset, config::Config, domain::Do
     (; land_routing) = config.model
 
     if land_routing == RoutingType.kinematic_wave
-        overland_flow = KinWaveOverlandFlowModel(dataset, config, domain.land)
+        overland_flow = init_kinematic_wave_overland_flow(dataset, config, domain.land)
     elseif land_routing == RoutingType.local_inertial
-        overland_flow = LocalInertialOverlandFlowModel(dataset, config, domain)
+        overland_flow = init_local_inertial_overland_flow(dataset, config, domain)
     end
     return overland_flow
 end
@@ -71,9 +71,11 @@ function initialize_river_flow(dataset::NCDataset, config::Config, domain::Domai
         ReservoirModel(dataset, config, domain.reservoir.network) : nothing
 
     if river_routing == RoutingType.kinematic_wave
-        river_flow = KinWaveRiverFlowModel(dataset, config, domain.river, reservoir)
-    elseif river_routing == RoutingType.local_inertial
-        river_flow = LocalInertialRiverFlowModel(dataset, config, domain.river, reservoir)
+        river_flow =
+            init_kinematic_wave_river_flow(dataset, config, domain.river, reservoir)
+    elseif river_routing == RoutingType.local_inertial ||
+           river_routing == RoutingType.manning_staggered
+        river_flow = init_staggered_river_flow(dataset, config, domain.river, reservoir)
     end
 end
 
diff --git a/Wflow/src/routing/routing.jl b/Wflow/src/routing/routing.jl
index 0259fd524..1ae8c710c 100644
--- a/Wflow/src/routing/routing.jl
+++ b/Wflow/src/routing/routing.jl
@@ -1,6 +1,14 @@
-abstract type AbstractSubsurfaceFlowModel end
-abstract type AbstractOverlandFlowModel end
 abstract type AbstractRiverFlowModel end
+abstract type AbstractOverlandFlowModel end
+abstract type AbstractSubsurfaceFlowModel end
+
+abstract type AbstractRoutingMethod end
+abstract type AbstractStaggeredRoutingMethod <: AbstractRoutingMethod end
+
+struct Manning <: AbstractRoutingMethod end
+struct KinematicWave <: AbstractRoutingMethod end
+struct ManningStaggered <: AbstractStaggeredRoutingMethod end
+struct LocalInertial <: AbstractStaggeredRoutingMethod end
 
 struct NoSubsurfaceFlow <: AbstractSubsurfaceFlowModel end
 struct NoOverlandFlow <: AbstractOverlandFlowModel end
diff --git a/Wflow/src/routing/surface/floodplain.jl b/Wflow/src/routing/surface/floodplain.jl
new file mode 100644
index 000000000..bf1f527a2
--- /dev/null
+++ b/Wflow/src/routing/surface/floodplain.jl
@@ -0,0 +1,402 @@
+abstract type AbstractFloodPlainModel end
+abstract type AbstractFloodPlainParameters end
+abstract type AbstractFloodPlainVariables end
+
+""""
+    FloodPlainModel
+
+Floodplain flow model for each river cell (as part of a `RiverFlowModel`). The
+`FloodPlainModel` can be part of three different river routing schemes:
+* Local inertial equation on a staggered grid.
+* Kinematic wave using Manning's equation on a staggered grid.
+* Kinematic wave approach by solving the kinematic wave equation using Newton's method.
+
+On a staggered grid the river and floodplain flow routing scheme are equal. For the
+kinematic wave river flow routing solved by using Newton's method, the floodplain flow
+routing is solved by using Manning's equation, as solving kinematic wave routing for near
+zero flows is computationally expensive.
+"""
+@with_kw struct FloodPlainModel{
+    T <: AbstractRoutingMethod,
+    P <: AbstractFloodPlainParameters,
+    V <: AbstractFloodPlainVariables,
+} <: AbstractFloodPlainModel
+    routing_method::T
+    parameters::P
+    variables::V
+end
+
+"""
+    FloodPlainProfile
+
+Floodplain `storage` is a function of `depth` (flood depth intervals). Based on the
+cumulative floodplain `storage` a floodplain profile as a function of `flood_depth` is
+derived with floodplain area `a` (cumulative) and wetted perimeter radius `p` (cumulative).
+"""
+@with_kw struct FloodPlainProfile
+    # Flood depth [m]
+    depth::Vector{Float64}
+    # Flood storage (cumulative) [m³]
+    storage::Matrix{Float64}
+    # Flood width [m]
+    width::Matrix{Float64}
+    # Flow area (cumulative) [m²]
+    flow_area::Matrix{Float64}
+    # Wetted perimeter (cumulative) [m]
+    wetted_perimeter::Matrix{Float64}
+end
+
+"Initialize floodplain profile `FloodPlainProfile`"
+function FloodPlainProfile(
+    dataset::NCDataset,
+    config::Config,
+    domain::DomainRiver;
+    index_pit::Vector{Int} = Int[],
+)
+    (; indices) = domain.network
+    (; flow_width, flow_length) = domain.parameters
+    storage = ncread(
+        dataset,
+        config,
+        "floodplain_water__sum_of_volume_per_depth",
+        Routing;
+        sel = indices,
+    )
+    n = length(indices)
+
+    # for convenience (interpolation) flood depth 0.0 m is added, with associated area (a),
+    # storage, width (river width) and wetted perimeter (p).
+    storage = vcat(fill(Float64(0), n)', storage)
+    start_storage = storage
+    flood_depths = Float64.(dataset["flood_depth"][:])
+    pushfirst!(flood_depths, 0.0)
+    n_depths = length(flood_depths)
+
+    wetted_perimeter = zeros(n_depths, n)
+    flow_area = zeros(n_depths, n)
+    segment_storage = zeros(n_depths, n)
+    width = zeros(n_depths, n)
+    width[1, :] = flow_width[1:n]
+
+    # determine flow area (a), width and wetted perimeter (p) FloodPlainModel
+    h = diff(flood_depths)
+    incorrect_vol = 0
+    riv_cells = 0
+    error_vol = 0
+    for i in 1:n
+        riv_cell = 0
+        diff_storage = diff(storage[:, i])
+
+        for j in 1:(n_depths - 1)
+            # assume rectangular shape of flood depth segment
+            width[j + 1, i] = diff_storage[j] / (h[j] * flow_length[i])
+            # check provided flood storage (floodplain width should be constant or increasing
+            # as a function of flood depth)
+            if width[j + 1, i] < width[j, i]
+                # raise warning only if difference is larger than rounding error of 0.01 m³
+                if ((width[j, i] - width[j + 1, i]) * h[j] * flow_length[i]) > 0.01
+                    incorrect_vol += 1
+                    riv_cell = 1
+                    error_vol =
+                        error_vol +
+                        ((width[j, i] - width[j + 1, i]) * h[j] * flow_length[i])
+                end
+                width[j + 1, i] = width[j, i]
+            end
+            flow_area[j + 1, i] = width[j + 1, i] * h[j]
+            wetted_perimeter[j + 1, i] = (width[j + 1, i] - width[j, i]) + 2.0 * h[j]
+            segment_storage[j + 1, i] = flow_area[j + 1, i] * flow_length[i]
+            if j == 1
+                # for interpolation wetted perimeter at flood depth 0.0 is required
+                wetted_perimeter[j, i] = wetted_perimeter[j + 1, i] - 2.0 * h[j]
+            end
+        end
+
+        wetted_perimeter[2:end, i] = cumsum(wetted_perimeter[2:end, i])
+        flow_area[:, i] = cumsum(flow_area[:, i])
+        storage[:, i] = cumsum(segment_storage[:, i])
+
+        riv_cells += riv_cell
+    end
+
+    if incorrect_vol > 0
+        perc_riv_cells = round(100.0 * (riv_cells / n); digits = 2)
+        perc_error_vol = round(100.0 * (error_vol / sum(start_storage[end, :])); digits = 2)
+        @warn string(
+            "The provided storage of $incorrect_vol rectangular floodplain schematization",
+            " segments for $riv_cells river cells ($perc_riv_cells % of total river cells)",
+            " is not correct and has been increased with $perc_error_vol % of provided storage.",
+        )
+    end
+
+    # set floodplain parameters for ghost points
+    storage = hcat(storage, storage[:, index_pit])
+    width = hcat(width, width[:, index_pit])
+    flow_area = hcat(flow_area, flow_area[:, index_pit])
+    wetted_perimeter = hcat(wetted_perimeter, wetted_perimeter[:, index_pit])
+
+    # initialize floodplain profile parameters
+    profile = FloodPlainProfile(;
+        storage,
+        width,
+        depth = flood_depths,
+        flow_area,
+        wetted_perimeter,
+    )
+    return profile
+end
+
+"Struct to store floodplain flow model parameters on a staggered grid"
+@with_kw struct FloodPlainStaggeredParameters <: AbstractFloodPlainParameters
+    # floodplain profile
+    profile::FloodPlainProfile
+    # manning's roughness [s m-1/3]
+    mannings_n::Vector{Float64} = Float64[]
+    # manning's roughness at edge [s m-1/3]
+    mannings_n_at_edge::Vector{Float64} = Float64[]
+    # manning's roughness squared at edge [(s m-1/3)2]
+    mannings_n_sq_at_edge::Vector{Float64} = Float64[]
+    # bankfull elevation at edge [m]
+    zb_at_edge::Vector{Float64} = Float64[]
+    # slope at edge [-]
+    slope_at_edge::Vector{Float64} = Float64[]
+end
+
+"Initialize floodplain flow model parameters on a staggered grid"
+function FloodPlainStaggeredParameters(
+    dataset::NCDataset,
+    config::Config,
+    domain::DomainRiver,
+    zb_floodplain::Vector{Float64},
+    index_pit::Vector{Int},
+)
+    (; indices, nodes_at_edge, graph) = domain.network
+    (; flow_length) = domain.parameters
+    (; river_routing) = config.model
+    n_edges = ne(graph)
+    profile = FloodPlainProfile(dataset, config, domain; index_pit)
+    mannings_n = ncread(
+        dataset,
+        config,
+        "floodplain_water_flow__manning_n_parameter",
+        Routing;
+        sel = indices,
+    )
+    # manning roughness at edges
+    append!(mannings_n, mannings_n[index_pit]) # copy to ghost nodes
+    mannings_n_at_edge =
+        compute_mannings_n_at_edge(mannings_n, flow_length, nodes_at_edge, n_edges)
+
+    if river_routing == RoutingType.local_inertial
+        zb_at_edge = compute_value_at_edge(zb_floodplain, nodes_at_edge, n_edges, maximum)
+        mannings_n_sq_at_edge = mannings_n_at_edge .* mannings_n_at_edge
+        slope_at_edge = []
+    elseif river_routing == RoutingType.manning_staggered
+        zb_at_edge = compute_value_at_edge(zb_floodplain, nodes_at_edge, n_edges, first)
+        mannings_n_sq_at_edge = []
+        flow_length_at_edge =
+            compute_value_at_edge(flow_length, nodes_at_edge, n_edges, mean)
+        slope_at_edge = compute_slope_at_edge(
+            zb_floodplain,
+            flow_length_at_edge,
+            nodes_at_edge,
+            n_edges,
+        )
+    end
+
+    parameters = FloodPlainStaggeredParameters(;
+        profile,
+        mannings_n,
+        mannings_n_at_edge,
+        mannings_n_sq_at_edge,
+        zb_at_edge,
+        slope_at_edge,
+    )
+    return parameters
+end
+
+"Struct to store floodplain flow model variables on a staggered grid"
+@with_kw struct FloodPlainStaggeredVariables <: AbstractFloodPlainVariables
+    n::Int
+    n_edges::Int
+    # water depth at edge [m]
+    water_depth_at_edge::Vector{Float64} = zeros(n_edges)
+    # discharge at edge at previous time step
+    q_previous::Vector{Float64} = zeros(n_edges)
+    # discharge at edge  [m³ s⁻¹]
+    q::Vector{Float64} = zeros(n_edges)
+    # cumulative discharge at edge [m³] for model timestep dt
+    q_cumulative::Vector{Float64} = zeros(n_edges)
+    # average discharge at edge [m³ s⁻¹] for model timestep dt
+    q_average::Vector{Float64} = zeros(n_edges)
+    # storage [m³]
+    storage::Vector{Float64} = zeros(n)
+    # water depth [m]
+    h::Vector{Float64}
+    # error storage [m³]
+    error::Vector{Float64} = zeros(n)
+end
+
+"Struct to store floodplain parameters"
+@with_kw struct FloodPlainParameters <: AbstractFloodPlainParameters
+    # floodplain profile
+    profile::FloodPlainProfile
+    # manning's roughness[s m-1/3]
+    mannings_n::Vector{Float64} = Float64[]
+    # slope  [m m⁻¹]
+    slope::Vector{Float64} = Float64[]
+end
+
+"Struct to store floodplain variables"
+@with_kw struct FloodPlainVariables <: AbstractFloodPlainVariables
+    n::Int
+    # discharge [m³ s⁻¹]
+    q::Vector{Float64} = zeros(n)
+    # cumulative discharge for model timestep dt
+    q_cumulative::Vector{Float64} = zeros(n)
+    # average discharge [m³ s⁻¹] for model timestep Δt
+    q_average::Vector{Float64} = zeros(n)
+    # inflow from upstream cells [m³ s⁻¹]
+    qin::Vector{Float64} = zeros(n)
+    # cumulative inflow from upstream cells for model timestep dt
+    qin_cumulative::Vector{Float64} = zeros(n)
+    # average inflow from upstream cells [m³ s⁻¹] for model timestep Δt
+    qin_average::Vector{Float64} = zeros(n)
+    # storage [m³]
+    storage::Vector{Float64} = zeros(n)
+    # water depth [m]
+    h::Vector{Float64} = zeros(n)
+    # flow capacity [m³ s⁻¹]
+    flow_capacity::Vector{Float64} = zeros(n)
+end
+
+"Determine the initial floodplain storage"
+function initialize_storage!(river, domain::Domain, nriv::Int)
+    (; flow_width, flow_length) = domain.river.parameters
+    (; floodplain) = river
+    (; profile) = floodplain.parameters
+    for i in 1:nriv
+        i1, i2 = interpolation_indices(floodplain.variables.h[i], profile.depth)
+        a = compute_floodplain_flow_area(profile, floodplain.variables.h[i], i, i1, i2)
+        floodplain.variables.storage[i] = flow_length[i] * a
+    end
+    return nothing
+end
+
+"helper function to get interpolation indices"
+function interpolation_indices(x, v::AbstractVector)
+    i1 = 1
+    for i in eachindex(v)
+        if v[i] <= x
+            i1 = i
+        end
+    end
+    if i1 == length(v)
+        i2 = i1
+    else
+        i2 = i1 + 1
+    end
+    return i1, i2
+end
+
+"""
+Compute flood flow area (including area above channel) based on flow depth `h` and
+floodplain `depth`, `flow_area` and `width` of a floodplain profile.
+"""
+function compute_flood_flow_area(
+    profile::FloodPlainProfile,
+    h::Float64,
+    idx::Int,
+    i1::Int,
+    i2::Int,
+)
+    delta_h = h - profile.depth[i1]  # depth at i1
+    flow_area = profile.flow_area[i1, idx] + (profile.width[i2, idx] * delta_h) # area at i1, width at i2
+    return flow_area
+end
+
+"""
+Compute floodplain wetted perimeter based on flow depth `h` and floodplain `depth` and
+`wetted_perimeter` of a floodplain profile.
+"""
+function compute_wetted_perimeter(profile::FloodPlainProfile, h::Float64, idx::Int, i1::Int)
+    delta_h = h - profile.depth[i1] # depth at i1
+    wetted_perimeter = profile.wetted_perimeter[i1, idx] + 2.0 * delta_h # p at i1
+    return wetted_perimeter
+end
+
+"Compute flood depth by interpolating flood storage `flood_storage` using flood depth intervals."
+function compute_flood_depth(
+    profile::FloodPlainProfile,
+    flood_storage::Float64,
+    flow_length::Float64,
+    i::Int,
+)
+    i1, i2 = interpolation_indices(flood_storage, @view profile.storage[:, i])
+    delta_A = (flood_storage - profile.storage[i1, i]) / flow_length
+    delta_h = delta_A / profile.width[i2, i]
+    flood_depth = profile.depth[i1] + delta_h
+    return flood_depth
+end
+
+"Compute floodplain flow area (excluding river channel area)"
+function compute_floodplain_flow_area(
+    profile::FloodPlainProfile,
+    h::Float64,
+    idx::Int,
+    i1::Int,
+    i2::Int,
+)
+    channel_area = profile.width[1, idx] * h
+    floodplain_flow_area = compute_flood_flow_area(profile, h, idx, i1, i2)
+    floodplain_flow_area = max(floodplain_flow_area - channel_area, 0.0)
+    return floodplain_flow_area
+end
+
+"Initialize floodplain geometry, model variables and parameters on staggered grid"
+function FloodPlainModel(
+    dataset::NCDataset,
+    config::Config,
+    domain::DomainRiver,
+    zb_floodplain::Vector{Float64},
+)
+    (; indices, local_drain_direction, graph) = domain.network
+    (; river_routing) = config.model
+    n = length(indices)
+    index_pit = findall(x -> x == LDD_PIT, local_drain_direction)
+    n_edges = ne(graph)
+    parameters =
+        FloodPlainStaggeredParameters(dataset, config, domain, zb_floodplain, index_pit)
+    if river_routing == RoutingType.local_inertial
+        routing_method = LocalInertial()
+    elseif river_routing == RoutingType.manning_staggered
+        routing_method = ManningStaggered()
+    end
+    variables = FloodPlainStaggeredVariables(; n, n_edges, h = zeros(n + length(index_pit)))
+
+    floodplain_model = FloodPlainModel(; routing_method, parameters, variables)
+    return floodplain_model
+end
+
+"""
+Initialize floodplain geometry, model variables and parameters for floodplain flow routing
+as part of kinematic wave river flow routing (solved using Newton's method).
+"""
+function FloodPlainModel(dataset::NCDataset, config::Config, domain::DomainRiver)
+    (; indices) = domain.network
+    n = length(indices)
+    profile = FloodPlainProfile(dataset, config, domain)
+    mannings_n = ncread(
+        dataset,
+        config,
+        "floodplain_water_flow__manning_n_parameter",
+        Routing;
+        sel = indices,
+    )
+    slope = ncread(dataset, config, "floodplain__slope", Routing; sel = indices)
+    clamp!(slope, 0.00001, Inf)
+    parameters = FloodPlainParameters(; profile, mannings_n, slope)
+    variables = FloodPlainVariables(; n)
+    floodplain_model = FloodPlainModel(; routing_method = Manning(), parameters, variables)
+    return floodplain_model
+end
diff --git a/Wflow/src/routing/surface/surface_flow.jl b/Wflow/src/routing/surface/surface_flow.jl
new file mode 100644
index 000000000..f9f25a8d8
--- /dev/null
+++ b/Wflow/src/routing/surface/surface_flow.jl
@@ -0,0 +1,71 @@
+abstract type AbstractRiverFlowParameters end
+abstract type AbstractRiverFlowVariables end
+
+abstract type AbstractOverlandFlowBC end
+abstract type AbstractOverlandFlowParameters end
+abstract type AbstractOverlandFlowVariables end
+
+"Struct for storing river flow model boundary conditions"
+@with_kw struct RiverFlowBC{R <: Union{ReservoirModel, Nothing}}
+    n::Int
+    # External inflow (abstraction/supply/demand) [m³ s⁻¹]
+    external_inflow::Vector{Float64} = zeros(n)
+    # Lateral inflow [m³ s⁻¹]
+    inwater::Vector{Float64} = zeros(n)
+    # Actual cumulative abstraction from external negative inflow [m³]
+    actual_external_abstraction_cumulative::Vector{Float64} = zeros(n)
+    # Actual average abstraction from external negative inflow [m³ s⁻¹]
+    actual_external_abstraction_average::Vector{Float64} = zeros(n)
+    # Abstraction (computed as part of water demand and allocation) [m³ s⁻¹]
+    abstraction::Vector{Float64} = zeros(n)
+    # Exchange with floodplain [m³ s⁻¹]
+    floodplain_water_exchange::Vector{Float64} = zeros(n)
+    # Reservoir model struct of arrays
+    reservoir::R
+end
+
+"Struct for storing Manning flow parameters"
+@with_kw struct ManningFlowParameters
+    slope::Vector{Float64}      # Slope [m m⁻¹]
+    mannings_n::Vector{Float64} # Manning's roughness [s m⁻⅓]
+    alpha_pow::Float64          # Used in the power part of alpha [-]
+    alpha_term::Vector{Float64} # Term used in computation of alpha [-]
+    alpha::Vector{Float64}      # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation [s3/5 m1/5]
+end
+
+"River flow model"
+@with_kw struct RiverFlowModel{
+    T <: AbstractRoutingMethod,
+    F <: Union{AbstractFloodPlainModel, Nothing},
+    P <: AbstractRiverFlowParameters,
+    V <: AbstractRiverFlowVariables,
+    A <: AbstractAllocationModel,
+} <: AbstractRiverFlowModel
+    routing_method::T
+    timestepping::TimeStepping
+    boundary_conditions::RiverFlowBC
+    parameters::P
+    variables::V
+    floodplain::F
+    allocation::A
+end
+
+"Overland flow model"
+@with_kw struct OverlandFlowModel{
+    T <: AbstractRoutingMethod,
+    B <: AbstractOverlandFlowBC,
+    P <: Union{ManningFlowParameters, AbstractOverlandFlowParameters},
+    V <: AbstractOverlandFlowVariables,
+} <: AbstractOverlandFlowModel
+    routing_method::T
+    timestepping::TimeStepping
+    boundary_conditions::B
+    parameters::P
+    variables::V
+end
+
+"Compute wetted perimeter for a rectangular river channel "
+function wetted_perimeter_channel(h::Float64, flow_width::Float64)
+    channel_perimeter = 2.0 * h + flow_width
+    return channel_perimeter
+end
diff --git a/Wflow/src/routing/surface/surface_kinwave.jl b/Wflow/src/routing/surface/surface_kinwave.jl
index ff561c466..0e21e426b 100644
--- a/Wflow/src/routing/surface/surface_kinwave.jl
+++ b/Wflow/src/routing/surface/surface_kinwave.jl
@@ -1,5 +1,8 @@
-"Struct for storing (shared) variables for river and overland flow models"
-@with_kw struct FlowVariables
+# constant in Manning's equation [-]
+const BETA_KINWAVE = 0.6
+
+"Struct for storing variables for river flow model"
+@with_kw struct RiverFlowVariables <: AbstractRiverFlowVariables
     n::Int
     # Discharge [m³ s⁻¹]
     q::Vector{Float64} = zeros(n)
@@ -11,40 +14,48 @@
     qin_cumulative::Vector{Float64} = zeros(n)
     # Average inflow from upstream cells  [m³ s⁻¹] for model time step dt
     qin_average::Vector{Float64} = zeros(n)
-    # Cumulative discharge [m³] for model timestep dt
+    # Cumulative river channel (+ floodplain) discharge [m³] for model timestep dt
     q_cumulative::Vector{Float64} = zeros(n)
-    # Average discharge [m³ s⁻¹] for model time step dt
+    # Average river channel (+ floodplain) discharge [m³ s⁻¹] for model time step dt
     q_average::Vector{Float64} = zeros(n)
+    # Cumulative river channel discharge [m³] (for model time step dt)
+    q_channel_cumulative::Vector{Float64} = q_cumulative
+    # Average river channel discharge [m³ s⁻¹] (for model time step dt)
+    q_channel_average::Vector{Float64} = q_average
     # Kinematic wave storage [m³] (based on water depth h)
     storage::Vector{Float64} = zeros(n)
     # Water depth [m]
     h::Vector{Float64} = zeros(n)
 end
 
-"Struct for storing Manning flow parameters"
-@with_kw struct ManningFlowParameters
-    # Slope [m m⁻¹]
-    slope::Vector{Float64}
-    # Manning's roughness [s m⁻⅓]
-    mannings_n::Vector{Float64}
-    # Used in the power part of alpha [-]
-    alpha_pow::Float64 = (2 // 3) * 0.6
-    # Term used in computation of alpha [s^3/5 m^-1/5]
-    alpha_term::Vector{Float64} = fill(MISSING_VALUE, length(slope))
-    # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation [s³ᐟ⁵ m¹ᐟ⁵]
-    alpha::Vector{Float64} = fill(MISSING_VALUE, length(slope))
+"Initialize Manning flow parameters"
+function ManningFlowParameters(
+    mannings_n::Vector{Float64},
+    slope::Vector{Float64},
+    wetted_perimeter::Vector{Float64},
+)
+    hydraulic_radius_pow = 2.0 / 3.0
+    alpha_term = @. pow(mannings_n / sqrt(slope), BETA_KINWAVE)
+    alpha_pow = hydraulic_radius_pow * BETA_KINWAVE
+    alpha = @. alpha_term * pow(wetted_perimeter, alpha_pow)
+
+    parameters = ManningFlowParameters(; slope, mannings_n, alpha_pow, alpha_term, alpha)
+    return parameters
 end
 
 "Struct for storing river flow model parameters"
-@with_kw struct RiverFlowParameters
+@with_kw struct RiverFlowParameters <: AbstractRiverFlowParameters
     flow::ManningFlowParameters
-    bankfull_depth::Vector{Float64} # Bankfull water level [m]
+    bankfull_depth::Vector{Float64}     # Bankfull water level [m]
+    bankfull_storage::Vector{Float64}   # Bankfull storage [m³]
 end
 
 "Overload `getproperty` for river flow model parameters"
 function Base.getproperty(v::RiverFlowParameters, s::Symbol)
     if s === :bankfull_depth
         getfield(v, s)
+    elseif s === :bankfull_storage
+        getfield(v, s)
     elseif s === :flow
         getfield(v, :flow)
     else
@@ -55,7 +66,7 @@ end
 "Initialize river flow model parameters"
 function RiverFlowParameters(dataset::NCDataset, config::Config, domain::DomainRiver)
     (; indices) = domain.network
-    (; slope) = domain.parameters
+    (; slope, flow_length, flow_width) = domain.parameters
     mannings_n = ncread(
         dataset,
         config,
@@ -66,28 +77,14 @@ function RiverFlowParameters(dataset::NCDataset, config::Config, domain::DomainR
     bankfull_depth =
         ncread(dataset, config, "river_bank_water__depth", Routing; sel = indices)
 
-    flow_params = ManningFlowParameters(; mannings_n, slope)
-    parameters = RiverFlowParameters(; flow = flow_params, bankfull_depth)
+    # use fixed wetted perimeter based on 0.5 * bankfull_depth
+    wetted_perimeter = flow_width + bankfull_depth
+    flow_params = ManningFlowParameters(mannings_n, slope, wetted_perimeter)
+    bankfull_storage = bankfull_depth .* flow_width .* flow_length
+    parameters = RiverFlowParameters(; flow = flow_params, bankfull_depth, bankfull_storage)
     return parameters
 end
 
-"Struct for storing river flow model boundary conditions"
-@with_kw struct RiverFlowBC{R <: Union{ReservoirModel, Nothing}}
-    n::Int
-    # External inflow (abstraction/supply/demand) [m³ s⁻¹]
-    external_inflow::Vector{Float64}
-    # Lateral inflow [m³ s⁻¹]
-    inwater::Vector{Float64} = zeros(n)
-    # Actual cumulative abstraction from external negative inflow [m³]
-    actual_external_abstraction_cumulative::Vector{Float64} = zeros(n)
-    # Actual average abstraction from external negative inflow [m³ s⁻¹]
-    actual_external_abstraction_average::Vector{Float64} = zeros(n)
-    # Abstraction (computed as part of water demand and allocation) [m³ s⁻¹]
-    abstraction::Vector{Float64} = zeros(n)
-    # Reservoir model struct of arrays
-    reservoir::R
-end
-
 "Initialize river flow model boundary conditions"
 function RiverFlowBC(
     dataset::NCDataset,
@@ -108,19 +105,8 @@ function RiverFlowBC(
     return bc
 end
 
-"River flow model using the kinematic wave method and the Manning flow equation"
-@with_kw struct KinWaveRiverFlowModel{R <: RiverFlowBC, A <: AbstractAllocationModel} <:
-                AbstractRiverFlowModel
-    timestepping::TimeStepping
-    boundary_conditions::R
-    parameters::RiverFlowParameters
-    variables::FlowVariables
-    # Water allocation
-    allocation::A
-end
-
-"Initialize river flow model `KinWaveRiverFlowModel`"
-function KinWaveRiverFlowModel(
+"Initialize kinematic wave river flow model"
+function init_kinematic_wave_river_flow(
     dataset::NCDataset,
     config::Config,
     domain::DomainRiver,
@@ -134,60 +120,78 @@ function KinWaveRiverFlowModel(
     allocation =
         do_water_demand(config) ? AllocationRiverModel(; n) : NoAllocationRiverModel(n)
 
-    variables = FlowVariables(; n)
+    variables = if config.model.floodplain_1d__flag
+        # When floodplain is enabled, q_channel fields must be separate from q fields
+        # because q_average will hold channel + floodplain combined discharge.
+        RiverFlowVariables(; n, q_channel_cumulative = zeros(n), q_channel_average = zeros(n))
+    else
+        # When floodplain is disabled, q_channel aliases q (they are identical)
+        RiverFlowVariables(; n)
+    end
     parameters = RiverFlowParameters(dataset, config, domain)
     boundary_conditions = RiverFlowBC(dataset, config, domain.network, reservoir)
+    if config.model.floodplain_1d__flag
+        floodplain = FloodPlainModel(dataset, config, domain)
+    else
+        floodplain = nothing
+    end
+    routing_method = KinematicWave()
 
-    river_flow = KinWaveRiverFlowModel(;
+    river_flow = RiverFlowModel(;
         timestepping,
         boundary_conditions,
         parameters,
         variables,
+        floodplain,
         allocation,
+        routing_method,
     )
 
     return river_flow
 end
 
 "Struct for storing overland flow model variables"
-@with_kw struct OverLandFlowVariables
+@with_kw struct OverLandFlowVariables <: AbstractOverlandFlowVariables
     n::Int
-    flow::FlowVariables = FlowVariables(; n)
+    # Overland discharge [m³ s⁻¹]
+    q::Vector{Float64} = zeros(n)
+    # Lateral inflow per unit length [m² s⁻¹]
+    qlat::Vector{Float64} = zeros(n)
+    # Overland inflow from upstream cells [m³ s⁻¹]
+    qin::Vector{Float64} = zeros(n)
+    # Cumulative overland inflow from upstream cells [m³] for model time step dt
+    qin_cumulative::Vector{Float64} = zeros(n)
+    # Average overland inflow from upstream cells [m³ s⁻¹] for model timestep Δt
+    qin_average::Vector{Float64} = zeros(n)
+    # Cumulative overland discharge [m³] for model timestep dt
+    q_cumulative::Vector{Float64} = zeros(n)
+    # Average overland discharge [m³ s⁻¹] for model timestep Δt
+    q_average::Vector{Float64} = zeros(n)
+    # Overland kinematic wave storage [m³] (based on water depth h)
+    storage::Vector{Float64} = zeros(n)
+    # Overland water depth [m]
+    h::Vector{Float64} = zeros(n)
     # Part of cumulative overland flow [m³ s⁻¹] that flows to the river
     to_river_cumulative::Vector{Float64} = zeros(n)
     # Part of average overland flow [m³ s⁻¹] that flows to the river
     to_river_average::Vector{Float64} = zeros(n)
 end
 
-"Overload `getproperty` for overland flow model variables"
-function Base.getproperty(v::OverLandFlowVariables, s::Symbol)
-    if s ∈ (:to_river_average, :to_river_cumulative, :flow)
-        getfield(v, s)
-    else
-        getfield(getfield(v, :flow), s)
-    end
-end
-
 "Struct for storing overland flow model boundary conditions"
-@with_kw struct LandFlowBC
+@with_kw struct LandFlowBC <: AbstractOverlandFlowBC
     n::Int
     # Lateral inflow [m³ s⁻¹]
     inwater::Vector{Float64} = zeros(n)
 end
 
-"Overland flow model using the kinematic wave method and the Manning flow{ equation"
-@with_kw struct KinWaveOverlandFlowModel <: AbstractOverlandFlowModel
-    n::Int
-    timestepping::TimeStepping
-    boundary_conditions::LandFlowBC = LandFlowBC(; n)
-    parameters::ManningFlowParameters
-    variables::OverLandFlowVariables = OverLandFlowVariables(; n)
-end
-
-"Initialize Overland flow model `KinWaveOverlandFlowModel`"
-function KinWaveOverlandFlowModel(dataset::NCDataset, config::Config, domain::DomainLand)
+"Initialize kinematic wave overland flow model"
+function init_kinematic_wave_overland_flow(
+    dataset::NCDataset,
+    config::Config,
+    domain::DomainLand,
+)
     (; indices) = domain.network
-    (; slope) = domain.parameters
+    (; slope, surface_flow_width) = domain.parameters
     mannings_n = ncread(
         dataset,
         config,
@@ -198,8 +202,19 @@ function KinWaveOverlandFlowModel(dataset::NCDataset, config::Config, domain::Do
 
     n = length(indices)
     timestepping = init_kinematic_wave_timestepping(config, n; domain = "land")
-    parameters = ManningFlowParameters(; mannings_n, slope)
-    overland_flow = KinWaveOverlandFlowModel(; n, timestepping, parameters)
+
+    variables = OverLandFlowVariables(; n)
+    parameters = ManningFlowParameters(mannings_n, slope, surface_flow_width)
+    boundary_conditions = LandFlowBC(; n)
+    routing_method = KinematicWave()
+
+    overland_flow = OverlandFlowModel(;
+        timestepping,
+        boundary_conditions,
+        variables,
+        parameters,
+        routing_method,
+    )
 
     return overland_flow
 end
@@ -274,9 +289,9 @@ function average_flow_vars!(river_flow_model::AbstractRiverFlowModel, dt::Float6
     return nothing
 end
 
-"Update overland flow model `KinWaveOverlandFlowModel` for a single timestep"
+"Update overland flow model `OverlandFlowModel{<:KinematicWave}` for a single timestep"
 function kinwave_land_update!(
-    overland_flow_model::KinWaveOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     domain::DomainLand,
     dt::Float64,
 )
@@ -308,12 +323,12 @@ function kinwave_land_update!(
                     )
                 end
 
-                q[v], crossarea =
+                q[v], flow_area =
                     kinematic_wave(qin[v], q[v], qlat[v], alpha[v], dt, flow_length[v])
 
                 # update h, only if flow width > 0.0
                 if surface_flow_width[v] > 0.0
-                    h[v] = crossarea / surface_flow_width[v]
+                    h[v] = flow_area / surface_flow_width[v]
                 end
                 storage[v] = flow_length[v] * surface_flow_width[v] * h[v]
 
@@ -326,17 +341,16 @@ function kinwave_land_update!(
 end
 
 """
-Update overland flow model `KinWaveOverlandFlowModel` for a single timestep `dt`. Timestepping within
-`dt` is either with a fixed timestep `dt_fixed` or adaptive.
+Update overland flow model `OverlandFlowModel{<:KinematicWave}` for a single timestep `dt`.
+Timestepping within `dt` is either with a fixed timestep `dt_fixed` or adaptive.
 """
 function update_overland_flow_model!(
-    overland_flow_model::KinWaveOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     domain::DomainLand,
     dt::Float64,
 )
     (; inwater) = overland_flow_model.boundary_conditions
-    (; alpha_term, mannings_n, alpha_pow, alpha) = overland_flow_model.parameters
-    (; surface_flow_width, flow_length, slope) = domain.parameters
+    (; flow_length) = domain.parameters
     (;
         q_average,
         q_cumulative,
@@ -348,9 +362,6 @@ function update_overland_flow_model!(
     ) = overland_flow_model.variables
     (; adaptive) = overland_flow_model.timestepping
 
-    @. alpha_term = pow(mannings_n / sqrt(slope), BETA_KINWAVE)
-    # use fixed alpha value based flow width
-    @. alpha = alpha_term * pow(surface_flow_width, alpha_pow)
     @. qlat = inwater / flow_length
 
     q_cumulative .= 0.0
@@ -373,13 +384,63 @@ function update_overland_flow_model!(
     return nothing
 end
 
-"""
-run reservoir model and copy reservoir outflow to inflow (qin) of
-downstream river cell
-"""
+function update_floodplain_model!(
+    river_flow_model::RiverFlowModel{T, F},
+    domain::DomainRiver,
+    dt::Float64,
+) where {T <: KinematicWave, F <: FloodPlainModel{<:Manning}}
+    (; floodplain) = river_flow_model
+    (; flow_length) = domain.parameters
+    (; graph, order) = domain.network
+    (; profile, mannings_n, slope) = floodplain.parameters
+    (; flow_capacity, h, q, q_cumulative, qin, qin_cumulative, storage) =
+        floodplain.variables
+
+    for (i, v) in enumerate(order)
+        i_ds = outneighbors(graph, v)
+        n = length(i_ds)
+        (n > 1) && error("bifurcations not supported")
+        to_pit = iszero(n)
+        if h[v] > 0.0
+            # check if downstream floodplain is active based on minimum flow area
+            i1, i2 = interpolation_indices(h[v], @view profile.depth[:])
+            flow_area = compute_floodplain_flow_area(profile, h[v], v, i1, i2)
+            if !to_pit
+                flow_area_ds =
+                    compute_floodplain_flow_area(profile, h[v], only(i_ds), i1, i2)
+            else
+                flow_area_ds = flow_area
+            end
+            if flow_area > 1.0e-05 && flow_area_ds > 1.0e-05
+                wetted_perimeter = compute_wetted_perimeter(profile, h[v], v, i1)
+                hydraulic_radius = flow_area/wetted_perimeter
+                flow_capacity[v] =
+                    manning_flow(mannings_n[v], hydraulic_radius, slope[v], flow_area)
+            else
+                flow_capacity[v] = 0.0
+            end
+        else
+            flow_capacity[v] = 0.0
+        end
+    end
+    q .= accucapacityflux(storage, domain.network, flow_capacity, dt)
+    @. q_cumulative += q * dt
+
+    flux_in!(qin, q, domain.network)
+
+    @. qin_cumulative += qin * dt
+end
+
+update_floodplain_model!(
+    river_flow_model::RiverFlowModel{T, F},
+    domain::DomainRiver,
+    dt::Float64,
+) where {T <: KinematicWave, F <: Nothing} = nothing
+
+"Run reservoir model and copy reservoir outflow to inflow (qin) of downstream river cell"
 function update_reservoir_model!(
     reservoir_model::ReservoirModel,
-    river_flow_vars::FlowVariables,
+    river_flow_vars::RiverFlowVariables,
     network::NetworkRiver,
     v::Int,
     dt::Float64,
@@ -427,21 +488,27 @@ function update_reservoir_model!(
     return nothing
 end
 
-"Update river flow model `KinWaveRiverFlowModel` for a single timestep"
+"Update river flow model `RiverFlowModel{<:KinematicWave}` for a single timestep"
 function kinwave_river_update!(
-    river_flow_model::KinWaveRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:KinematicWave},
     domain::DomainRiver,
     dt::Float64,
 )
     (; order_of_subdomains, order_subdomain, subdomain_indices, upstream_nodes) =
         domain.network
 
-    (; reservoir, external_inflow, actual_external_abstraction_cumulative, abstraction) =
-        river_flow_model.boundary_conditions
+    (;
+        reservoir,
+        external_inflow,
+        actual_external_abstraction_cumulative,
+        abstraction,
+        floodplain_water_exchange,
+    ) = river_flow_model.boundary_conditions
 
     (; alpha) = river_flow_model.parameters
     (; flow_width, flow_length) = domain.parameters
     (; h, q, q_cumulative, storage, qin, qin_cumulative, qlat) = river_flow_model.variables
+    (; floodplain) = river_flow_model
 
     ns = length(order_of_subdomains)
     qin .= 0.0
@@ -464,7 +531,11 @@ function kinwave_river_update!(
                 # negative external inflow as part of water allocation computations.
                 _inflow -= abstraction[v] / flow_length[v]
 
-                q[v], crossarea = kinematic_wave(
+                if !isnothing(floodplain)
+                    _inflow += floodplain_water_exchange[v] / flow_length[v]
+                end
+
+                q[v], flow_area = kinematic_wave(
                     qin[v],
                     q[v],
                     qlat[v] + _inflow,
@@ -483,8 +554,8 @@ function kinwave_river_update!(
                     )
                 end
                 # update h and storage
-                h[v] = crossarea / flow_width[v]
-                storage[v] = flow_length[v] * flow_width[v] * h[v]
+                h[v] = flow_area / flow_width[v]
+                storage[v] = flow_length[v] * flow_area
 
                 # average variables
                 q_cumulative[v] += q[v] * dt
@@ -495,31 +566,71 @@ function kinwave_river_update!(
 end
 
 """
-Update river flow model `KinWaveRiverFlowModel` for a single timestep `dt`. Timestepping within
-`dt` is either with a fixed timestep `dt_fixed` or adaptive.
+Exchange of channel-floodpain water (assumed instantaneously) based on total storage and
+river channel capacity (`bankfull_storage`).
+"""
+function river_channel_floodplain_exchange!(
+    river_flow_model::RiverFlowModel{T, F},
+    parameters::RiverParameters,
+    dt::Float64,
+) where {T <: KinematicWave, F <: FloodPlainModel{<:Manning}}
+    river_v = river_flow_model.variables
+    river_p = river_flow_model.parameters
+    floodplain_v = river_flow_model.floodplain.variables
+    floodplain_p = river_flow_model.floodplain.parameters
+
+    (; flow_length, flow_width) = parameters
+    (; floodplain_water_exchange) = river_flow_model.boundary_conditions
+
+    for i in eachindex(river_v.qlat)
+        storage_total = river_v.storage[i] + floodplain_v.storage[i]
+        if storage_total > river_p.bankfull_storage[i]
+            flood_storage = storage_total - river_p.bankfull_storage[i]
+            h = compute_flood_depth(floodplain_p.profile, flood_storage, flow_length[i], i)
+            river_storage = (river_p.bankfull_depth[i] + h) * flow_width[i] * flow_length[i]
+            delta_river_storage = river_storage - river_v.storage[i]
+            floodplain_v.storage[i] = max(storage_total - river_storage, 0.0)
+            floodplain_v.h[i] = floodplain_v.storage[i] > 0.0 ? h : 0.0
+        else
+            delta_river_storage = max(storage_total - river_v.storage[i], 0.0)
+            floodplain_v.h[i] = 0.0
+            floodplain_v.storage[i] = 0.0
+        end
+        floodplain_water_exchange[i] = delta_river_storage/dt
+    end
+end
+
+river_channel_floodplain_exchange!(
+    river_flow_model::RiverFlowModel{T, F},
+    parameters::RiverParameters,
+    dt::Float64,
+) where {T <: KinematicWave, F <: Nothing} = nothing
+
+"""
+Update river flow model `RiverFlowModel{<:KinematicWave}` for a single timestep `dt`.
+Timestepping within `dt` is either with a fixed timestep `dt_fixed` or adaptive.
 """
 function update_river_flow_model!(
-    river_flow_model::KinWaveRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:KinematicWave},
     domain::Domain,
     clock::Clock,
     dt::Float64,
 )
+    (; floodplain) = river_flow_model
     (; reservoir, inwater) = river_flow_model.boundary_conditions
-    (; alpha_term, mannings_n, alpha_pow, alpha, bankfull_depth) =
-        river_flow_model.parameters
-    (; slope, flow_width, flow_length) = domain.river.parameters
+    (; flow_length) = domain.river.parameters
     (; qlat, qin_average, qin_cumulative) = river_flow_model.variables
     (; adaptive) = river_flow_model.timestepping
 
-    @. alpha_term = pow(mannings_n / sqrt(slope), BETA_KINWAVE)
-    # use fixed alpha value based on 0.5 * bankfull_depth
-    @. alpha = alpha_term * pow(flow_width + bankfull_depth, alpha_pow)
     @. qlat = inwater / flow_length
-
     set_flow_vars!(river_flow_model)
     qin_cumulative .= 0.0
     set_reservoir_vars!(reservoir)
     update_index_hq!(reservoir, clock)
+    if !isnothing(floodplain)
+        floodplain.variables.q_cumulative .= 0.0
+        floodplain.variables.qin_cumulative .= 0.0
+    end
 
     t = 0.0
     while t < dt
@@ -527,20 +638,29 @@ function update_river_flow_model!(
             adaptive ? stable_timestep(river_flow_model, flow_length, 0.05) :
             river_flow_model.timestepping.dt_fixed
         dt_s = check_timestepsize(dt_s, t, dt)
+        river_channel_floodplain_exchange!(river_flow_model, domain.river.parameters, dt_s)
         kinwave_river_update!(river_flow_model, domain.river, dt_s)
+        update_floodplain_model!(river_flow_model, domain.river, dt_s)
         t += dt_s
     end
 
     average_reservoir_vars!(reservoir, dt)
     average_flow_vars!(river_flow_model, dt)
     @. qin_average = qin_cumulative / dt
+    if !isnothing(floodplain)
+        v_floodplain = river_flow_model.floodplain.variables
+        v_river = river_flow_model.variables
+
+        @. v_floodplain.q_average = v_floodplain.q_cumulative / dt
+        @. v_river.q_channel_average = v_river.q_average
+        @. v_river.q_average = v_river.q_channel_average + v_floodplain.q_average
+        @. v_floodplain.qin_average = v_floodplain.qin_cumulative / dt
+        @. v_river.qin_average = v_river.qin_average + v_floodplain.qin_average
+    end
 
     return nothing
 end
 
-# constant in Manning's equation [-]
-const BETA_KINWAVE = 0.6
-
 """
 Compute a stable timestep size for the kinematice wave method for a river or overland flow
 model using a nonlinear scheme (Chow et al., 1988).
@@ -555,13 +675,14 @@ function stable_timestep(
     flow_model::S,
     flow_length::Vector{Float64},
     p::Float64,
-) where {S <: Union{KinWaveOverlandFlowModel, KinWaveRiverFlowModel}}
+) where {S <: Union{RiverFlowModel{<:KinematicWave}, OverlandFlowModel{<:KinematicWave}}}
     (; q) = flow_model.variables
     (; alpha) = flow_model.parameters
     (; stable_timesteps) = flow_model.timestepping
 
     n = length(q)
     stable_timesteps .= Inf
+    dt_min_default = 600.0
     k = 0
     for i in 1:n
         if q[i] > KIN_WAVE_MIN_FLOW
@@ -576,7 +697,7 @@ function stable_timestep(
     elseif k > 0
         quantile!(@view(stable_timesteps[1:k]), p)
     else
-        600.0
+        dt_min_default
     end
 
     return dt_min
@@ -614,10 +735,10 @@ end
 
 """
 Update boundary condition lateral inflow `inwater` of a kinematic wave overland flow model
-`KinWaveOverlandFlowModel` for a single timestep.
+`OverlandFlowModel{<:KinematicWave}` for a single timestep.
 """
 function update_lateral_inflow!(
-    overland_flow_model::KinWaveOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     external_models::NamedTuple,
     domain::Domain,
     config::Config,
@@ -669,12 +790,13 @@ end
 # For the river kinematic wave, the variable `to_river` can be excluded, because this part
 # is added to the river kinematic wave.
 get_inflow_reservoir(
-    river_flow_model::KinWaveRiverFlowModel,
-    overland_flow_model::KinWaveOverlandFlowModel,
+    river_flow_model::RiverFlowModel{<:KinematicWave},
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     inds::Vector{Int},
 ) = overland_flow_model.variables.q_average[inds]
+
 get_inflow_reservoir(
-    river_flow_model::KinWaveRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:KinematicWave},
     subsurface_flow_model::LateralSSFModel,
     inds::Vector{Int},
 ) = subsurface_flow_model.variables.q_average[inds]
diff --git a/Wflow/src/routing/surface/surface_process.jl b/Wflow/src/routing/surface/surface_process.jl
index 9cd70779e..f84bbe029 100644
--- a/Wflow/src/routing/surface/surface_process.jl
+++ b/Wflow/src/routing/surface/surface_process.jl
@@ -161,3 +161,9 @@ function local_inertial_flow(
 
     return q
 end
+
+"Compute flow using Manning's equation"
+function manning_flow(mannings_n, hydraulic_radius, slope, area)
+    q = cbrt(hydraulic_radius^2) * sqrt(slope) * area / mannings_n
+    return q
+end
diff --git a/Wflow/src/routing/surface/surface_routing.jl b/Wflow/src/routing/surface/surface_routing.jl
index 9d3f52f39..7fb160eca 100644
--- a/Wflow/src/routing/surface/surface_routing.jl
+++ b/Wflow/src/routing/surface/surface_routing.jl
@@ -47,15 +47,23 @@ end
 
 """
     surface_routing!(
-        model::Model{R}
-    ) where {R <: Routing{<:LocalInertialOverlandFlowModel, <:LocalInertialRiverFlowModel}}
+    model::Model{R},
+) where {
+    R <: Routing{
+        <:OverlandFlowModel{<:LocalInertial},
+        <:RiverFlowModel{<:LocalInertial},
+    },
+}
 
 Run surface routing (land and river) for a model type that contains the routing components
-`LocalInertialOverlandFlowModel` and `LocalInertialRiverFlowModel` for a single timestep.
+`OverlandFlowModel{<:LocalInertial}` and `RiverFlowModel{<:LocalInertial}` for a single
+timestep.
 """
 function surface_routing!(
     model::Model{R},
-) where {R <: Routing{<:LocalInertialOverlandFlowModel, <:LocalInertialRiverFlowModel}}
+) where {
+    R <: Routing{<:OverlandFlowModel{<:LocalInertial}, <:RiverFlowModel{<:LocalInertial}},
+}
     (; routing, land, domain, clock, config) = model
     (; soil, runoff) = land
     (; overland_flow, river_flow, subsurface_flow) = routing
diff --git a/Wflow/src/routing/surface/surface_local_inertial.jl b/Wflow/src/routing/surface/surface_staggered_scheme.jl
similarity index 55%
rename from Wflow/src/routing/surface/surface_local_inertial.jl
rename to Wflow/src/routing/surface/surface_staggered_scheme.jl
index d9e8b7d4d..a9fbd325d 100644
--- a/Wflow/src/routing/surface/surface_local_inertial.jl
+++ b/Wflow/src/routing/surface/surface_staggered_scheme.jl
@@ -1,67 +1,53 @@
-abstract type AbstractFloodPlainModel end
-
-"Struct for storing local inertial river flow model parameters"
-@with_kw struct LocalInertialRiverFlowParameters
+"Struct for storing river flow model parameters on a staggered grid"
+@with_kw struct RiverFlowStaggeredParameters <: AbstractRiverFlowParameters
     # number of cells [-]
     n::Int
     # number of edges [-]
-    ne::Int
+    n_edges::Int
     # active nodes [-]
     active_n::Vector{Int}
     # active edges [-]
     active_e::Vector{Int}
     # if true a check is performed if froude number > 1.0 (algorithm is modified) [-]
-    froude_limit::Bool
+    froude_limit::Bool = false
     # depth threshold for calculating flow [m]
-    h_thresh::Float64
+    h_thresh::Float64 = 0.0
     # river bed elevation [m]
-    zb::Vector{Float64}
-    # maximum channel bed elevation [m]
-    zb_max::Vector{Float64}
+    zb::Vector{Float64} = Float64[]
     # bankfull storage [m³]
-    bankfull_storage::Vector{Float64}
+    bankfull_storage::Vector{Float64} = Float64[]
     # bankfull depth [m]
-    bankfull_depth::Vector{Float64}
-    # Manning's roughness squared at edge [(s m-1/3)²]
-    mannings_n_sq::Vector{Float64}
-    # Manning's roughness [s m-1/3] at node
-    mannings_n::Vector{Float64}
+    bankfull_depth::Vector{Float64} = Float64[]
+    # channel bed elevation at edge [m]
+    zb_at_edge::Vector{Float64} = Float64[]
+    # Manning's roughness [s m-1/3]
+    mannings_n::Vector{Float64} = Float64[]
+    # Manning's roughness at edge [s m-1/3]
+    mannings_n_at_edge::Vector{Float64} = Float64[]
+    # Manning's roughness squared at edge [(s m-1/3)2]
+    mannings_n_sq_at_edge::Vector{Float64} = Float64[]
+    # slope at edge [-]
+    slope_at_edge::Vector{Float64} = Float64[]
     # flow (river) length at edge [m]
-    flow_length_at_edge::Vector{Float64}
+    flow_length_at_edge::Vector{Float64} = Float64[]
     # flow (river) width at edge [m]
-    flow_width_at_edge::Vector{Float64}
+    flow_width_at_edge::Vector{Float64} = Float64[]
 end
 
-"Initialize local inertial river flow model parameters"
-function LocalInertialRiverFlowParameters(
+"Get river parameters for a river flow model on a staggered grid"
+function get_river_parameters(
     dataset::NCDataset,
     config::Config,
     domain::DomainRiver,
+    index_pit::Vector{Int},
 )
-    alpha = config.model.river_local_inertial_flow__alpha_coefficient # stability coefficient for model time step (0.2-0.7)
-    waterdepth_threshold = config.model.river_water_flow_threshold__depth # depth threshold for flow at edge
-    froude_limit = config.model.river_water_flow__froude_limit_flag # limit flow to subcritical according to Froude number
-    floodplain_1d = config.model.floodplain_1d__flag
-
-    @info "Local inertial approach is used for river flow." alpha waterdepth_threshold froude_limit floodplain_1d
-
-    (; pit_indices, indices, graph, local_drain_direction, nodes_at_edge) = domain.network
-    (; flow_width, flow_length, reservoir_outlet) = domain.parameters
+    (; pit_indices, indices) = domain.network
+    (; flow_length, flow_width) = domain.parameters
 
-    riverlength_bc = ncread(
-        dataset,
-        config,
-        "model_boundary_condition_river__length",
-        Routing;
-        sel = pit_indices,
-    )
     bankfull_elevation =
         ncread(dataset, config, "river_bank_water__elevation", Routing; sel = indices)
     bankfull_depth =
         ncread(dataset, config, "river_bank_water__depth", Routing; sel = indices)
-    zb = bankfull_elevation - bankfull_depth # river bed elevation
-
-    bankfull_storage = bankfull_depth .* flow_width .* flow_length
     mannings_n = ncread(
         dataset,
         config,
@@ -69,59 +55,107 @@ function LocalInertialRiverFlowParameters(
         Routing;
         sel = indices,
     )
-
-    n = length(indices)
-    index_pit = findall(==(5), local_drain_direction)
-    # set ghost points for boundary condition (downstream river outlet): river width, bed
-    # elevation, manning n is copied from the upstream cell.
+    riverlength_bc = ncread(
+        dataset,
+        config,
+        "model_boundary_condition_river__length",
+        Routing;
+        sel = pit_indices,
+    )
     append!(flow_length, riverlength_bc)
-    append!(zb, zb[index_pit])
+
+    # set ghost points for boundary condition (downstream river outlet): flow width,
+    # mannings n, bankfull depth and elevation is copied from the upstream cell.
     append!(flow_width, flow_width[index_pit])
     append!(mannings_n, mannings_n[index_pit])
     append!(bankfull_depth, bankfull_depth[index_pit])
+    append!(bankfull_elevation, bankfull_elevation[index_pit])
+
+    return mannings_n, bankfull_depth, bankfull_elevation, flow_length, flow_width
+end
+
+"""
+Generate info log message for river flow routing on a staggered grid, solving the local
+inertial or Manning's equation.
+"""
+function log_message_staggered_flow(config::Config)
+    (; river_routing) = config.model
+    waterdepth_threshold = config.model.river_water_flow_threshold__depth # depth threshold for flow at edge
+    floodplain_1d = config.model.floodplain_1d__flag
+
+    if river_routing == RoutingType.local_inertial
+        alpha = config.model.river_local_inertial_flow__alpha_coefficient # stability coefficient for model time step (0.2-0.7)
+        froude_limit = config.model.river_water_flow__froude_limit_flag # limit flow to subcritical according to Froude number
+        @info "Local inertial approach is used for river flow." alpha waterdepth_threshold froude_limit floodplain_1d
+    elseif river_routing == RoutingType.manning_staggered
+        alpha = config.model.river_staggered_manning_flow__alpha_coefficient
+        froude_limit = false
+        @info "Manning's equation on a staggered grid is used for river flow." alpha floodplain_1d
+    end
+    return waterdepth_threshold, froude_limit
+end
+
+"Initialize river flow model parameters on a staggered grid."
+function RiverFlowStaggeredParameters(
+    dataset::NCDataset,
+    config::Config,
+    domain::DomainRiver,
+)
+    (; river_routing) = config.model
+    (; graph, indices, local_drain_direction, nodes_at_edge) = domain.network
+    (; reservoir_outlet) = domain.parameters
 
-    # determine z, width, length and manning's n at edges
+    waterdepth_threshold, froude_limit = log_message_staggered_flow(config)
+
+    n = length(indices)
     n_edges = ne(graph)
-    zb_max = fill(Float64(0), n_edges)
-    width_at_edge = fill(Float64(0), n_edges)
-    length_at_edge = fill(Float64(0), n_edges)
-    mannings_n_sq = fill(Float64(0), n_edges)
-    for i in 1:n_edges
-        src_node = nodes_at_edge.src[i]
-        dst_node = nodes_at_edge.dst[i]
-        zb_max[i] = max(zb[src_node], zb[dst_node])
-        width_at_edge[i] = min(flow_width[src_node], flow_width[dst_node])
-        length_at_edge[i] = (flow_length[dst_node] + flow_length[src_node]) / 2
-        mannings_n_i =
-            (
-                mannings_n[dst_node] * flow_length[dst_node] +
-                mannings_n[src_node] * flow_length[src_node]
-            ) / (flow_length[dst_node] + flow_length[src_node])
-        mannings_n_sq[i] = mannings_n_i * mannings_n_i
+    active_index = findall(x -> x == 0, reservoir_outlet)
+    index_pit = findall(x -> x == LDD_PIT, local_drain_direction)
+
+    mannings_n, bankfull_depth, bankfull_elevation, flow_length, flow_width =
+        get_river_parameters(dataset, config, domain, index_pit)
+    zb = bankfull_elevation - bankfull_depth # river bed elevation
+    bankfull_storage = bankfull_depth .* flow_width .* flow_length
+
+    # determine parameters at edges
+    flow_width_at_edge = compute_value_at_edge(flow_width, nodes_at_edge, n_edges, minimum)
+    flow_length_at_edge = compute_value_at_edge(flow_length, nodes_at_edge, n_edges, mean)
+    mannings_n_at_edge =
+        compute_mannings_n_at_edge(mannings_n, flow_length, nodes_at_edge, n_edges)
+    if river_routing == RoutingType.local_inertial
+        zb_at_edge = compute_value_at_edge(zb, nodes_at_edge, n_edges, maximum)
+        mannings_n_sq_at_edge = mannings_n_at_edge .* mannings_n_at_edge
+        slope_at_edge = []
+    elseif river_routing == RoutingType.manning_staggered
+        zb_at_edge = compute_value_at_edge(zb, nodes_at_edge, n_edges, first)
+        mannings_n_sq_at_edge = []
+        slope_at_edge =
+            compute_slope_at_edge(zb, flow_length_at_edge, nodes_at_edge, n_edges)
     end
-    active_index = findall(==(0), reservoir_outlet)
 
-    parameters = LocalInertialRiverFlowParameters(;
+    parameters = RiverFlowStaggeredParameters(;
         n,
-        ne = n_edges,
+        n_edges,
         active_n = active_index,
         active_e = active_index,
         froude_limit,
         h_thresh = waterdepth_threshold,
         zb,
-        zb_max,
         bankfull_storage,
         bankfull_depth,
+        zb_at_edge,
         mannings_n,
-        mannings_n_sq,
-        flow_length_at_edge = length_at_edge,
-        flow_width_at_edge = width_at_edge,
+        mannings_n_at_edge,
+        mannings_n_sq_at_edge,
+        slope_at_edge,
+        flow_length_at_edge,
+        flow_width_at_edge,
     )
     return parameters
 end
 
-"Struct for storing local inertial river flow model variables"
-@with_kw struct LocalInertialRiverFlowVariables
+"Struct for storing river flow model variables on a staggered grid."
+@with_kw struct RiverFlowStaggeredVariables <: AbstractRiverFlowVariables
     n_cells::Int
     n_edges::Int
     # river discharge at edge (subgrid channel) [m³ s⁻¹]
@@ -138,31 +172,28 @@ end
     q_channel_average::Vector{Float64} = q_average
     # water depth [m]
     h::Vector{Float64}
-    # maximum water elevation at edge [m]
-    zs_max::Vector{Float64} = zeros(n_edges)
+    # water elevation at edge [m]
+    zs_at_edge::Vector{Float64} = zeros(n_edges)
     # water elevation of source node of edge [m]
     zs_src::Vector{Float64} = zeros(n_edges)
     # water elevation of downstream node of edge [m]
     zs_dst::Vector{Float64} = zeros(n_edges)
     # water depth at edge [m]
     water_depth_at_edge::Vector{Float64} = zeros(n_edges)
-    # flow area at edge [m²]
-    flow_area::Vector{Float64} = zeros(n_edges)
-    # wetted perimeter at edge [m]
-    hydraulic_radius::Vector{Float64} = zeros(n_edges)
     # river storage [m³]
     storage::Vector{Float64} = zeros(n_cells)
     # error storage [m³]
     error::Vector{Float64} = zeros(n_cells)
 end
 
-"Initialize shallow water river flow model variables"
-function LocalInertialRiverFlowVariables(
+"Initialize river flow model variables on a staggered grid."
+function RiverFlowStaggeredVariables(
     dataset::NCDataset,
     config::Config,
     network::NetworkRiver,
 )
     (; pit_indices, indices, graph) = network
+    (; river_routing) = config.model
 
     riverdepth_bc = ncread(
         dataset,
@@ -177,10 +208,11 @@ function LocalInertialRiverFlowVariables(
     # set ghost points for boundary condition (downstream river outlet): river depth `h`
     h = zeros(n_cells)
     append!(h, riverdepth_bc)
+
     variables = if config.model.floodplain_1d__flag
         # When floodplain is enabled, q_channel fields must be separate from q fields
         # because q_average will hold channel + floodplain combined discharge.
-        LocalInertialRiverFlowVariables(;
+        RiverFlowStaggeredVariables(;
             n_cells,
             n_edges,
             h,
@@ -189,38 +221,22 @@ function LocalInertialRiverFlowVariables(
         )
     else
         # When floodplain is disabled, q_channel aliases q (they are identical)
-        LocalInertialRiverFlowVariables(; n_cells, n_edges, h)
+        RiverFlowStaggeredVariables(; n_cells, n_edges, h)
     end
     return variables
 end
 
-"Shallow water river flow model using the local inertial method"
-@with_kw struct LocalInertialRiverFlowModel{
-    R <: RiverFlowBC,
-    F <: Union{AbstractFloodPlainModel, Nothing},
-    A <: AbstractAllocationModel,
-} <: AbstractRiverFlowModel
-    timestepping::TimeStepping
-    boundary_conditions::R
-    parameters::LocalInertialRiverFlowParameters
-    variables::LocalInertialRiverFlowVariables
-    # Floodplain (1D) schematization
-    floodplain::F
-    # Water allocation
-    allocation::A
-end
-
-"Initialize shallow water river flow model `LocalInertialRiverFlowModel`"
-function LocalInertialRiverFlowModel(
+"Initialize river flow model on a staggered grid."
+function init_staggered_river_flow(
     dataset::NCDataset,
     config::Config,
     domain::DomainRiver,
     reservoir_model::Union{ReservoirModel, Nothing},
 )
-    # The local inertial approach makes use of a staggered grid (Bates et al. (2010)),
-    # with nodes and edges. This information is extracted from the directed graph of the
-    # river. Discharge q is calculated at edges between nodes and mapped to the source
-    # nodes for gridded output (index of edge is equal to source node index, e.g.:
+    # This river flow model makes use of a staggered grid (Bates et al. (2010)), with nodes
+    # and edges. This information is extracted from the directed graph of the river.
+    # Discharge q is calculated at edges between nodes and mapped to the source nodes for
+    # gridded output (index of edge is equal to source node index, e.g.:
     # Edge 1 => 5
     # Edge 2 => 1
     # Edge 3 => 2
@@ -229,11 +245,18 @@ function LocalInertialRiverFlowModel(
 
     # The following boundary conditions can be set at ghost nodes, downstream of river
     # outlets (pits): river length and river depth
-    alpha_coefficient = config.model.river_local_inertial_flow__alpha_coefficient # stability coefficient for model time step (0.2-0.7)
+
+    # stability coefficient for model time step (0.2-0.7)
+    (; river_routing) = config.model
+    alpha_coefficient = if river_routing == RoutingType.local_inertial
+        config.model.river_local_inertial_flow__alpha_coefficient
+    elseif river_routing == RoutingType.manning_staggered
+        config.model.river_staggered_manning_flow__alpha_coefficient
+    end
     timestepping = TimeStepping(; alpha_coefficient)
 
-    parameters = LocalInertialRiverFlowParameters(dataset, config, domain)
-    variables = LocalInertialRiverFlowVariables(dataset, config, domain.network)
+    parameters = RiverFlowStaggeredParameters(dataset, config, domain)
+    variables = RiverFlowStaggeredVariables(dataset, config, domain.network)
     boundary_conditions = RiverFlowBC(dataset, config, domain.network, reservoir_model)
 
     if config.model.floodplain_1d__flag
@@ -243,8 +266,15 @@ function LocalInertialRiverFlowModel(
         floodplain = nothing
     end
 
+    (; river_routing) = config.model
+    routing_method = if river_routing == RoutingType.local_inertial
+        LocalInertial()
+    elseif river_routing == RoutingType.manning_staggered
+        ManningStaggered()
+    end
+
     n = length(domain.network.indices)
-    river_flow = LocalInertialRiverFlowModel(;
+    river_flow = RiverFlowModel(;
         timestepping,
         boundary_conditions,
         parameters,
@@ -252,13 +282,16 @@ function LocalInertialRiverFlowModel(
         floodplain,
         allocation = do_water_demand(config) ? AllocationRiverModel(n) :
                      NoAllocationRiverModel(n),
+        routing_method,
     )
     return river_flow
 end
 
-"Return the upstream inflow for flow_area reservoir in `LocalInertialRiverFlowModel`"
+"""
+Return the upstream inflow for a reservoir in a river flow model on a staggered grid.
+"""
 function get_inflow_reservoir(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     src_edge::Vector{Int},
 )
     q_in = sum_at(river_flow_model.variables.q, src_edge)
@@ -271,15 +304,15 @@ end
 # For local inertial river routing, `to_river` is included, as reservoir cells are excluded
 # (boundary condition).
 get_inflow_reservoir(
-    river_flow_model::LocalInertialRiverFlowModel,
-    overland_flow_model::KinWaveOverlandFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
+    overland_flow_model::OverlandFlowModel{<:KinematicWave},
     inds::Vector{Int},
 ) =
     overland_flow_model.variables.q_average[inds] .+
     overland_flow_model.variables.to_river_average[inds]
 
 get_inflow_reservoir(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     subsurface_flow_model::LateralSSFModel,
     inds::Vector{Int},
 ) = (
@@ -291,39 +324,7 @@ get_inflow_reservoir(
 Update river channel flow for the local inertial river flow model.
 """
 function update_river_channel_flow!(
-    river_flow_model::LocalInertialRiverFlowModel,
-    domain::DomainRiver,
-    dt::Float64,
-    update_h::Bool,
-)
-    # Update river channel flow
-    update_river_channel_flow!(river_flow_model, domain.river, dt)
-
-    # Update floodplain flow if present
-    update_floodplain_flow!(river_flow_model, domain.river, dt)
-
-    # Handle reservoir boundary conditions
-    update_bc_reservoir_model!(river_flow_model, domain, dt)
-
-    # Update water depth and storage if requested
-    if update_h
-        update_water_depth_and_storage!(river_flow_model, domain.river, dt)
-        update_water_depth_and_storage!(
-            river_flow_model.floodplain,
-            river_flow_model,
-            domain.river,
-            dt,
-        )
-    end
-
-    return nothing
-end
-
-"""
-Update river channel flow for the local inertial river flow model.
-"""
-function update_river_channel_flow!(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     domain::DomainRiver,
     dt::Float64,
 )
@@ -344,14 +345,16 @@ function update_river_channel_flow!(
         river_v.zs_src[i] = river_p.zb[i_src] + river_v.h[i_src]
         river_v.zs_dst[i] = river_p.zb[i_dst] + river_v.h[i_dst]
 
-        river_v.zs_max[i] = max(river_v.zs_src[i], river_v.zs_dst[i])
-        river_v.water_depth_at_edge[i] = (river_v.zs_max[i] - river_p.zb_max[i])
+        river_v.zs_at_edge[i] = max(river_v.zs_src[i], river_v.zs_dst[i])
+        river_v.water_depth_at_edge[i] = (river_v.zs_at_edge[i] - river_p.zb_at_edge[i])
 
-        river_v.flow_area[i] =
-            river_p.flow_width_at_edge[i] * river_v.water_depth_at_edge[i] # flow area (rectangular channel)
-        river_v.hydraulic_radius[i] =
-            river_v.flow_area[i] /
-            (river_p.flow_width_at_edge[i] + 2 * river_v.water_depth_at_edge[i]) # hydraulic radius (rectangular channel)
+        # rectangular channel
+        flow_area_at_edge = river_p.flow_width_at_edge[i] * river_v.water_depth_at_edge[i]
+        hydraulic_radius_at_edge =
+            flow_area_at_edge / wetted_perimeter_channel(
+                river_v.water_depth_at_edge[i],
+                river_p.flow_width_at_edge[i],
+            )
 
         river_v.q[i] = ifelse(
             river_v.water_depth_at_edge[i] > river_p.h_thresh,
@@ -360,10 +363,10 @@ function update_river_channel_flow!(
                 river_v.zs_src[i],
                 river_v.zs_dst[i],
                 river_v.water_depth_at_edge[i],
-                river_v.flow_area[i],
-                river_v.hydraulic_radius[i],
+                flow_area_at_edge,
+                hydraulic_radius_at_edge,
                 river_p.flow_length_at_edge[i],
-                river_p.mannings_n_sq[i],
+                river_p.mannings_n_sq_at_edge[i],
                 river_p.froude_limit,
                 dt,
             ),
@@ -379,92 +382,135 @@ function update_river_channel_flow!(
     return nothing
 end
 
+"""
+Update river channel flow using Manning's equation on a staggered grid.
+"""
+function update_river_channel_flow!(
+    river_flow_model::RiverFlowModel{<:ManningStaggered},
+    domain::DomainRiver,
+    dt::Float64,
+)
+    (; nodes_at_edge) = domain.network
+    river_v = river_flow_model.variables
+    river_p = river_flow_model.parameters
+    (; inwater, abstraction) = river_flow_model.boundary_conditions
+
+    @batch per = thread minbatch = 1000 for j in eachindex(river_p.active_e)
+        i = river_p.active_e[j]
+        i_src = nodes_at_edge.src[i]
+        i_dst = nodes_at_edge.dst[i]
+        river_v.zs_src[i] = river_p.zb[i_src] + river_v.h[i_src]
+        river_v.zs_dst[i] = river_p.zb[i_dst] + river_v.h[i_dst]
+
+        # use "upwind" option for water level and depth estimation at edge to avoid large q
+        # fluctuations (compared to "max" option)
+        river_v.zs_at_edge[i] = river_v.zs_src[i]
+        river_v.water_depth_at_edge[i] = river_v.h[i_src]
+
+        # rectangular channel
+        flow_area_at_edge = river_p.flow_width_at_edge[i] * river_v.water_depth_at_edge[i]
+        hydraulic_radius_at_edge =
+            flow_area_at_edge / wetted_perimeter_channel(
+                river_v.water_depth_at_edge[i],
+                river_p.flow_width_at_edge[i],
+            )
+
+        river_v.q[i] = ifelse(
+            river_v.water_depth_at_edge[i] > river_p.h_thresh,
+            manning_flow(
+                river_p.mannings_n_at_edge[i],
+                hydraulic_radius_at_edge,
+                river_p.slope_at_edge[i],
+                flow_area_at_edge,
+            ),
+            0.0,
+        )
+
+        # limit q in case water is not available
+        river_v.q[i] = min(river_v.q[i], river_v.storage[i_src]/dt)
+        # average river discharge (here accumulated for model timestep Δt)
+        river_v.q_cumulative[i] += river_v.q[i] * dt
+    end
+    return nothing
+end
+
 """
 Update floodplain flow for the local inertial river flow model.
 """
 function update_floodplain_flow!(
-    river_flow_model::LocalInertialRiverFlowModel{R, F},
+    river_flow_model::RiverFlowModel{T, F},
     domain::DomainRiver,
     dt::Float64,
-) where {R, F <: AbstractFloodPlainModel}
+) where {T <: LocalInertial, F <: FloodPlainModel{<:LocalInertial}}
     (; nodes_at_edge) = domain.network
     (; flow_width) = domain.parameters
 
     river_v = river_flow_model.variables
     river_p = river_flow_model.parameters
     floodplain_p = river_flow_model.floodplain.parameters
+    (; profile) = floodplain_p
     floodplain_v = river_flow_model.floodplain.variables
 
-    @batch per = thread minbatch = 1000 for i in 1:length(floodplain_v.water_depth_at_edge)
-        floodplain_v.water_depth_at_edge[i] =
-            max(river_v.zs_max[i] - floodplain_p.zb_max[i], 0.0)
-    end
-
-    n = 0
-    @inbounds for i in river_p.active_e
-        @inbounds if river_v.water_depth_at_edge[i] > river_p.h_thresh
-            n += 1
-            floodplain_v.hf_index[n] = i
-        else
-            floodplain_v.q[i] = 0.0
-        end
-    end
-
-    get_area(i, i1, i2, idx) = flow_area(
-        floodplain_p.profile.width[i2, idx],
-        floodplain_p.profile.flow_area[i1, idx],
-        floodplain_p.profile.depth[i1],
-        floodplain_v.water_depth_at_edge[i],
-    )
-
-    get_wetted_perimeter(i, i1, idx) = wetted_perimeter(
-        floodplain_p.profile.wetted_perimeter[i1, idx],
-        floodplain_p.profile.depth[i1],
-        floodplain_v.water_depth_at_edge[i],
-    )
-
-    @batch per = thread minbatch = 1000 for j in 1:n
-        i = floodplain_v.hf_index[j]
+    @batch per = thread minbatch = 1000 for j in eachindex(river_p.active_e)
+        i = river_p.active_e[j]
         i_src = nodes_at_edge.src[i]
         i_dst = nodes_at_edge.dst[i]
 
-        i0 = 0
-        for k in eachindex(floodplain_p.profile.depth)
-            i0 += 1 * (floodplain_p.profile.depth[k] <= floodplain_v.water_depth_at_edge[i])
-        end
-        i1 = max(i0, 1)
-        i2 = ifelse(i1 == length(floodplain_p.profile.depth), i1, i1 + 1)
-
-        a_src = get_area(i, i1, i2, i_src)
-        a_src = max(a_src - (floodplain_v.water_depth_at_edge[i] * flow_width[i_src]), 0.0)
-
-        a_dst = get_area(i, i1, i2, i_dst)
-        a_dst = max(a_dst - (floodplain_v.water_depth_at_edge[i] * flow_width[i_dst]), 0.0)
+        floodplain_v.water_depth_at_edge[i] =
+            max(river_v.zs_at_edge[i] - floodplain_p.zb_at_edge[i], 0.0)
 
-        floodplain_v.flow_area[i] = min(a_src, a_dst)
+        i1, i2 = interpolation_indices(
+            floodplain_v.water_depth_at_edge[i],
+            @view profile.depth[:]
+        )
+        a_src = compute_floodplain_flow_area(
+            profile,
+            floodplain_v.water_depth_at_edge[i],
+            i_src,
+            i1,
+            i2,
+        )
+        a_dst = compute_floodplain_flow_area(
+            profile,
+            floodplain_v.water_depth_at_edge[i],
+            i_dst,
+            i1,
+            i2,
+        )
+        flow_area_at_edge = min(a_src, a_dst)
 
-        floodplain_v.hydraulic_radius[i] = if a_src < a_dst
-            a_src / get_wetted_perimeter(i, i1, i_src)
-        else
-            a_dst / get_wetted_perimeter(i, i1, i_dst)
-        end
+        hydraulic_radius_at_edge = ifelse(
+            a_src < a_dst,
+            a_src / compute_wetted_perimeter(
+                profile,
+                floodplain_v.water_depth_at_edge[i],
+                i_src,
+                i1,
+            ),
+            a_dst / compute_wetted_perimeter(
+                profile,
+                floodplain_v.water_depth_at_edge[i],
+                i_dst,
+                i1,
+            ),
+        )
 
-        floodplain_v.q[i] = if floodplain_v.flow_area[i] > 1.0e-05
+        floodplain_v.q[i] = ifelse(
+            flow_area_at_edge > 1.0e-05,
             local_inertial_flow(
                 floodplain_v.q_previous[i],
                 river_v.zs_src[i],
                 river_v.zs_dst[i],
                 floodplain_v.water_depth_at_edge[i],
-                floodplain_v.flow_area[i],
-                floodplain_v.hydraulic_radius[i],
+                flow_area_at_edge,
+                hydraulic_radius_at_edge,
                 river_p.flow_length_at_edge[i],
-                floodplain_p.mannings_n_sq[i],
+                floodplain_p.mannings_n_sq_at_edge[i],
                 river_p.froude_limit,
                 dt,
-            )
-        else
-            0.0
-        end
+            ),
+            0.0,
+        )
 
         # limit floodplain q in case water is not available
         if floodplain_v.h[i_src] <= 0.0
@@ -485,55 +531,102 @@ function update_floodplain_flow!(
     end
     return nothing
 end
-update_floodplain_flow!(
-    river_flow_model::LocalInertialRiverFlowModel{R, Nothing},
-    domain::DomainRiver,
-    dt::Float64,
-) where {R} = nothing
 
-function update_bc_reservoir_model!(
-    river_flow_model::LocalInertialRiverFlowModel,
-    domain::Domain,
+"""
+Update floodplain flow for a river flow model using Manning's equation on a staggered grid.
+"""
+function update_floodplain_flow!(
+    river_flow_model::RiverFlowModel{T, F},
+    domain::DomainRiver,
     dt::Float64,
-)
-    (; edges_at_node) = domain.river.network
-    (; reservoir) = river_flow_model.boundary_conditions
-    inds_reservoir = domain.reservoir.network.river_indices
-    isnothing(reservoir) && return nothing
+) where {T <: ManningStaggered, F <: FloodPlainModel{<:ManningStaggered}}
+    (; nodes_at_edge) = domain.network
+    (; flow_width) = domain.parameters
 
     river_v = river_flow_model.variables
-    res_bc = reservoir.boundary_conditions
+    river_p = river_flow_model.parameters
+    floodplain_p = river_flow_model.floodplain.parameters
+    (; profile) = floodplain_p
+    floodplain_v = river_flow_model.floodplain.variables
 
-    for v in eachindex(inds_reservoir)
-        i = inds_reservoir[v]
+    @batch per = thread minbatch = 1000 for j in eachindex(river_p.active_e)
+        i = river_p.active_e[j]
+        i_src = nodes_at_edge.src[i]
+        i_dst = nodes_at_edge.dst[i]
 
-        q_in = get_inflow_reservoir(river_flow_model, edges_at_node.src[i])
-        # If external_inflow < 0, abstraction is limited
-        if res_bc.external_inflow[v] < 0.0
-            abstraction = min(
-                -res_bc.external_inflow[v],
-                (reservoir.variables.storage[v] / dt) * 0.98,
-            )
-            res_bc.actual_external_abstraction_cumulative[v] += abstraction * dt
-            inflow = -abstraction
-        else
-            inflow = res_bc.external_inflow[v]
-        end
-        net_inflow = q_in + res_bc.inflow_overland[v] + res_bc.inflow_subsurface[v] + inflow
-        update_reservoir_model!(reservoir, v, net_inflow, dt)
-        river_v.q[i] = reservoir.variables.outflow[v]
-        # average river discharge (here accumulated for model timestep dt)
-        river_v.q_cumulative[i] += river_v.q[i] * dt
+        # use "upwind" option for water depth estimation at edge to avoid large q
+        # fluctuations (compared to "max" option)
+        floodplain_v.water_depth_at_edge[i] =
+            max(river_v.zs_at_edge[i] - floodplain_p.zb_at_edge[i], 0.0)
+
+        i1, i2 = interpolation_indices(
+            floodplain_v.water_depth_at_edge[i],
+            @view profile.depth[:]
+        )
+        a_src = compute_floodplain_flow_area(
+            profile,
+            floodplain_v.water_depth_at_edge[i],
+            i_src,
+            i1,
+            i2,
+        )
+        a_dst = compute_floodplain_flow_area(
+            profile,
+            floodplain_v.water_depth_at_edge[i],
+            i_dst,
+            i1,
+            i2,
+        )
+        flow_area_at_edge = min(a_src, a_dst)
+
+        hydraulic_radius_at_edge = ifelse(
+            a_src < a_dst,
+            a_src / compute_wetted_perimeter(
+                profile,
+                floodplain_v.water_depth_at_edge[i],
+                i_src,
+                i1,
+            ),
+            a_dst / compute_wetted_perimeter(
+                profile,
+                floodplain_v.water_depth_at_edge[i],
+                i_dst,
+                i1,
+            ),
+        )
+
+        floodplain_v.q[i] = ifelse(
+            flow_area_at_edge > 1.0e-05,
+            manning_flow(
+                floodplain_p.mannings_n_at_edge[i],
+                hydraulic_radius_at_edge,
+                floodplain_p.slope_at_edge[i],
+                flow_area_at_edge,
+            ),
+            0.0,
+        )
+
+        # limit floodplain q in case water is not available
+        floodplain_v.q[i] = min(floodplain_v.q[i], floodplain_v.storage[i_src]/dt)
+
+        # average floodplain discharge (here accumulated for model timestep Δt)
+        floodplain_v.q_cumulative[i] += floodplain_v.q[i] * dt
     end
     return nothing
 end
 
+update_floodplain_flow!(
+    model::RiverFlowModel{T, F},
+    domain::DomainRiver,
+    dt::Float64,
+) where {T <: AbstractStaggeredRoutingMethod, F <: Nothing} = nothing
+
 """
-Update reservoir boundary conditions for the local inertial river flow model.
+Update reservoir boundary conditions for a river flow model on a staggered grid.
 """
 function update_bc_reservoir_model!(
     reservoir_model::ReservoirModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::Domain,
     dt::Float64,
 )
@@ -566,19 +659,21 @@ function update_bc_reservoir_model!(
     end
     return nothing
 end
+
 update_bc_reservoir_model!(
     reservoir_model::Nothing,
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::Domain,
     dt::Float64,
 ) = nothing
 
 """
-Update floodplain water depth and storage.
+Update river and floodplain water depth and storage for river and floodplain models on a
+staggered grid.
 """
 function update_water_depth_and_storage!(
     floodplain_model::AbstractFloodPlainModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::DomainRiver,
     dt::Float64,
 )
@@ -601,7 +696,7 @@ function update_water_depth_and_storage!(
         storage_total = river_v.storage[i] + floodplain_v.storage[i]
         if storage_total > river_p.bankfull_storage[i]
             flood_storage = storage_total - river_p.bankfull_storage[i]
-            h = flood_depth(floodplain_p.profile, flood_storage, flow_length[i], i)
+            h = compute_flood_depth(floodplain_p.profile, flood_storage, flow_length[i], i)
             river_v.h[i] = river_p.bankfull_depth[i] + h
             river_v.storage[i] = river_v.h[i] * flow_width[i] * flow_length[i]
             floodplain_v.storage[i] = max(storage_total - river_v.storage[i], 0.0)
@@ -617,16 +712,16 @@ function update_water_depth_and_storage!(
 end
 update_water_depth_and_storage!(
     floodplain_model::Nothing,
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::DomainRiver,
     dt::Float64,
 ) = nothing
 
 """
-Update water depth and storage for river.
+Update water depth and storage for river flow model on a staggered grid.
 """
 function update_water_depth_and_storage!(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::DomainRiver,
     dt::Float64,
 )
@@ -663,9 +758,9 @@ function update_water_depth_and_storage!(
     return nothing
 end
 
-"Update local inertial river flow model `LocalInertialRiverFlowModel` for flow_area single timestep"
-function local_inertial_river_update!(
-    river_flow_model::LocalInertialRiverFlowModel,
+"Update river flow model on a staggered grid for a single timestep"
+function staggered_scheme_river_update!(
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::Domain,
     dt::Float64,
     update_h::Bool,
@@ -699,11 +794,11 @@ function local_inertial_river_update!(
 end
 
 """
-Update local inertial river flow model `LocalInertialRiverFlow` for flow_area single timestep `dt`. An adaptive
-timestepping method is used (computing flow_area sub timestep `dt_s`).
+Update river flow model on a staggered grid for a single timestep `dt`. An adaptive
+timestepping method is used (computing a sub timestep `dt_s`).
 """
 function update_river_flow_model!(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     domain::Domain,
     clock::Clock,
     dt::Float64;
@@ -724,7 +819,7 @@ function update_river_flow_model!(
     while t < dt
         dt_s = stable_timestep(river_flow_model, flow_length)
         dt_s = check_timestepsize(dt_s, t, dt)
-        local_inertial_river_update!(river_flow_model, domain, dt_s, update_h)
+        staggered_scheme_river_update!(river_flow_model, domain, dt_s, update_h)
         t += dt_s
     end
     average_flow_vars!(river_flow_model, dt)
@@ -743,7 +838,7 @@ function update_river_flow_model!(
 end
 
 "Struct to store local inertial overland flow model variables"
-@with_kw struct LocalInertialOverlandFlowVariables
+@with_kw struct LocalInertialOverlandFlowVariables <: AbstractOverlandFlowVariables
     n::Int
     # flow in y direction at edge at previous time step [m³ s⁻¹]
     qy0::Vector{Float64} = zeros(n + 1)
@@ -770,13 +865,13 @@ end
 end
 
 "Struct to store local inertial overland flow model parameters"
-@with_kw struct LocalInertialOverlandFlowParameters
+@with_kw struct LocalInertialOverlandFlowParameters <: AbstractOverlandFlowParameters
     # number of cells [-]
     n::Int
     # effective flow width x direction at edge (floodplain) [m]
-    xwidth::Vector{Float64}
+    xwidth_at_edge::Vector{Float64}
     # effective flow width y direction at edge (floodplain) [m]
-    ywidth::Vector{Float64}
+    ywidth_at_edge::Vector{Float64}
     # acceleration due to gravity [m s⁻²]
     g::Float64 = 9.80665
     # weighting factor (de Almeida et al., 2012) [-]
@@ -784,18 +879,18 @@ end
     # depth threshold for calculating flow [m]
     h_thresh::Float64
     # maximum cell elevation at edge [m] (x direction)
-    zx_max::Vector{Float64}
+    zx_max_at_edge::Vector{Float64}
     # maximum cell elevation at edge [m] (y direction)
-    zy_max::Vector{Float64}
+    zy_max_at_edge::Vector{Float64}
     # Manning's roughness squared at edge [(s m-1/3)2]
-    mannings_n_sq::Vector{Float64}
+    mannings_n_sq_at_edge::Vector{Float64}
     # elevation [m] of cell
     z::Vector{Float64}
     # if true a check is performed if froude number > 1.0 (algorithm is modified) [-]
     froude_limit::Bool
 end
 
-"Initialize shallow water overland flow model parameters"
+"Initialize local inertial overland flow model parameters"
 function LocalInertialOverlandFlowParameters(
     dataset::NCDataset,
     config::Config,
@@ -831,16 +926,16 @@ function LocalInertialOverlandFlowParameters(
     )
     n = length(domain.land.network.indices)
 
-    zx_max = zeros(n)
-    zy_max = zeros(n)
+    zx_max_at_edge = zeros(n)
+    zy_max_at_edge = zeros(n)
     for i in 1:n
-        idx_right = edge_indices.idx_right[i]
-        if idx_right <= n
-            zx_max[i] = max(elevation[i], elevation[idx_right])
+        ind_x_up = edge_indices.ind_x_up[i]
+        if ind_x_up <= n
+            zx_max_at_edge[i] = max(elevation[i], elevation[ind_x_up])
         end
-        idx_up = edge_indices.idx_up[i]
-        if idx_up <= n
-            zy_max[i] = max(elevation[i], elevation[idx_up])
+        ind_y_up = edge_indices.ind_y_up[i]
+        if ind_y_up <= n
+            zy_max_at_edge[i] = max(elevation[i], elevation[ind_y_up])
         end
     end
 
@@ -852,13 +947,13 @@ function LocalInertialOverlandFlowParameters(
     set_effective_flowwidth!(we_x, we_y, domain)
     parameters = LocalInertialOverlandFlowParameters(;
         n,
-        xwidth = we_x,
-        ywidth = we_y,
+        xwidth_at_edge = we_x,
+        ywidth_at_edge = we_y,
         theta,
         h_thresh = waterdepth_threshold,
-        zx_max,
-        zy_max,
-        mannings_n_sq = mannings_n .* mannings_n,
+        zx_max_at_edge,
+        zy_max_at_edge,
+        mannings_n_sq_at_edge = mannings_n .* mannings_n,
         z = elevation,
         froude_limit,
     )
@@ -866,22 +961,18 @@ function LocalInertialOverlandFlowParameters(
 end
 
 "Struct to store local inertial overland flow model boundary conditions"
-@with_kw struct LocalInertialOverlandFlowBC
+@with_kw struct LocalInertialOverlandFlowBC <: AbstractOverlandFlowBC
     n::Int
     # runoff from hydrological model [m³ s⁻¹]
     runoff::Vector{Float64} = zeros(n)
 end
 
-"Local inertial overland flow model using the local inertial method"
-@with_kw struct LocalInertialOverlandFlowModel <: AbstractOverlandFlowModel
-    timestepping::TimeStepping
-    boundary_conditions::LocalInertialOverlandFlowBC
-    parameters::LocalInertialOverlandFlowParameters
-    variables::LocalInertialOverlandFlowVariables
-end
-
 "Initialize local inertial overland flow model"
-function LocalInertialOverlandFlowModel(dataset::NCDataset, config::Config, domain::Domain)
+function init_local_inertial_overland_flow(
+    dataset::NCDataset,
+    config::Config,
+    domain::Domain,
+)
     alpha_coefficient = config.model.land_local_inertial_flow__alpha_coefficient # stability coefficient for model time step (0.2-0.7)
     timestepping = TimeStepping(; alpha_coefficient)
 
@@ -889,30 +980,33 @@ function LocalInertialOverlandFlowModel(dataset::NCDataset, config::Config, doma
     boundary_conditions = LocalInertialOverlandFlowBC(; n)
     parameters = LocalInertialOverlandFlowParameters(dataset, config, domain)
     variables = LocalInertialOverlandFlowVariables(; n)
+    routing_method = LocalInertial()
 
-    overland_flow_model = LocalInertialOverlandFlowModel(;
+    overland_flow_model = OverlandFlowModel(;
         timestepping,
         boundary_conditions,
         parameters,
         variables,
+        routing_method,
     )
 
     return overland_flow_model
 end
 
 """
-    stable_timestep(river_flow_model::LocalInertialRiverFlowModel, flow_length::Vector{Float64})
-    stable_timestep(overland_flow_model::LocalInertialOverlandFlowModel, parameters::LandParameters)
+    stable_timestep(river_flow_model::RiverFlowModel{<:LocalInertial}, flow_length::Vector{Float64})
+    stable_timestep(overland_flow_model::OverlandFlowModel{<:LocalInertial}, parameters::LandParameters)
 
-Compute flow_area stable timestep size for the local inertial approach, based on Bates et al. (2010).
+Compute stable timestep size for the local inertial approach, based on Bates et al. (2010).
 
 dt = α * (Δx / sqrt(g max(h))
 """
 function stable_timestep(
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     flow_length::Vector{Float64},
 )
     dt_min = Inf
+    dt_min_default = 60.0
     (; alpha_coefficient) = river_flow_model.timestepping
     (; n) = river_flow_model.parameters
     (; h) = river_flow_model.variables
@@ -921,15 +1015,16 @@ function stable_timestep(
             alpha_coefficient * flow_length[i] / sqrt(GRAVITATIONAL_ACCELERATION * h[i])
         dt_min = min(dt, dt_min)
     end
-    dt_min = isinf(dt_min) ? 60.0 : dt_min
+    dt_min = isinf(dt_min) ? dt_min_default : dt_min
     return dt_min
 end
 
 function stable_timestep(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     parameters::LandParameters,
 )
     dt_min = Inf
+    dt_min_default = 60.0
     (; alpha_coefficient) = overland_flow_model.timestepping
     (; n) = overland_flow_model.parameters
     (; x_length, y_length, river_location) = parameters
@@ -943,16 +1038,47 @@ function stable_timestep(
         end
         dt_min = min(dt, dt_min)
     end
-    dt_min = isinf(dt_min) ? 60.0 : dt_min
+    dt_min = isinf(dt_min) ? dt_min_default : dt_min
+    return dt_min
+end
+
+"""
+    stable_timestep(river_flow_model::RiverFlowModel{<:ManningStaggered}, _)
+
+Compute stable timestep size for a river flow model using the Manning's equation on staggered grid.
+"""
+function stable_timestep(river_flow_model::RiverFlowModel{<:ManningStaggered}, _)
+    dt_min = Inf
+    dt_min_default = 60.0
+    (; alpha_coefficient) = river_flow_model.timestepping
+    (;
+        active_e,
+        mannings_n_at_edge,
+        flow_length_at_edge,
+        flow_width_at_edge,
+        slope_at_edge,
+    ) = river_flow_model.parameters
+    (; water_depth_at_edge) = river_flow_model.variables
+
+    @batch per = thread reduction = ((min, dt_min),) for i in active_e
+        @fastmath @inbounds h_r =
+            (flow_width_at_edge[i] * water_depth_at_edge[i]) /
+            wetted_perimeter_channel(water_depth_at_edge[i], flow_width_at_edge[i])
+        # estimate wave celerity for wide rectangular channel
+        celerity =
+            cbrt(h_r^2) * sqrt(slope_at_edge[i]) / mannings_n_at_edge[i] / BETA_KINWAVE
+        dt = alpha_coefficient * flow_length_at_edge[i] / celerity
+        dt_min = min(dt, dt_min)
+    end
+    dt_min = isinf(dt_min) ? dt_min_default : dt_min
     return dt_min
 end
 
 """
-Update boundary condition `runoff` overland flow model `LocalInertialOverlandFlowModel` for flow_area
-single timestep.
+Update boundary condition `runoff` local inertial overland flow model for a single timestep.
 """
 function update_bc_overland_flow_model!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     external_models::NamedTuple,
     domain::Domain,
     dt::Float64,
@@ -971,12 +1097,12 @@ function update_bc_overland_flow_model!(
 end
 
 """
-Update subsurface flow contribution to inflow of flow_area reservoir model for flow_area river flow model
-`LocalInertialRiverFlowModel` for flow_area single timestep.
+Update subsurface flow contribution to inflow of a reservoir model for a river flow model on
+a staggered grid for a single timestep.
 """
 function update_inflow!(
     reservoir_model::ReservoirModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    river_flow_model::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     subsurface_flow_model::AbstractSubsurfaceFlowModel,
     network::NetworkReservoir,
 )
@@ -988,17 +1114,17 @@ function update_inflow!(
 end
 update_inflow!(
     ::Nothing,
-    ::LocalInertialRiverFlowModel,
+    ::RiverFlowModel{<:AbstractStaggeredRoutingMethod},
     ::AbstractSubsurfaceFlowModel,
     ::NetworkReservoir,
 ) = nothing
 
 """
-Helper function to set flow variables of the `LocalInertialOverlandFlowModel` model to zero. This
-is done at the start of each simulation timestep, during the timestep the total (weighted)
-sum is computed from values at each sub timestep.
+Helper function to set flow variables of the local inertial overland flow model to zero.
+This is done at the start of each simulation timestep, during the timestep the total
+(weighted) sum is computed from values at each sub timestep.
 """
-function set_flow_vars!(overland_flow_model::LocalInertialOverlandFlowModel)
+function set_flow_vars!(overland_flow_model::OverlandFlowModel{<:LocalInertial})
     (; qx_cumulative, qy_cumulative) = overland_flow_model.variables
     qx_cumulative .= 0.0
     qy_cumulative .= 0.0
@@ -1006,11 +1132,11 @@ function set_flow_vars!(overland_flow_model::LocalInertialOverlandFlowModel)
 end
 
 """
-Helper function to compute average flow variables of the `LocalInertialOverlandFlowModel` model.
+Helper function to compute average flow variables of the local inertial overland flow model.
 This is done at the end of each simulation timestep.
 """
 function average_flow_vars!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     dt::Float64,
 )
     (; qx_average, qx_cumulative, qy_average, qy_cumulative) = overland_flow_model.variables
@@ -1021,13 +1147,12 @@ function average_flow_vars!(
 end
 
 """
-Update combined river `LocalInertialRiverFlowModel` and overland flow `LocalInertialOverlandFlowModel`
-models for flow_area single timestep `dt`. An adaptive timestepping method is used (computing flow_area sub
-timestep `dt_s`).
+Update local inertial river and overland flow model for a single timestep `dt`. An adaptive
+timestepping method is used (computing a sub timestep `dt_s`).
 """
 function update_overland_flow_model!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     domain::Domain,
     clock::Clock,
     dt::Float64;
@@ -1035,7 +1160,7 @@ function update_overland_flow_model!(
 )
     (; reservoir) = river_flow_model.boundary_conditions
     (; flow_length) = domain.river.parameters
-    (; parameters) = domain.land
+    land_parameters = domain.land.parameters
 
     set_reservoir_vars!(reservoir)
     update_index_hq!(reservoir, clock)
@@ -1045,13 +1170,13 @@ function update_overland_flow_model!(
     t = 0.0
     while t < dt
         dt_river = stable_timestep(river_flow_model, flow_length)
-        dt_land = stable_timestep(overland_flow_model, parameters)
+        dt_land = stable_timestep(overland_flow_model, land_parameters)
         dt_s = min(dt_river, dt_land)
         dt_s = check_timestepsize(dt_s, t, dt)
 
         local_inertial_update_fluxes!(overland_flow_model, domain, dt_s)
         update_inflow_reservoir!(overland_flow_model, reservoir, domain)
-        local_inertial_river_update!(river_flow_model, domain, dt_s, update_h)
+        staggered_scheme_river_update!(river_flow_model, domain, dt_s, update_h)
         local_inertial_update_water_depth!(
             overland_flow_model,
             river_flow_model,
@@ -1069,11 +1194,12 @@ function update_overland_flow_model!(
 end
 
 """
-Update flow for flow_area single direction in the local inertial overland flow model.
-`is_x_direction`: true for x-direction (idx_right/idx_left), false for y-direction (idx_up/idx_down)
+Update flow for the local inertial overland flow model at edge `i` in a single direction.
+`is_x_direction`: true for x-direction (ind_x_up/ind_x_down), false for y-direction
+(ind_y_up/ind_y_down)
 """
 @inline function update_directional_flow!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     domain::Domain,
     i::Int,
     dt::Float64,
@@ -1086,19 +1212,19 @@ Update flow for flow_area single direction in the local inertial overland flow m
 
     # Select direction-specific parameters based on the boolean flag
     if is_x_direction
-        upstream_idx = indices.idx_right[i]
-        downstream_idx = indices.idx_left[i]
-        width = land_p.ywidth[i]
-        z_max = land_p.zx_max[i]
+        upstream_idx = indices.ind_x_up[i]
+        downstream_idx = indices.ind_x_down[i]
+        width_at_edge = land_p.ywidth_at_edge[i]
+        z_max_at_edge = land_p.zx_max_at_edge[i]
         length_vec = x_length
         q_current = land_v.qx
         q_prev = land_v.qx0
         q_cumulative = land_v.qx_cumulative
     else
-        upstream_idx = indices.idx_up[i]
-        downstream_idx = indices.idx_down[i]
-        width = land_p.xwidth[i]
-        z_max = land_p.zy_max[i]
+        upstream_idx = indices.ind_y_up[i]
+        downstream_idx = indices.ind_y_down[i]
+        width_at_edge = land_p.xwidth_at_edge[i]
+        z_max_at_edge = land_p.zy_max_at_edge[i]
         length_vec = y_length
         q_current = land_v.qy
         q_prev = land_v.qy0
@@ -1107,14 +1233,14 @@ Update flow for flow_area single direction in the local inertial overland flow m
 
     # the effective flow width is zero when the river width exceeds the cell width and
     # floodplain flow is not calculated.
-    if upstream_idx <= land_p.n && width != 0.0
+    if upstream_idx <= land_p.n && width_at_edge != 0.0
         zs_current = land_p.z[i] + land_v.h[i]
         zs_upstream = land_p.z[upstream_idx] + land_v.h[upstream_idx]
-        zs_max = max(zs_current, zs_upstream)
-        water_depth_at_edge = (zs_max - z_max)
+        zs_max_at_edge = max(zs_current, zs_upstream)
+        water_depth_at_edge = (zs_max_at_edge - z_max_at_edge)
 
         if water_depth_at_edge > land_p.h_thresh
-            length = 0.5 * (length_vec[i] + length_vec[upstream_idx]) # can be precalculated
+            length_at_edge = 0.5 * (length_vec[i] + length_vec[upstream_idx]) # can be precalculated
             q_current[i] = local_inertial_flow(
                 land_p.theta,
                 q_prev[i],
@@ -1123,9 +1249,9 @@ Update flow for flow_area single direction in the local inertial overland flow m
                 zs_current,
                 zs_upstream,
                 water_depth_at_edge,
-                width,
-                length,
-                land_p.mannings_n_sq[i],
+                width_at_edge,
+                length_at_edge,
+                land_p.mannings_n_sq_at_edge[i],
                 land_p.froude_limit,
                 dt,
             )
@@ -1145,11 +1271,10 @@ Update flow for flow_area single direction in the local inertial overland flow m
 end
 
 """
-Update fluxes for overland flow `LocalInertialOverlandFlowModel` model for flow_area single timestep
-`dt`.
+Update fluxes for local inertial overland flow model for a single timestep `dt`.
 """
 function local_inertial_update_fluxes!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     domain::Domain,
     dt::Float64,
 )
@@ -1170,12 +1295,11 @@ function local_inertial_update_fluxes!(
 end
 
 """
-Update boundary condition inflow to flow_area reservoir from land `inflow_reservoir` of combined
-river `LocalInertialRiverFlowModel`and overland flow `LocalInertialOverlandFlowModel` models for flow_area
-single timestep.
+Update boundary condition inflow to a reservoir from land of the local inertial river and
+overland flow model for a single timestep.
 """
 function update_inflow_reservoir!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     reservoir_model::Union{ReservoirModel, Nothing},
     domain::Domain,
 )
@@ -1185,21 +1309,22 @@ function update_inflow_reservoir!(
     land_v = overland_flow_model.variables
 
     for (i, j) in enumerate(reservoir_indices)
-        idx_down = indices.idx_down[j]
-        idx_left = indices.idx_left[j]
+        ind_y_down = indices.ind_y_down[j]
+        ind_x_down = indices.ind_x_down[j]
         reservoir_model.boundary_conditions.inflow_overland[i] =
             land_bc.runoff[j] +
-            (land_v.qx[idx_left] - land_v.qx[j] + land_v.qy[idx_down] - land_v.qy[j])
+            (land_v.qx[ind_x_down] - land_v.qx[j] + land_v.qy[ind_y_down] - land_v.qy[j])
     end
     return nothing
 end
 
 """
-Compute storage change for flow_area river cell from fluxes.
+Compute storage change at node index `i` containing a river, from river and overland flow
+fluxes of the local inertial river and overland flow model.
 """
 @inline function compute_river_storage_change(
-    overland_flow_model::LocalInertialOverlandFlowModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     domain::Domain,
     i::Int,
     dt::Float64,
@@ -1209,15 +1334,15 @@ Compute storage change for flow_area river cell from fluxes.
     edges_at_node = domain.river.network.edges_at_node
     river_idx = inds_river[i]
 
-    idx_down = indices.idx_down[i]
-    idx_left = indices.idx_left[i]
+    ind_y_down = indices.ind_y_down[i]
+    ind_x_down = indices.ind_x_down[i]
 
     net_river_flow =
         sum_at(river_flow_model.variables.q, edges_at_node.src[river_idx]) -
         sum_at(river_flow_model.variables.q, edges_at_node.dst[river_idx])
     net_land_flow =
-        overland_flow_model.variables.qx[idx_left] - overland_flow_model.variables.qx[i] +
-        overland_flow_model.variables.qy[idx_down] - overland_flow_model.variables.qy[i]
+        overland_flow_model.variables.qx[ind_x_down] - overland_flow_model.variables.qx[i] +
+        overland_flow_model.variables.qy[ind_y_down] - overland_flow_model.variables.qy[i]
     net_flow =
         net_river_flow + net_land_flow + overland_flow_model.boundary_conditions.runoff[i] -
         river_flow_model.boundary_conditions.abstraction[river_idx]
@@ -1227,12 +1352,13 @@ Compute storage change for flow_area river cell from fluxes.
 end
 
 """
-Compute external inflow for river cells, including negative inflow (abstraction).
+Compute external inflow at river node index `river_idx`, including negative inflow
+(abstraction), of the local inertial river and overland flow model.
 Returns tuple: (inflow, abstraction_to_add)
 """
 @inline function compute_external_inflow(
-    river_flow_model::LocalInertialRiverFlowModel,
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    river_flow_model::RiverFlowModel{<:LocalInertial},
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     i::Int,
     river_idx::Int,
     dt::Float64,
@@ -1263,7 +1389,7 @@ Returns tuple: (river_h, land_h, river_storage)
     total_storage::Float64,
     river_idx::Int,
     i::Int,
-    river::LocalInertialRiverFlowModel,
+    river::RiverFlowModel{<:LocalInertial},
     domain::Domain,
 )
     if total_storage >= river.parameters.bankfull_storage[river_idx]
@@ -1290,31 +1416,33 @@ Returns tuple: (river_h, land_h, river_storage)
 end
 
 """
-Compute storage change for flow_area land cell from horizontal fluxes and runoff.
+Compute storage change for local inertial overland flow model at node index `i` (non-river)
+from horizontal fluxes and runoff.
 """
 @inline function compute_land_storage_change(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     network::NetworkLand,
     i::Int,
     dt::Float64,
 )
     indices = network.edge_indices
-    idx_down = indices.idx_down[i]
-    idx_left = indices.idx_left[i]
+    ind_y_down = indices.ind_y_down[i]
+    ind_x_down = indices.ind_x_down[i]
 
     return (
-        overland_flow_model.variables.qx[idx_left] - overland_flow_model.variables.qx[i] +
-        overland_flow_model.variables.qy[idx_down] - overland_flow_model.variables.qy[i] +
+        overland_flow_model.variables.qx[ind_x_down] - overland_flow_model.variables.qx[i] +
+        overland_flow_model.variables.qy[ind_y_down] - overland_flow_model.variables.qy[i] +
         overland_flow_model.boundary_conditions.runoff[i]
     ) * dt
 end
 
 """
-Update storage and water depth for flow_area single river cell.
+Update storage and water depth for the local inertial overland flow model at node index `i`
+containing a river.
 """
-@inline function update_river_cell_storage_and_depth!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+@inline function update_river_and_land_storage_and_depth!(
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     domain::Domain,
     i::Int,
     dt::Float64,
@@ -1357,10 +1485,11 @@ Update storage and water depth for flow_area single river cell.
 end
 
 """
-Update storage and water depth for flow_area single land cell (non-river).
+Update storage and water depth for local inertial overland flow model at node index `i`
+(non-river).
 """
-@inline function update_land_cell_storage_and_depth!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
+@inline function update_land_storage_and_depth!(
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
     domain::DomainLand,
     i::Int,
     dt::Float64,
@@ -1385,12 +1514,12 @@ Update storage and water depth for flow_area single land cell (non-river).
 end
 
 """
-Update storage and water depth for combined river `LocalInertialRiverFlowModel` and overland flow
-`LocalInertialOverlandFlowModel` models for flow_area single timestep `dt`.
+Update storage and water depth for local inertial river and overland flow model for a single
+timestep `dt`.
 """
 function local_inertial_update_water_depth!(
-    overland_flow_model::LocalInertialOverlandFlowModel,
-    river_flow_model::LocalInertialRiverFlowModel,
+    overland_flow_model::OverlandFlowModel{<:LocalInertial},
+    river_flow_model::RiverFlowModel{<:LocalInertial},
     domain::Domain,
     dt::Float64,
 )
@@ -1400,7 +1529,7 @@ function local_inertial_update_water_depth!(
         if river_location[i]
             # Process river cells (excluding reservoir outlets)
             if !reservoir_outlet[i]
-                update_river_cell_storage_and_depth!(
+                update_river_and_land_storage_and_depth!(
                     overland_flow_model,
                     river_flow_model,
                     domain,
@@ -1410,296 +1539,8 @@ function local_inertial_update_water_depth!(
             end
         else
             # Process land cells (non-river)
-            update_land_cell_storage_and_depth!(overland_flow_model, domain.land, i, dt)
-        end
-    end
-    return nothing
-end
-
-"""
-    FloodPlainProfile
-
-Floodplain `storage` is flow_area function of `depth` (flood depth intervals). Based on the
-cumulative floodplain `storage` flow_area floodplain profile as flow_area function of `flood_depth` is
-derived with floodplain area `flow_area` (cumulative) and wetted perimeter radius `wetted_perimeter` (cumulative).
-"""
-@with_kw struct FloodPlainProfile
-    depth::Vector{Float64}        # Flood depth [m]
-    storage::Matrix{Float64}      # Flood storage (cumulative) [m³]
-    width::Matrix{Float64}        # Flood width [m]
-    flow_area::Matrix{Float64}            # Flow area (cumulative) [m²]
-    wetted_perimeter::Matrix{Float64}            # Wetted perimeter (cumulative) [m]
-end
-
-"Initialize floodplain profile `FloodPlainProfile`"
-function FloodPlainProfile(
-    dataset::NCDataset,
-    config::Config,
-    domain::DomainRiver,
-    index_pit::Vector{Int},
-)
-    (; indices) = domain.network
-    (; flow_width, flow_length) = domain.parameters
-    storage = ncread(
-        dataset,
-        config,
-        "floodplain_water__sum_of_volume_per_depth",
-        Routing;
-        sel = indices,
-    )
-    n = length(indices)
-
-    # for convenience (interpolation) flood depth 0.0 m is added, with associated area (a),
-    # storage, width (river width) and wetted perimeter (p).
-    storage = vcat(fill(Float64(0), n)', storage)
-    start_storage = storage
-    flood_depths = Float64.(dataset["flood_depth"][:])
-    pushfirst!(flood_depths, 0.0)
-    n_depths = length(flood_depths)
-
-    wetted_perimeter = zeros(n_depths, n)
-    flow_area = zeros(n_depths, n)
-    segment_storage = zeros(n_depths, n)
-    width = zeros(n_depths, n)
-    width[1, :] = flow_width[1:n]
-
-    # determine flow area (a), width and wetted perimeter (p) FloodPlainModel
-    h = diff(flood_depths)
-    incorrect_vol = 0
-    riv_cells = 0
-    error_vol = 0
-    for i in 1:n
-        riv_cell = 0
-        diff_storage = diff(storage[:, i])
-
-        for j in 1:(n_depths - 1)
-            # assume rectangular shape of flood depth segment
-            width[j + 1, i] = diff_storage[j] / (h[j] * flow_length[i])
-            # check provided flood storage (floodplain width should be constant or increasing
-            # as a function of flood depth)
-            if width[j + 1, i] < width[j, i]
-                # raise warning only if difference is larger than rounding error of 0.01 m³
-                if ((width[j, i] - width[j + 1, i]) * h[j] * flow_length[i]) > 0.01
-                    incorrect_vol += 1
-                    riv_cell = 1
-                    error_vol =
-                        error_vol +
-                        ((width[j, i] - width[j + 1, i]) * h[j] * flow_length[i])
-                end
-                width[j + 1, i] = width[j, i]
-            end
-            flow_area[j + 1, i] = width[j + 1, i] * h[j]
-            wetted_perimeter[j + 1, i] = (width[j + 1, i] - width[j, i]) + 2.0 * h[j]
-            segment_storage[j + 1, i] = flow_area[j + 1, i] * flow_length[i]
-            if j == 1
-                # for interpolation wetted perimeter at flood depth 0.0 is required
-                wetted_perimeter[j, i] = wetted_perimeter[j + 1, i] - 2.0 * h[j]
-            end
+            update_land_storage_and_depth!(overland_flow_model, domain.land, i, dt)
         end
-
-        wetted_perimeter[2:end, i] = cumsum(wetted_perimeter[2:end, i])
-        flow_area[:, i] = cumsum(flow_area[:, i])
-        storage[:, i] = cumsum(segment_storage[:, i])
-
-        riv_cells += riv_cell
-    end
-
-    if incorrect_vol > 0
-        perc_riv_cells = round(100.0 * (riv_cells / n); digits = 2)
-        perc_error_vol = round(100.0 * (error_vol / sum(start_storage[end, :])); digits = 2)
-        @warn string(
-            "The provided storage of $incorrect_vol rectangular floodplain schematization",
-            " segments for $riv_cells river cells ($perc_riv_cells % of total river cells)",
-            " is not correct and has been increased with $perc_error_vol % of provided storage.",
-        )
-    end
-
-    # set floodplain parameters for ghost points
-    storage = hcat(storage, storage[:, index_pit])
-    width = hcat(width, width[:, index_pit])
-    flow_area = hcat(flow_area, flow_area[:, index_pit])
-    wetted_perimeter = hcat(wetted_perimeter, wetted_perimeter[:, index_pit])
-
-    # initialize floodplain profile parameters
-    profile = FloodPlainProfile(;
-        storage,
-        width,
-        depth = flood_depths,
-        flow_area,
-        wetted_perimeter,
-    )
-    return profile
-end
-
-"Struct to store floodplain flow model parameters"
-@with_kw struct FloodPlainParameters
-    profile::FloodPlainProfile          # floodplain profile
-    mannings_n::Vector{Float64}         # manning's roughness [s m-1/3]
-    mannings_n_sq::Vector{Float64}      # manning's roughness squared at edge [(s m-1/3)2]
-    zb_max::Vector{Float64}             # maximum bankfull elevation at edge [m]
-end
-
-"Initialize floodplain flow model parameters"
-function FloodPlainParameters(
-    dataset::NCDataset,
-    config::Config,
-    domain::DomainRiver,
-    zb_floodplain::Vector{Float64},
-    index_pit::Vector{Int},
-)
-    (; indices, nodes_at_edge, graph) = domain.network
-    (; flow_length) = domain.parameters
-    n_edges = ne(graph)
-    profile = FloodPlainProfile(dataset, config, domain, index_pit)
-
-    mannings_n = ncread(
-        dataset,
-        config,
-        "floodplain_water_flow__manning_n_parameter",
-        Routing;
-        sel = indices,
-    )
-    # manning roughness at edges
-    append!(mannings_n, mannings_n[index_pit]) # copy to ghost nodes
-    mannings_n_sq = fill(Float64(0), n_edges)
-    zb_max = fill(Float64(0), n_edges)
-    for i in 1:n_edges
-        src_node = nodes_at_edge.src[i]
-        dst_node = nodes_at_edge.dst[i]
-        mannings_n_i =
-            (
-                mannings_n[dst_node] * flow_length[dst_node] +
-                mannings_n[src_node] * flow_length[src_node]
-            ) / (flow_length[dst_node] + flow_length[src_node])
-        mannings_n_sq[i] = mannings_n_i * mannings_n_i
-        zb_max[i] = max(zb_floodplain[src_node], zb_floodplain[dst_node])
-    end
-    parameters = FloodPlainParameters(profile, mannings_n, mannings_n_sq, zb_max)
-    return parameters
-end
-
-"Struct to store floodplain flow model variables"
-@with_kw struct FloodPlainVariables
-    n::Int
-    n_edges::Int
-    # storage [m³]
-    storage::Vector{Float64} = zeros(n)
-    # water depth [m]
-    h::Vector{Float64}
-    # error storage [m³]
-    error::Vector{Float64} = zeros(n)
-    # flow area at edge [m²]
-    flow_area::Vector{Float64} = zeros(n_edges)
-    # hydraulic radius at edge [m]
-    hydraulic_radius::Vector{Float64} = zeros(n_edges)
-    # water depth at edge [m]
-    water_depth_at_edge::Vector{Float64} = zeros(n_edges)
-    # discharge at edge at previous time step
-    q_previous::Vector{Float64} = zeros(n_edges)
-    # discharge at edge  [m³ s⁻¹]
-    q::Vector{Float64} = zeros(n_edges)
-    # cumulative river discharge at edge [m³] for model timestep dt
-    q_cumulative::Vector{Float64} = zeros(n_edges)
-    # average river discharge at edge [m³ s⁻¹] for model timestep dt
-    q_average::Vector{Float64} = zeros(n_edges)
-    # edge index with `hf` [-] above depth threshold
-    hf_index::Vector{Int} = zeros(Int, n_edges)
-end
-
-"Floodplain flow model"
-@with_kw struct FloodPlainModel <: AbstractFloodPlainModel
-    parameters::FloodPlainParameters
-    variables::FloodPlainVariables
-end
-
-"Determine the initial floodplain storage"
-function initialize_storage!(river, domain::Domain, nriv::Int)
-    (; flow_width, flow_length) = domain.river.parameters
-    (; floodplain) = river
-    (; profile) = floodplain.parameters
-    for i in 1:nriv
-        i1, i2 = interpolation_indices(floodplain.variables.h[i], profile.depth)
-        flow_area = flow_area(
-            profile.width[i2, i],
-            profile.flow_area[i1, i],
-            profile.depth[i1],
-            floodplain.variables.h[i],
-        )
-        flow_area = max(flow_area - (flow_width[i] * floodplain.variables.h[i]), 0.0)
-        floodplain.variables.storage[i] = flow_length[i] * flow_area
     end
     return nothing
 end
-
-"helper function to get interpolation indices"
-function interpolation_indices(x, v::AbstractVector)
-    i1 = 1
-    for i in eachindex(v)
-        if v[i] <= x
-            i1 = i
-        end
-    end
-    if i1 == length(v)
-        i2 = i1
-    else
-        i2 = i1 + 1
-    end
-    return i1, i2
-end
-
-"""
-    flow_area(width, area, depth, h)
-
-Compute floodplain flow area based on flow depth `h` and floodplain `depth`, `area` and
-`width` of flow_area floodplain profile.
-"""
-function flow_area(width, area, depth, h)
-    dh = h - depth  # depth at i1
-    area = area + (width * dh) # area at i1, width at i2
-    return area
-end
-
-"""
-    function wetted_perimeter(wetted_perimeter, depth, h)
-
-Compute floodplain wetted perimeter based on flow depth `h` and floodplain `depth` and
-wetted perimeter `wetted_perimeter` of flow_area floodplain profile.
-"""
-function wetted_perimeter(wetted_perimeter, depth, h)
-    dh = h - depth # depth at i1
-    wetted_perimeter += 2.0 * dh # p at i1
-    return wetted_perimeter
-end
-
-"Compute flood depth by interpolating flood storage `flood_storage` using flood depth intervals."
-function flood_depth(
-    profile::FloodPlainProfile,
-    flood_storage::Float64,
-    flow_length::Float64,
-    i::Int,
-)
-    i1, i2 = interpolation_indices(flood_storage, @view profile.storage[:, i])
-    ΔA = (flood_storage - profile.storage[i1, i]) / flow_length
-    dh = ΔA / profile.width[i2, i]
-    flood_depth = profile.depth[i1] + dh
-    return flood_depth
-end
-
-"Initialize floodplain geometry and `FloodPlainModel` variables and parameters"
-function FloodPlainModel(
-    dataset::NCDataset,
-    config::Config,
-    domain::DomainRiver,
-    zb_floodplain::Vector{Float64},
-)
-    (; indices, local_drain_direction, graph) = domain.network
-    n = length(indices)
-    index_pit = findall(x -> x == 5, local_drain_direction)
-    parameters = FloodPlainParameters(dataset, config, domain, zb_floodplain, index_pit)
-    h = zeros(n + length(index_pit))
-    n_edges = ne(graph)
-    variables = FloodPlainVariables(; n, n_edges, h)
-
-    floodplain = FloodPlainModel(; parameters, variables)
-    return floodplain
-end
diff --git a/Wflow/src/routing/utils.jl b/Wflow/src/routing/utils.jl
index fdbffb660..b590473a0 100644
--- a/Wflow/src/routing/utils.jl
+++ b/Wflow/src/routing/utils.jl
@@ -12,7 +12,7 @@ function flowgraph(ldd::AbstractVector, indices::AbstractVector, PCR_DIR::Abstra
     for (from_node, from_index) in enumerate(indices)
         ldd_val = ldd[from_node]
         # skip pits to prevent cycles
-        ldd_val == 5 && continue
+        ldd_val == LDD_PIT && continue
         to_index = from_index + PCR_DIR[ldd_val]
         # find the node id of the downstream cell
         to_node = searchsortedfirst(indices, to_index)
@@ -166,3 +166,37 @@ function flux_in!(flux_in, flux, network)
     end
     return nothing
 end
+
+function compute_value_at_edge(v, nodes_at_edge, n_edges, func::Function)
+    x = zeros(n_edges)
+    for i in 1:n_edges
+        src_node = nodes_at_edge.src[i]
+        dst_node = nodes_at_edge.dst[i]
+        x[i] = func((v[src_node], v[dst_node]))
+    end
+    return x
+end
+
+function compute_mannings_n_at_edge(mannings_n, flow_length, nodes_at_edge, n_edges)
+    mannings_n_at_edge = zeros(n_edges)
+    for i in 1:n_edges
+        src_node = nodes_at_edge.src[i]
+        dst_node = nodes_at_edge.dst[i]
+        mannings_n_at_edge[i] =
+            (
+                mannings_n[dst_node] * flow_length[dst_node] +
+                mannings_n[src_node] * flow_length[src_node]
+            ) / (flow_length[dst_node] + flow_length[src_node])
+    end
+    return mannings_n_at_edge
+end
+
+function compute_slope_at_edge(elev, length_at_edge, nodes_at_edge, n_edges)
+    slope = zeros(n_edges)
+    for i in 1:n_edges
+        src_node = nodes_at_edge.src[i]
+        dst_node = nodes_at_edge.dst[i]
+        slope[i] = max((elev[src_node] - elev[dst_node]) / length_at_edge[i], 0.00001)
+    end
+    return slope
+end
diff --git a/Wflow/src/sbm_model.jl b/Wflow/src/sbm_model.jl
index 68db24817..d70015abc 100644
--- a/Wflow/src/sbm_model.jl
+++ b/Wflow/src/sbm_model.jl
@@ -106,7 +106,13 @@ function set_states!(model::AbstractModel{<:Union{SbmModel, SbmGwfModel}})
     land_v = routing.overland_flow.variables
     river_v = routing.river_flow.variables
 
-    (; land_routing, cold_start__flag, reservoir__flag, floodplain_1d__flag) = config.model
+    (;
+        river_routing,
+        land_routing,
+        cold_start__flag,
+        reservoir__flag,
+        floodplain_1d__flag,
+    ) = config.model
 
     # read and set states in model object if cold_start=false
     if !cold_start__flag
@@ -165,7 +171,7 @@ function set_states!(model::AbstractModel{<:Union{SbmModel, SbmGwfModel}})
         river_v.storage[1:nriv] .=
             river_v.h[1:nriv] .* flow_width[1:nriv] .* flow_length[1:nriv]
 
-        if floodplain_1d__flag
+        if floodplain_1d__flag && river_routing != RoutingType.kinematic_wave
             initialize_storage!(routing.river_flow, domain, nriv)
         end
 
diff --git a/Wflow/src/standard_name/standard_name_routing.jl b/Wflow/src/standard_name/standard_name_routing.jl
index 45568408b..bcb33c98b 100644
--- a/Wflow/src/standard_name/standard_name_routing.jl
+++ b/Wflow/src/standard_name/standard_name_routing.jl
@@ -233,8 +233,8 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         allow_dynamic_input = true,
         tags = [
             :kinematic_wave_river_static_cyclic_forcing_input,
-            :local_inertial_river_static_cyclic_forcing_input,
-            :local_inertial_river_input,
+            :staggered_grid_river_static_cyclic_forcing_input,
+            :staggered_grid_river_input,
         ],
     ),
     "river__slope" => ParameterMetadata(;
@@ -251,8 +251,8 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         tags = [
             :kinematic_wave_river_state,
             :kinematic_wave_river_output,
-            :local_inertial_river_state,
-            :local_inertial_river_output,
+            :staggered_grid_river_state,
+            :staggered_grid_river_output,
         ],
     ),
     "river_water__depth" => ParameterMetadata(;
@@ -381,26 +381,26 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         description = "Part of subsurface flow that flows to the river",
         tags = [:kinematic_lateral_subsurface_output],
     ),
-    ## Local inertial
+    ## Staggered grid (local inertial and Manning's equation)
     ### River flow
     #### Input
     "model_boundary_condition_river__length" => ParameterMetadata(;
         unit = Unit(; m = 1),
         default = 1.0e4,
         description = "Boundary condition river length downstream river outlets",
-        tags = [:local_inertial_river_input],
+        tags = [:staggered_grid_river_input],
     ),
     "model_boundary_condition_river_bank_water__depth" => ParameterMetadata(;
         default = 0,
         unit = Unit(; m = 1),
         description = "Boundary condition bankfull depth downstream river outlets",
-        tags = [:local_inertial_river_input],
+        tags = [:staggered_grid_river_input],
     ),
     "river_bank_water__elevation" => ParameterMetadata(;
         unit = Unit(; m = 1),
         fill = 0.0,
         description = "Bankfull elevation of the river",
-        tags = [:local_inertial_river_input],
+        tags = [:staggered_grid_river_input],
     ),
     ### 1D floodplain flow
     #### Input
@@ -408,14 +408,26 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         unit = Unit(; m = 3),
         dimname = :flood_depth,
         description = "Floodplain profile (cumulative volume per flood depth)",
-        tags = [:local_inertial_floodplain_1D_flow_input],
+        tags = [
+            :staggered_grid_floodplain_1D_flow_input,
+            :kinematic_wave_floodplain_1D_flow_input,
+        ],
     ),
     "floodplain_water_flow__manning_n_parameter" => ParameterMetadata(;
         lens = @optic(_.routing.river_flow.floodplain.parameters.mannings_n),
         unit = Unit(; s = 1, m = -1 // 3),
         default = 0.072,
         description = "Manning's roughness",
-        tags = [:local_inertial_floodplain_1D_flow_input],
+        tags = [
+            :staggered_grid_floodplain_1D_flow_input,
+            :kinematic_wave_floodplain_1D_flow_input,
+        ],
+    ),
+    "floodplain__slope" => ParameterMetadata(;
+        lens = @optic(_.routing.river_flow.floodplain.parameters.slope),
+        unit = Unit(; m = (1, 1)),
+        description = "Floodplain slope",
+        tags = [:kinematic_wave_floodplain_1D_flow_input],
     ),
     #### States
     "floodplain_water__instantaneous_volume_flow_rate" => ParameterMetadata(;
@@ -423,8 +435,9 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         unit = Unit(; m = 3, s = -1),
         description = "Floodplain discharge",
         tags = [
-            :local_inertial_floodplain_1D_flow_state,
-            :local_inertial_floodplain_1D_flow_output,
+            :staggered_grid_floodplain_1D_flow_state,
+            :staggered_grid_floodplain_1D_flow_output,
+            :kinematic_wave_floodplain_1D_flow_output,
         ],
     ),
     "floodplain_water__depth" => ParameterMetadata(;
@@ -432,22 +445,29 @@ const routing_standard_name_map = OrderedDict{String, ParameterMetadata}(
         unit = Unit(; m = 1),
         description = "Floodplain water depth",
         tags = [
-            :local_inertial_floodplain_1D_flow_state,
-            :local_inertial_floodplain_1D_flow_output,
+            :staggered_grid_floodplain_1D_flow_state,
+            :staggered_grid_floodplain_1D_flow_output,
         ],
     ),
-    #### Output
     "floodplain_water__volume" => ParameterMetadata(;
         lens = @optic(_.routing.river_flow.floodplain.variables.storage),
         unit = Unit(; m = 3),
         description = "Floodplain water volume",
-        tags = [:local_inertial_floodplain_1D_flow_output],
+        tags = [
+            :staggered_grid_floodplain_1D_flow_output,
+            :kinematic_wave_floodplain_1D_flow_state,
+            :kinematic_wave_floodplain_1D_flow_output,
+        ],
     ),
+    #### Output
     "floodplain_water__volume_flow_rate" => ParameterMetadata(;
         lens = @optic(_.routing.river_flow.floodplain.variables.q_average),
         unit = Unit(; m = 3, s = -1),
         description = "Floodplain discharge",
-        tags = [:local_inertial_floodplain_1D_flow_output],
+        tags = [
+            :staggered_grid_floodplain_1D_flow_output,
+            :kinematic_wave_floodplain_1D_flow_output,
+        ],
     ),
     ### Overland flow
     #### Input
diff --git a/Wflow/src/states.jl b/Wflow/src/states.jl
index 8fcefc7d0..a7cb2f087 100644
--- a/Wflow/src/states.jl
+++ b/Wflow/src/states.jl
@@ -124,9 +124,16 @@ function extract_required_states(config::Config)
     end
 
     # Floodplain states
-    floodplain_states =
-        do_floodplains ?
-        get_states_by_tag(Routing, :local_inertial_floodplain_1D_flow_state) : String[]
+    floodplain_states = if !do_floodplains
+        String[]
+    else
+        if config.model.river_routing == RoutingType.local_inertial ||
+           config.model.river_routing == RoutingType.manning_staggered
+            get_states_by_tag(Routing, :staggered_grid_floodplain_1D_flow_state)
+        else
+            get_states_by_tag(Routing, :kinematic_wave_floodplain_1D_flow_state)
+        end
+    end
 
     # Reservoir states
     reservoir_states =
diff --git a/Wflow/src/utils.jl b/Wflow/src/utils.jl
index dc2adeab6..03db893b3 100644
--- a/Wflow/src/utils.jl
+++ b/Wflow/src/utils.jl
@@ -625,44 +625,44 @@ function set_effective_flowwidth!(
             we_x[idx] = reservoir_outlet[v] ? 0.0 : max(we_x[idx] - 0.5 * w, 0.0)
             we_y[idx] = reservoir_outlet[v] ? 0.0 : max(we_y[idx] - 0.5 * w, 0.0)
         elseif dir == CartesianIndex(-1, -1)
-            if edge_indices.idx_left[idx] <= n
-                we_y[edge_indices.idx_left[idx]] =
+            if edge_indices.ind_x_down[idx] <= n
+                we_y[edge_indices.ind_x_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_y[edge_indices.idx_left[idx]] - 0.5 * w, 0.0)
+                    max(we_y[edge_indices.ind_x_down[idx]] - 0.5 * w, 0.0)
             end
-            if edge_indices.idx_down[idx] <= n
-                we_x[edge_indices.idx_down[idx]] =
+            if edge_indices.ind_y_down[idx] <= n
+                we_x[edge_indices.ind_y_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_x[edge_indices.idx_down[idx]] - 0.5 * w, 0.0)
+                    max(we_x[edge_indices.ind_y_down[idx]] - 0.5 * w, 0.0)
             end
         elseif dir == CartesianIndex(1, 0)
             we_y[idx] = reservoir_outlet[v] ? 0.0 : max(we_y[idx] - w, 0.0)
         elseif dir == CartesianIndex(0, 1)
             we_x[idx] = reservoir_outlet[v] ? 0.0 : max(we_x[idx] - w, 0.0)
         elseif dir == CartesianIndex(-1, 0)
-            if edge_indices.idx_left[idx] <= n
-                we_y[edge_indices.idx_left[idx]] =
+            if edge_indices.ind_x_down[idx] <= n
+                we_y[edge_indices.ind_x_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_y[edge_indices.idx_left[idx]] - w, 0.0)
+                    max(we_y[edge_indices.ind_x_down[idx]] - w, 0.0)
             end
         elseif dir == CartesianIndex(0, -1)
-            if edge_indices.idx_down[idx] <= n
-                we_x[edge_indices.idx_down[idx]] =
+            if edge_indices.ind_y_down[idx] <= n
+                we_x[edge_indices.ind_y_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_x[edge_indices.idx_down[idx]] - w, 0.0)
+                    max(we_x[edge_indices.ind_y_down[idx]] - w, 0.0)
             end
         elseif dir == CartesianIndex(1, -1)
             we_y[idx] = max(we_y[idx] - 0.5 * w, 0.0)
-            if edge_indices.idx_down[idx] <= n
-                we_x[edge_indices.idx_down[idx]] =
+            if edge_indices.ind_y_down[idx] <= n
+                we_x[edge_indices.ind_y_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_x[edge_indices.idx_down[idx]] - 0.5 * w, 0.0)
+                    max(we_x[edge_indices.ind_y_down[idx]] - 0.5 * w, 0.0)
             end
         elseif dir == CartesianIndex(-1, 1)
-            if edge_indices.idx_left[idx] <= n
-                we_y[edge_indices.idx_left[idx]] =
+            if edge_indices.ind_x_down[idx] <= n
+                we_y[edge_indices.ind_x_down[idx]] =
                     reservoir_outlet[v] ? 0.0 :
-                    max(we_y[edge_indices.idx_left[idx]] - 0.5 * w, 0.0)
+                    max(we_y[edge_indices.ind_x_down[idx]] - 0.5 * w, 0.0)
             end
             we_x[idx] = reservoir_outlet[v] ? 0.0 : max(we_x[idx] - 0.5 * w, 0.0)
         end
diff --git a/Wflow/test/reservoir.jl b/Wflow/test/reservoir.jl
index 0a06daf44..0297bc9c6 100644
--- a/Wflow/test/reservoir.jl
+++ b/Wflow/test/reservoir.jl
@@ -149,7 +149,7 @@ end
         ),
     )
 
-    river_flow_vars = Wflow.FlowVariables(; n = 2, q = [0.04, 0.04])
+    river_flow_vars = Wflow.RiverFlowVariables(; n = 2, q = [0.04, 0.04])
 
     graph = DiGraph(2)
     add_edge!(graph, 1, 2)
@@ -198,8 +198,11 @@ end
         ),
     )
 
-    river_flow_vars =
-        Wflow.FlowVariables(; n = 2, q = [0.00012002923701686638, 0.21747539140212965])
+    n_river = 2
+    river_flow_vars = Wflow.RiverFlowVariables(;
+        n = n_river,
+        q = [0.00012002923701686638, 0.21747539140212965],
+    )
 
     graph = DiGraph(2)
     add_edge!(graph, 1, 2)
diff --git a/Wflow/test/routing_process.jl b/Wflow/test/routing_process.jl
index bea022e48..01bfcea4e 100644
--- a/Wflow/test/routing_process.jl
+++ b/Wflow/test/routing_process.jl
@@ -78,7 +78,7 @@ end
     # a topological sort is used for visiting nodes in order from upstream to downstream
     toposort = topological_sort_by_dfs(graph)
     sink = toposort[end]
-    @test ldd[sink] == 5  # the most downstream node must be a sink
+    @test ldd[sink] == Wflow.LDD_PIT  # the most downstream node must be a sink
 
     # calculate parameters of kinematic wave
     q = 0.000001
@@ -287,383 +287,6 @@ end
     @test flux == flux_
 end
 
-@testitem "unit: kinwave_river_update!" begin
-    using Graphs: DiGraph, add_edge!
-    n = 2
-    model = Wflow.KinWaveRiverFlowModel(;
-        timestepping = Wflow.TimeStepping(),
-        boundary_conditions = Wflow.RiverFlowBC(;
-            n,
-            external_inflow = [-0.1],
-            reservoir = Wflow.ReservoirModel(;
-                boundary_conditions = Wflow.ReservoirBC(; n, external_inflow = [0.02]),
-                parameters = Wflow.ReservoirParameters(;
-                    id = [1],
-                    storage_curve_type = [Wflow.ReservoirProfileType.linear],
-                    outflow_curve_type = [Wflow.ReservoirOutflowType.simple],
-                    area = [2500.0],
-                ),
-                variables = Wflow.ReservoirVariables(;
-                    waterlevel = [3.0],
-                    storage = [7500.0],
-                ),
-            ),
-        ),
-        parameters = Wflow.RiverFlowParameters(;
-            flow = Wflow.ManningFlowParameters(;
-                slope = [0.01],
-                mannings_n = [0.03],
-                alpha_pow = 0.4,
-                alpha_term = [0.5],
-                alpha = [5.0],
-            ),
-            bankfull_depth = [10.0],
-        ),
-        variables = Wflow.FlowVariables(; n, q = [0.2]),
-        allocation = Wflow.NoAllocationRiverModel(n),
-    )
-    graph = DiGraph(2)
-    add_edge!(graph, 1, 2)
-    domain = Wflow.DomainRiver(;
-        network = Wflow.NetworkRiver(;
-            graph,
-            order_of_subdomains = [[1]],
-            order_subdomain = [[1]],
-            subdomain_indices = [[1]],
-            upstream_nodes = [[]],
-            reservoir_indices = [1],
-        ),
-        parameters = Wflow.RiverParameters(; flow_width = [30.0], flow_length = [800.0]),
-    )
-    dt = 1200.0
-
-    Wflow.kinwave_river_update!(model, domain, dt)
-
-    @test model.variables.q ≈ [0.1598124775930105]
-    @test model.variables.h[1] ≈ 0.055464507410878765
-    @test model.variables.storage[1] ≈ 1331.1481778610903
-    @test model.variables.q_cumulative[1] ≈ 191.7749731116126
-end
-
-@testitem "unit: local_inertial_river_update!" begin
-    dt = 86400.0
-    n = 2
-    river_flow_model = Wflow.LocalInertialRiverFlowModel(;
-        timestepping = Wflow.TimeStepping(),
-        boundary_conditions = Wflow.RiverFlowBC(;
-            n,
-            external_inflow = [-1.0, -1.0],
-            inwater = [100.0, 100.0],
-            reservoir = Wflow.ReservoirModel(;
-                boundary_conditions = Wflow.ReservoirBC(;
-                    n,
-                    external_inflow = [-1.0],
-                    inflow_overland = [3000.0],
-                    inflow_subsurface = [5000.0],
-                    precipitation = [2.3148148148148148e-8],
-                    evaporation = [1.1574074074074074e-8],
-                ),
-                parameters = Wflow.ReservoirParameters(;
-                    id = [1, 2],
-                    storage_curve_type = [Wflow.ReservoirProfileType.linear],
-                    outflow_curve_type = [Wflow.ReservoirOutflowType.simple],
-                    area = [5.0e6, 4.0e6],
-                    maximum_release = [10.0, 10.0],
-                    demand = [1.5, 1.5],
-                    target_minimum_fraction = [0.3, 0.3],
-                    target_full_fraction = [0.1, 0.1],
-                    maximum_storage = [Inf, Inf],
-                    threshold = [0.0, 0.0],
-                ),
-                variables = Wflow.ReservoirVariables(;
-                    waterlevel = [1.0, 1.0],
-                    storage = [2.5e5, 2.5e5],
-                    outflow = [1.8, 1.8],
-                ),
-            ),
-        ),
-        parameters = Wflow.LocalInertialRiverFlowParameters(;
-            n = n,
-            ne = 2,
-            active_n = [2],
-            active_e = [1],
-            froude_limit = true,
-            h_thresh = 0.0,
-            zb = [0.0, 0.0],
-            zb_max = [0.5, 0.5],
-            bankfull_storage = [1e3, 1e3],
-            bankfull_depth = [0.5, 0.5],
-            mannings_n_sq = [1e-3, 1e-3],
-            mannings_n = [1e-2, 1e-2],
-            flow_length_at_edge = [1000.0, 1000.0],
-            flow_width_at_edge = [100.0, 100.0],
-        ),
-        variables = Wflow.LocalInertialRiverFlowVariables(;
-            n_cells = n,
-            n_edges = 2,
-            h = [1.0, 2.0],
-            q = [0.0, 1e-4],
-        ),
-        floodplain = Wflow.FloodPlainModel(;
-            parameters = Wflow.FloodPlainParameters(;
-                profile = Wflow.FloodPlainProfile(;
-                    depth = [10.0, 10.0],
-                    storage = [1e5 1e6; 1e5 1e6],
-                    width = [100.0 100.0; 100.0 100.0],
-                    flow_area = [1e3 1e3; 1e3 1e3],
-                    wetted_perimeter = [400.0 400.0; 400.0 400.0],
-                ),
-                mannings_n = [0.04],
-                mannings_n_sq = [1.2e-4],
-                zb_max = [1.0],
-            ),
-            variables = Wflow.FloodPlainVariables(; n, n_edges = 1, h = [0.1, 0.2]),
-        ),
-        allocation = Wflow.AllocationRiverModel(; n),
-    )
-    domain = Wflow.Domain(;
-        river = Wflow.DomainRiver(;
-            network = Wflow.NetworkRiver(;
-                nodes_at_edge = Wflow.NodesAtEdge(; src = [1], dst = [2]),
-                edges_at_node = Wflow.EdgesAtNode(; src = [[1], [1]], dst = [[1], [1]]),
-            ),
-            parameters = Wflow.RiverParameters(;
-                flow_width = [10.0, 10.0],
-                flow_length = [100.0, 100.0],
-            ),
-        ),
-        reservoir = Wflow.DomainReservoir(;
-            network = Wflow.NetworkReservoir(; river_indices = [1]),
-        ),
-    )
-    dt = 1000.0
-
-    Wflow.update_river_channel_flow!(river_flow_model, domain.river, dt)
-
-    @test river_flow_model.variables.zs_src[1] ≈ 1.0
-    @test river_flow_model.variables.zs_dst[1] ≈ 2.0
-    @test river_flow_model.variables.zs_max[1] ≈ 2.0
-    @test river_flow_model.variables.water_depth_at_edge[1] ≈ 1.5
-    @test river_flow_model.variables.flow_area[1] ≈ 150.0
-    @test river_flow_model.variables.hydraulic_radius[1] ≈ 1.4563106796116505
-    @test river_flow_model.variables.q[1] ≈ -575.3037784510024
-    @test river_flow_model.variables.q_cumulative[1] ≈ -575303.7784510023
-
-    Wflow.update_floodplain_flow!(river_flow_model, domain.river, dt)
-
-    @test river_flow_model.floodplain.variables.water_depth_at_edge[1] ≈ 1.0
-    @test river_flow_model.floodplain.variables.hf_index[1] == 1
-    @test river_flow_model.floodplain.variables.flow_area[1] ≈ 90.0
-    @test river_flow_model.floodplain.variables.hydraulic_radius[1] ≈ 0.2356020942408377
-    @test river_flow_model.floodplain.variables.q[1] ≈ -281.84014086002725
-    @test river_flow_model.floodplain.variables.q_cumulative[1] ≈ -281840.1408600272
-
-    Wflow.update_bc_reservoir_model!(river_flow_model, domain, dt)
-
-    @test river_flow_model.boundary_conditions.reservoir.variables.storage[1] ≈
-          7.391913765967477e6
-    @test river_flow_model.boundary_conditions.reservoir.variables.waterlevel[1] ≈
-          2.428382753193495
-    @test river_flow_model.boundary_conditions.reservoir.variables.outflow[1] ≈
-          0.00018509186397934759
-    @test river_flow_model.boundary_conditions.reservoir.boundary_conditions.inflow_cumulative[1] ≈
-          7.141856080688971e6
-    @test river_flow_model.boundary_conditions.reservoir.variables.outflow_cumulative[1] ≈
-          0.1850918639793476
-    @test river_flow_model.boundary_conditions.reservoir.variables.actevap_cumulative[1] ≈
-          1.1574074074074073e-5
-    @test river_flow_model.variables.q[1] ≈ 0.00018509186397934759
-    @test river_flow_model.variables.q_cumulative[1] ≈ -575303.5933591384
-
-    Wflow.update_water_depth_and_storage!(river_flow_model, domain.river, dt)
-
-    @test river_flow_model.variables.storage[2] ≈ 99000.0
-    @test river_flow_model.variables.h[2] ≈ 99.0
-
-    Wflow.update_water_depth_and_storage!(
-        river_flow_model.floodplain,
-        river_flow_model,
-        domain.river,
-        dt,
-    )
-
-    @test river_flow_model.variables.h ≈ [1.0, -79.7]
-    @test river_flow_model.variables.storage ≈ [0.0, -79700.0]
-    @test river_flow_model.floodplain.variables.storage ≈ [0.0, 178700.0]
-end
-
-@testitem "unit: update_directional_flow!" begin
-    n = 3
-    overland_flow_model = Wflow.LocalInertialOverlandFlowModel(;
-        timestepping = Wflow.TimeStepping(),
-        boundary_conditions = Wflow.LocalInertialOverlandFlowBC(; n),
-        parameters = Wflow.LocalInertialOverlandFlowParameters(;
-            n,
-            ywidth = fill(900.0, n),
-            xwidth = [250.0, 300.0, 450.0],
-            zx_max = [750.0, 900.0, 800.0],
-            theta = 1.0,
-            h_thresh = 1e-3,
-            zy_max = [800.0, 800.0, 800.0],
-            mannings_n_sq = [0.06, 0.06, 0.06],
-            z = [800.0, 800.0, 800.0],
-            froude_limit = true,
-        ),
-        variables = Wflow.LocalInertialOverlandFlowVariables(;
-            n,
-            qx0 = [1e-3, 2e-3, 3e-3],
-            h = [0.03, 0.02, 0.05],
-        ),
-    )
-    domain = Wflow.Domain(;
-        land = Wflow.DomainLand(;
-            network = Wflow.NetworkLand(;
-                edge_indices = Wflow.EdgeConnectivity(;
-                    n = 1,
-                    idx_right = [2],
-                    idx_left = [3],
-                ),
-            ),
-            parameters = Wflow.LandParameters(;
-                x_length = fill(600.0, n),
-                y_length = fill(900.0, n),
-            ),
-        ),
-    )
-    i = 1
-    dt = 60.0
-    is_x_direction = true
-
-    Wflow.update_directional_flow!(overland_flow_model, domain, i, dt, is_x_direction)
-    @test overland_flow_model.variables.qx_cumulative[1] ≈ 26493.90166029366
-end
-
-@testitem "unit: local_inertial_update_water_depth!" begin
-    n = 2
-    overland_flow_model = Wflow.LocalInertialOverlandFlowModel(;
-        timestepping = Wflow.TimeStepping(),
-        variables = Wflow.LocalInertialOverlandFlowVariables(;
-            n,
-            qx = [0.1, 0.3],
-            qy = [0.25, 0.15],
-            storage = [1000.0, 1250.0],
-        ),
-        boundary_conditions = Wflow.LocalInertialOverlandFlowBC(; n, runoff = [0.2, 0.3]),
-        parameters = Wflow.LocalInertialOverlandFlowParameters(;
-            n,
-            xwidth = [600.0],
-            ywidth = [900.0],
-            theta = 1.0,
-            h_thresh = 1e-3,
-            zx_max = [800.0],
-            zy_max = [800.0],
-            mannings_n_sq = [0.06],
-            z = [700.0],
-            froude_limit = true,
-        ),
-    )
-    river_flow_model = Wflow.LocalInertialRiverFlowModel(;
-        timestepping = Wflow.TimeStepping(),
-        boundary_conditions = Wflow.RiverFlowBC(;
-            n,
-            external_inflow = [-0.2, -0.1],
-            reservoir = nothing,
-        ),
-        parameters = Wflow.LocalInertialRiverFlowParameters(;
-            n,
-            ne = 2,
-            active_n = [1, 1],
-            active_e = [1, 1],
-            froude_limit = true,
-            h_thresh = 1e-3,
-            zb = [500.0, 500.0],
-            zb_max = [400.0, 400.0],
-            bankfull_storage = [500.0, 500.0],
-            bankfull_depth = [1.0, 1.0],
-            mannings_n_sq = [9.0e-3, 9.0e-3],
-            mannings_n = [0.03, 0.03],
-            flow_length_at_edge = [800.0, 800.0],
-            flow_width_at_edge = [30.0, 30.0],
-        ),
-        variables = Wflow.LocalInertialRiverFlowVariables(;
-            n_cells = n,
-            n_edges = 2,
-            q = [0.03, 0.04],
-            h = [1.0, 1.0],
-            storage = [30.0e3, 25e3],
-        ),
-        floodplain = nothing,
-        allocation = Wflow.NoAllocationRiverModel(1),
-    )
-    domain = Wflow.Domain(;
-        land = Wflow.DomainLand(;
-            network = Wflow.NetworkLand(;
-                river_indices = [1],
-                edge_indices = Wflow.EdgeConnectivity(;
-                    n,
-                    idx_left = [2, 1],
-                    idx_down = [2, 1],
-                ),
-            ),
-            parameters = Wflow.LandParameters(;
-                x_length = [600.0, 600.0],
-                y_length = [900.0, 900.0],
-            ),
-        ),
-        river = Wflow.DomainRiver(;
-            network = Wflow.NetworkRiver(;
-                edges_at_node = Wflow.EdgesAtNode(; src = [[1]], dst = [[2]]),
-            ),
-            parameters = Wflow.RiverParameters(;
-                flow_width = [30.0],
-                flow_length = [800.0],
-            ),
-        ),
-    )
-
-    dt = 1000.0
-
-    # update_river_cell_storage_and_depth!
-    @test Wflow.compute_river_storage_change(
-        overland_flow_model,
-        river_flow_model,
-        domain,
-        1,
-        dt,
-    ) ≈ 290.0
-    @test Wflow.compute_external_inflow(river_flow_model, overland_flow_model, 1, 1, dt) |>
-          collect ≈ [-0.2, 0.2]
-
-    river_h_expected = 1.0014614814814815
-    land_h_expected = 0.001461481481481508
-    river_storage_expected = 24035.07555555556
-    @test Wflow.compute_water_depths(1289.2, 1, 1, river_flow_model, domain) |> collect ≈
-          [river_h_expected, land_h_expected, river_storage_expected]
-
-    Wflow.update_river_cell_storage_and_depth!(
-        overland_flow_model,
-        river_flow_model,
-        domain,
-        1,
-        dt,
-    )
-    @test river_flow_model.variables.h[1] ≈ 1.0010925925925926
-    @test overland_flow_model.variables.h[1] ≈ 0.001092592592592645
-    @test river_flow_model.variables.storage[1] ≈ 24026.222222222223
-
-    # update_land_cell_storage_and_depth!
-    @test Wflow.compute_land_storage_change(
-        overland_flow_model,
-        domain.land.network,
-        2,
-        dt,
-    ) ≈ 200.0
-
-    Wflow.update_land_cell_storage_and_depth!(overland_flow_model, domain.land, 2, dt)
-    @test overland_flow_model.variables.storage[2] ≈ 1450.0
-    @test overland_flow_model.variables.h[2] ≈ 0.002685185185185185
-end
-
 @testitem "unit: kinwave_river_update!" begin
     # Test river kinematic wave routing on a 2-node graph (1 → 2).
     # Node 1 has a simple reservoir (outflow_curve_type = simple) and a negative external inflow
@@ -672,10 +295,14 @@ end
     # evaporation for the reservoir.
     using Graphs: DiGraph, add_edge!
 
-    dt = 86400.0
     n = 2
+    bankfull_depth = [1.0, 1.0]
+    flow_width = [94.73094177246094, 94.73094177246094]
+    flow_length = [1059.8125, 951.96875]
+    bankfull_storage = @. flow_length * flow_width * bankfull_depth
 
-    river_flow_model = Wflow.KinWaveRiverFlowModel(;
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.KinematicWave(),
         timestepping = Wflow.TimeStepping(; stable_timesteps = zeros(n)),
         boundary_conditions = Wflow.RiverFlowBC(;
             n,
@@ -715,10 +342,15 @@ end
                 alpha_term = [0.41038937516728013, 0.4113871693777356],
                 alpha = [2.544585458995107, 2.5507721996678145],
             ),
-            bankfull_depth = [1.0, 1.0],
+            bankfull_depth,
+            bankfull_storage,
+        ),
+        variables = Wflow.RiverFlowVariables(;
+            n,
+            q = [0.5499295110293246, 3.0005238507869465],
         ),
-        variables = Wflow.FlowVariables(; n, q = [0.5499295110293246, 3.0005238507869465]),
         allocation = Wflow.NoAllocationRiverModel(n),
+        floodplain = nothing,
     )
     graph = DiGraph(2)
     add_edge!(graph, 1, 2)
@@ -731,10 +363,7 @@ end
             upstream_nodes = [[], [1]],
             reservoir_indices = [1, 0],
         ),
-        parameters = Wflow.RiverParameters(;
-            flow_width = [94.73094177246094, 94.73094177246094],
-            flow_length = [1059.8125, 951.96875],
-        ),
+        parameters = Wflow.RiverParameters(; flow_width, flow_length),
     )
     dt = Wflow.stable_timestep(river_flow_model, domain.parameters.flow_length, 0.05)
 
@@ -763,7 +392,8 @@ end
     #   4. update_water_depth_and_storage!
     model_dt = 86400.0
     n = 3
-    river_flow_model = Wflow.LocalInertialRiverFlowModel(;
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping = Wflow.TimeStepping(),
         boundary_conditions = Wflow.RiverFlowBC(;
             n,
@@ -795,23 +425,27 @@ end
                 ),
             ),
         ),
-        parameters = Wflow.LocalInertialRiverFlowParameters(;
+        parameters = Wflow.RiverFlowStaggeredParameters(;
             n,
-            ne = 2,
+            n_edges = 2,
             active_n = [1, 3],
             active_e = [1],
             froude_limit = true,
             h_thresh = 0.001,
             zb = [315.1000061035156, 314.3999938964844, 278.8000183105469],
-            zb_max = [315.1000061035156, 314.3999938964844],
+            zb_at_edge = [315.1000061035156, 314.3999938964844],
             bankfull_storage = [57155.32165002823, 100397.03622722626, 90180.89622545242],
             bankfull_depth = [1.0, 1.0, 1.0],
-            mannings_n_sq = [0.0008999999597668652, 0.0008999999597668652],
-            mannings_n = [0.029999999329447746, 0.029999999329447746, 0.029999999329447746],
+            mannings_n_sq_at_edge = [0.0008999999597668652, 0.0008999999597668652],
+            mannings_n_at_edge = [
+                0.029999999329447746,
+                0.029999999329447746,
+                0.029999999329447746,
+            ],
             flow_length_at_edge = [831.578125, 1005.890625],
             flow_width_at_edge = [94.73094177246094, 94.73094177246094],
         ),
-        variables = Wflow.LocalInertialRiverFlowVariables(;
+        variables = Wflow.RiverFlowStaggeredVariables(;
             n_cells = n,
             n_edges = 2,
             h = [0.04484241735240722, 0.0, 0.07939389691400389],
@@ -847,10 +481,8 @@ end
     @test dt ≈ 909.829412320351
     @test river_flow_model.variables.zs_src ≈ [315.14484852086804, 0.0]
     @test river_flow_model.variables.zs_dst ≈ [314.3999938964844, 0.0]
-    @test river_flow_model.variables.zs_max ≈ [315.14484852086804, 0.0]
+    @test river_flow_model.variables.zs_at_edge ≈ [315.14484852086804, 0.0]
     @test river_flow_model.variables.water_depth_at_edge ≈ [0.04484241735241312, 0.0]
-    @test river_flow_model.variables.flow_area ≈ [4.247964427147839, 0.0]
-    @test river_flow_model.variables.hydraulic_radius ≈ [0.04480000374545946, 0.0]
     @test river_flow_model.variables.q ≈ [0.534558444239785, 0.0]
     @test river_flow_model.variables.q_cumulative ≈ [486.35699517356477, 0.0]
 
@@ -895,7 +527,7 @@ end
     @test river_flow_model.variables.h ≈ [0.03653704705843743, 0.0, 0.10966599406601261]
 end
 
-@testitem "unit: local inertial river flow with flow_area floodplain" begin
+@testitem "unit: local inertial river flow with floodplain" begin
     # Test local inertial river routing with a floodplain on a 3-node graph (1 → 2 → 3).
     # Each sub-step of the local inertial update is called and verified individually:
     #   1. update_river_channel_flow!  — edge discharge, water surface elevations
@@ -903,7 +535,8 @@ end
     #   3. update_bc_reservoir_model!  — no reservoir
     #   4. update_water_depth_and_storage!
     n = 3
-    river_flow_model = Wflow.LocalInertialRiverFlowModel(;
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7),
         boundary_conditions = Wflow.RiverFlowBC(;
             n,
@@ -911,23 +544,24 @@ end
             inwater = [0.001214603164946035, 0.0069799723162954465, 0.00022016439899281862],
             reservoir = nothing,
         ),
-        parameters = Wflow.LocalInertialRiverFlowParameters(;
+        parameters = Wflow.RiverFlowStaggeredParameters(;
             n,
-            ne = 2,
+            n_edges = 2,
             active_n = [1, 2, 3],
             active_e = [1, 2],
             froude_limit = true,
             h_thresh = 0.001,
             zb = [165.10784077644348, 165.10784077644348, 165.10784077644348],
-            zb_max = [165.10784077644348, 165.10784077644348],
+            zb_at_edge = [165.10784077644348, 165.10784077644348],
             bankfull_storage = [107921.00118967971, 200292.17449130033, 69688.46493603183],
             bankfull_depth = [1.592156171798706, 1.592156171798706, 1.592156171798706],
-            mannings_n_sq = [0.0008999999597668652, 0.0008999999597668652],
             mannings_n = [0.029999999329447746, 0.029999999329447746, 0.029999999329447746],
+            mannings_n_sq_at_edge = [0.0008999999597668652, 0.0008999999597668652],
+            mannings_n_at_edge = [0.029999999329447746, 0.029999999329447746],
             flow_length_at_edge = [648.828125, 568.34375],
             flow_width_at_edge = [149.17837524414062, 149.17837524414062],
         ),
-        variables = Wflow.LocalInertialRiverFlowVariables(;
+        variables = Wflow.RiverFlowStaggeredVariables(;
             n_cells = n,
             n_edges = 2,
             h = [1.8817912224982847, 1.8197068314233162, 1.7619620034455687],
@@ -935,7 +569,8 @@ end
             q = [137.1750567107639, 133.74797838974442],
         ),
         floodplain = Wflow.FloodPlainModel(;
-            parameters = Wflow.FloodPlainParameters(;
+            routing_method = Wflow.LocalInertial(),
+            parameters = Wflow.FloodPlainStaggeredParameters(;
                 profile = Wflow.FloodPlainProfile(;
                     depth = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5],
                     storage = [
@@ -947,35 +582,36 @@ end
                         747343.0 1.87239e6 814606.0
                     ],
                     width = [
-                        149.178 149.178 149.178
-                        341.577 666.779 758.764
-                        492.562 778.794 988.691
-                        693.057 881.476 1118.98
-                        789.594 988.161 1237.77
-                        972.752 1125.52 1448.54
+                        149.17837524414062 149.17837524414062 149.17837524414062;
+                        341.5768913342503 666.778728923476 758.7636596016615;
+                        492.562310866575 778.7935519733185 988.6905953775695;
+                        693.0574965612104 881.4757828423199 1118.9809351368624;
+                        789.5944979367263 988.1614230127849 1237.7718606880392;
+                        972.7515818431912 1125.5153603853994 1448.5444669293854;
                     ],
                     flow_area = [
-                        0.0 0.0 0.0
-                        170.788 333.389 379.382
-                        417.07 722.786 873.727
-                        763.598 1163.52 1433.22
-                        1158.4 1657.6 2052.1
-                        1644.77 2220.36 2776.38
+                        0.0 0.0 0.0;
+                        170.78844566712516 333.389364461738 379.38182980083076;
+                        417.0696011004127 722.7861404483972 873.7271274896154;
+                        763.5983493810179 1163.5240318695571 1433.2175950580468;
+                        1158.395598349381 1657.6047433759495 2052.103525402066;
+                        1644.7713892709767 2220.3624235686493 2776.375758866759;
                     ],
                     wetted_perimeter = [
-                        192.399 517.6 609.585
-                        193.399 518.6 610.585
-                        345.384 631.615 841.512
-                        546.879 735.297 972.803
-                        644.416 842.983 1092.59
-                        828.573 981.337 1304.37
+                        192.3985160901097 517.6003536793354 609.5852843575209;
+                        193.3985160901097 518.6003536793354 610.5852843575209;
+                        345.38393562243436 631.6151767291778 841.5122201334289;
+                        546.8791213170698 735.2974075981792 972.8025598927218;
+                        644.4161226925856 842.9830477686443 1092.5934854438985;
+                        828.5732065990505 981.3369851412588 1304.3660916852448;
                     ],
                 ),
                 mannings_n = [0.072, 0.072, 0.072],
-                mannings_n_sq = [0.005184, 0.005184],
-                zb_max = [166.6999969482422, 166.6999969482422],
+                mannings_n_at_edge = [0.072, 0.072],
+                mannings_n_sq_at_edge = [0.005184, 0.005184],
+                zb_at_edge = [166.6999969482422, 166.6999969482422],
             ),
-            variables = Wflow.FloodPlainVariables(;
+            variables = Wflow.FloodPlainStaggeredVariables(;
                 n,
                 n_edges = 2,
                 q = [3.306660222819796, 6.14041420989245],
@@ -1008,12 +644,9 @@ end
     @test dt ≈ 49.40931052556788
     @test river_flow_model.variables.zs_src ≈ [166.98963199894177, 166.92754760786679]
     @test river_flow_model.variables.zs_dst ≈ [166.92754760786679, 166.86980277988906]
-    @test river_flow_model.variables.zs_max ≈ [166.98963199894177, 166.92754760786679]
+    @test river_flow_model.variables.zs_at_edge ≈ [166.98963199894177, 166.92754760786679]
     @test river_flow_model.variables.water_depth_at_edge ≈
           [1.8817912224982933, 1.8197068314233036]
-    @test river_flow_model.variables.flow_area ≈ [280.7225571209805, 271.4609085323917]
-    @test river_flow_model.variables.hydraulic_radius ≈
-          [1.8354842671202385, 1.7763698223484754]
     @test river_flow_model.variables.q ≈ [137.1827776559179, 133.7538757670657]
     @test river_flow_model.variables.q_cumulative ≈ [6778.106459961183, 6608.686781773178]
 
@@ -1021,14 +654,10 @@ end
 
     @test river_flow_model.floodplain.variables.water_depth_at_edge ≈
           [0.2896350506995873, 0.22755065962459753]
-    @test river_flow_model.floodplain.variables.hf_index == [1, 2]
-    @test river_flow_model.floodplain.variables.flow_area ≈
-          [55.72538543569419, 117.7803635852996]
-    @test river_flow_model.floodplain.variables.hydraulic_radius ≈
-          [0.28876507912737354, 0.22735103521877678]
-    @test river_flow_model.floodplain.variables.q ≈ [3.3074651032168534, 6.14213116671929]
+    @test river_flow_model.floodplain.variables.q ≈ [3.3074672215578524, 6.1421232455587536]
+
     @test river_flow_model.floodplain.variables.q_cumulative ≈
-          [163.41957033732095, 303.47846610520196]
+          [163.41967500308917, 303.4780747261213]
 
     Wflow.update_bc_reservoir_model!(
         river_flow_model.boundary_conditions.reservoir,
@@ -1051,20 +680,318 @@ end
     )
 
     @test river_flow_model.variables.storage ≈
-          [124521.72968438723, 228924.52867849727, 78479.82315286293]
+          [124521.73491986065, 228924.54042985922, 78479.82701991481]
     @test river_flow_model.variables.h ≈
-          [1.8370663564507992, 1.8197595694254676, 1.7930102910029269]
+          [1.8370664336896243, 1.8197596628390198, 1.793010379352563]
     @test river_flow_model.floodplain.variables.storage ≈
-          [21410.323896476315, 99344.9914658879, 35924.01098543289]
+          [21410.318556337137, 99344.98021057082, 35924.006727001935]
     @test river_flow_model.floodplain.variables.h ≈
-          [0.2449101846520932, 0.2276033976267614, 0.2008541192042207]
+          [0.2449102618909183, 0.22760349104031377, 0.20085420755385677]
+end
+
+@testitem "unit: river flow on a staggered grid using Manning's equation with floodplain" begin
+    # Test river flow on a staggered grid using Manning's equation with a floodplain on a
+    # 3-node graph (1 → 2 → 3).
+    # Each sub-step of the update of Manning's flow on staggered grid  is called and
+    # verified individually:
+    #   1. update_river_channel_flow!  — edge discharge, water surface elevations
+    #   2. update_floodplain_flow!     — floodplain discharge and geometry
+    #   3. update_bc_reservoir_model!  — no reservoir
+    #   4. update_water_depth_and_storage!
+    n = 3
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.ManningStaggered(),
+        timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7),
+        boundary_conditions = Wflow.RiverFlowBC(;
+            n,
+            external_inflow = zeros(n),
+            inwater = [0.5469041025637776, 1.2789538388011605, 0.1723191462506401],
+            reservoir = nothing,
+        ),
+        parameters = Wflow.RiverFlowStaggeredParameters(;
+            n,
+            n_edges = 2,
+            active_n = [1, 2, 3],
+            active_e = [1, 2],
+            h_thresh = 0.001,
+            zb = [415.6000061035156, 415.6000061035156, 415.6000061035156],
+            zb_at_edge = [415.6000061035156, 415.6000061035156],
+            bankfull_storage = [36208.125, 30283.125, 25717.5],
+            bankfull_depth = [1.0, 1.0, 1.0],
+            mannings_n = [0.029999999329447746, 0.029999999329447746, 0.029999999329447746],
+            mannings_n_at_edge = [0.029999999329447746, 0.029999999329447746],
+            flow_length_at_edge = [1108.1875, 933.34375],
+            flow_width_at_edge = [30.0, 30.0],
+            slope_at_edge = [1.0e-5, 1.0e-5],
+        ),
+        variables = Wflow.RiverFlowStaggeredVariables(;
+            n_cells = n,
+            n_edges = 2,
+            h = [1.0029472866230495, 2.432593354336019, 0.40220945912422745],
+            storage = [36314.840722458204, 73666.52862352695, 10343.82176502732],
+            q = [12.586919865724491, 12.586919865724491],
+            water_depth_at_edge = [2.4325367293470777, 2.4325367293470777],
+        ),
+        floodplain = Wflow.FloodPlainModel(;
+            routing_method = Wflow.ManningStaggered(),
+            parameters = Wflow.FloodPlainStaggeredParameters(;
+                profile = Wflow.FloodPlainProfile(;
+                    depth = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5],
+                    storage = [
+                        0.0 0.0 0.0
+                        37399.0 36242.0 34514.0;
+                        79402.0 73635.0 69028.0;
+                        150174.0 111029.0 103542.0;
+                        233029.0 155325.0 138056.0;
+                        341777.0 212852.0 172570.0;
+                    ],
+                    width = [
+                        30.0 30.0 30.0;
+                        61.97338304593237 71.8063277815615 80.52260134149898;
+                        69.60260991144943 74.08680577054052 80.52260134149898;
+                        117.27533530112372 74.08878707200793 80.52260134149898;
+                        137.2979131065196 87.76372980001238 80.52260134149898;
+                        180.20485733519754 113.9783295152003 80.52260134149898;
+                    ],
+                    flow_area = [
+                        0.0 0.0 0.0;
+                        30.986691522966186 35.90316389078075 40.26130067074949;
+                        65.7879964786909 72.94656677605101 80.52260134149898;
+                        124.42566412925277 109.99096031205498 120.78390201224848;
+                        193.07462068251255 153.87282521206117 161.04520268299797;
+                        283.1770493501113 210.86198996966132 201.30650335374747;
+                    ],
+                    wetted_perimeter = [
+                        31.973383045932373 41.80632778156151 50.522601341498984;
+                        32.97338304593237 42.80632778156151 51.522601341498984;
+                        41.602609911449434 46.08680577054052 52.522601341498984;
+                        90.27533530112372 47.08878707200793 53.522601341498984;
+                        111.29791310651959 61.763729800012385 54.522601341498984;
+                        155.20485733519754 88.9783295152003 55.522601341498984;
+                    ],
+                ),
+                mannings_n = [0.072, 0.072, 0.072],
+                mannings_n_at_edge = [0.072, 0.072],
+                zb_at_edge = [416.6000061035156, 416.6000061035156],
+                slope_at_edge = [1.0e-5, 1.0e-5],
+            ),
+            variables = Wflow.FloodPlainStaggeredVariables(;
+                n,
+                n_edges = 2,
+                q = [2.0759155671111853, 3.279029180032925],
+                water_depth_at_edge = [1.4325367293470777, 1.4325367293470777],
+                h = [0.0029472866230494782, 1.432593354336019, 0.0],
+                storage = [113.73542237265065, 62604.388160555245, 0.0],
+            ),
+        ),
+        allocation = Wflow.AllocationRiverModel(; n),
+    )
+    domain = Wflow.Domain(;
+        river = Wflow.DomainRiver(;
+            network = Wflow.NetworkRiver(;
+                nodes_at_edge = Wflow.NodesAtEdge(; src = [1, 2], dst = [2, 3]),
+                edges_at_node = Wflow.EdgesAtNode(;
+                    src = [[], [1], [2]],
+                    dst = [[1], [2], []],
+                ),
+            ),
+            parameters = Wflow.RiverParameters(;
+                flow_width = [30.0, 30.0, 30.0],
+                flow_length = [1206.9375, 1009.4375, 857.25],
+            ),
+        ),
+        reservoir = Wflow.DomainReservoir(; network = Wflow.NetworkReservoir()),
+    )
+    dt = Wflow.stable_timestep(river_flow_model, domain.river.parameters.flow_length)
+
+    Wflow.update_river_channel_flow!(river_flow_model, domain.river, dt)
+
+    @test dt ≈ 2272.7523105527266
+    @test river_flow_model.variables.zs_src ≈ [416.6029533901387, 418.03259945785163]
+    @test river_flow_model.variables.zs_dst ≈ [418.03259945785163, 416.00221556263983]
+    @test river_flow_model.variables.zs_at_edge ≈ [416.6029533901387, 418.03259945785163]
+    @test river_flow_model.variables.water_depth_at_edge ≈
+          [1.0029472866230495, 2.432593354336019]
+    @test river_flow_model.variables.q ≈ [3.0436243193280137, 12.587380945649572]
+    @test river_flow_model.variables.q_cumulative ≈ [6917.404204207213, 28607.999128032432]
+
+    Wflow.update_floodplain_flow!(river_flow_model, domain.river, dt)
+
+    @test river_flow_model.floodplain.variables.water_depth_at_edge ≈
+          [0.0029472866230494782, 1.432593354336019]
+    @test river_flow_model.floodplain.variables.q ≈ [8.50693948416018e-5, 3.279243909752732]
+
+    @test river_flow_model.floodplain.variables.q_cumulative ≈
+          [0.1933416636835727, 7452.909172756478]
+
+    Wflow.update_bc_reservoir_model!(
+        river_flow_model.boundary_conditions.reservoir,
+        river_flow_model,
+        domain,
+        dt,
+    )
+
+    Wflow.update_water_depth_and_storage!(river_flow_model, domain.river, dt)
+
+    @test river_flow_model.variables.storage ≈
+          [30640.414081003582, 54882.67899192734, 39343.459630853366]
+    @test river_flow_model.variables.h ≈
+          [0.8462303441838974, 1.8123188736937599, 1.529832201063609]
+
+    Wflow.update_water_depth_and_storage!(
+        river_flow_model.floodplain,
+        river_flow_model,
+        domain.river,
+        dt,
+    )
+
+    @test river_flow_model.variables.storage ≈
+          [30753.95616171255, 63042.829748341144, 33570.77415623857]
+    @test river_flow_model.variables.h ≈
+          [0.849366161923948, 2.0817808514920815, 1.3053669352090433]
+    @test river_flow_model.floodplain.variables.storage ≈
+          [0.0, 46991.52157304865, 13225.594647371276]
+    @test river_flow_model.floodplain.variables.h ≈
+          [0.0, 1.0817808514920815, 0.3053669352090434]
+end
+
+@testitem "unit: kinematic river flow including 1D floodplain schematization" begin
+    # Test river kinematic wave routing with floodplain on a 2-node graph (1 → 2).
+    # Each sub-step of the update of kinematic wave routing with floodplain is called and
+    # verified individually:
+    #   1. river_channel_floodplain_exchange!
+    #   2. kinwave_river_update!
+    #   3. update_floodplain_model!
+    using Graphs: DiGraph, add_edge!
+
+    n = 2
+    flow_length = [750.953125, 851.8125]
+    flow_width = [229.91920471191406, 229.91920471191406]
+
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.KinematicWave(),
+        timestepping = Wflow.TimeStepping(; stable_timesteps = zeros(n)),
+        boundary_conditions = Wflow.RiverFlowBC(;
+            n,
+            inwater = [0.02953203136486251, 0.0001849569866329713],
+            reservoir = nothing,
+        ),
+        parameters = Wflow.RiverFlowParameters(;
+            flow = Wflow.ManningFlowParameters(;
+                slope = [1.0e-5, 1.0e-5],
+                mannings_n = [0.03, 0.03],
+                alpha_pow = 0.4,
+                alpha_term = [3.8572052304976667, 3.8572052304976667],
+                alpha = [34.0789466790827, 34.0789466790827],
+            ),
+            bankfull_depth = [2.1051321029663086, 2.1051321029663086],
+            bankfull_storage = [363469.04651181493, 412286.02275520907],
+        ),
+        variables = Wflow.RiverFlowVariables(;
+            n,
+            q = [296.52948301601174, 192.8313119108856],
+            qlat = [3.9326064945614956e-5, 2.1713344971219756e-7],
+            h = [4.509741437854894, 3.4835130995322534],
+            storage = [778645.3962305915, 682239.2566234164],
+        ),
+        allocation = Wflow.NoAllocationRiverModel(n),
+        floodplain = Wflow.FloodPlainModel(;
+            routing_method = Wflow.Manning(),
+            parameters = Wflow.FloodPlainParameters(;
+                profile = Wflow.FloodPlainProfile(;
+                    depth = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5],
+
+                    storage = [
+                        0.0 0.0;
+                        86329.2726379633 97924.02628183365;
+                        172659.2726379633 207762.02628183365;
+                        258989.2726379633 386716.02628183365;
+                        369518.2726379633 605009.0262818336;
+                        724648.2726379633 869843.0262818336;
+                    ],
+                    width = [
+                        229.91920471191406 229.91920471191406;
+                        229.91920471191406 229.91920471191406;
+                        229.9211418821914 257.89243524836746;
+                        229.9211418821914 420.1722796977034;
+                        294.369904912507 512.5376770122533;
+                        945.8113647239966 621.8128989654414;
+                    ],
+                    flow_area = [
+                        0.0 0.0;
+                        114.95960235595703 114.95960235595703;
+                        229.9201732970527 243.90581998014076;
+                        344.8807442381484 453.99195982899244;
+                        492.06569669440194 710.260798335119;
+                        964.9713790564002 1021.1672478178398;
+                    ],
+                    wetted_perimeter = [
+                        0.0 0.0;
+                        1.0 1.0;
+                        2.00193717027733 29.9732305364534;
+                        3.00193717027733 193.25307498578934;
+                        68.45070020059296 286.61847230033925;
+                        720.8921600120825 396.8936942535273;
+                    ],
+                ),
+                mannings_n = [0.072, 0.072],
+                slope = [1.0e-5, 1.0e-5],
+            ),
+            variables = Wflow.FloodPlainVariables(;
+                n,
+                q = [1.2861909826521447, 1.9846650910027395],
+                h = [0.0, 0.0],
+                storage = [0.0, 0.0],
+            ),
+        ),
+    )
+    graph = DiGraph(2)
+    add_edge!(graph, 1, 2)
+    domain = Wflow.DomainRiver(;
+        network = Wflow.NetworkRiver(;
+            graph,
+            order = [1, 2],
+            order_of_subdomains = [[1]],
+            order_subdomain = [[1, 2]],
+            subdomain_indices = [[1, 2]],
+            upstream_nodes = [[], [1]],
+            reservoir_indices = [0, 0],
+        ),
+        parameters = Wflow.RiverParameters(; flow_width, flow_length),
+    )
+    dt = Wflow.stable_timestep(river_flow_model, domain.parameters.flow_length, 0.05)
+
+    Wflow.river_channel_floodplain_exchange!(river_flow_model, domain.parameters, dt)
+
+    @test dt ≈ 1602.881460805217
+    @test river_flow_model.variables.h ≈ [4.509741437854894, 3.4835130995322534]
+    @test river_flow_model.variables.storage ≈ [778645.3962305915, 682239.2566234164]
+    @test river_flow_model.floodplain.variables.h ≈ [2.0642836103410205, 1.1737631111525133]
+    @test river_flow_model.floodplain.variables.storage ≈
+          [58760.1445203583, 40074.01437791623]
+
+    Wflow.kinwave_river_update!(river_flow_model, domain, dt)
+
+    @test river_flow_model.variables.h ≈ [2.872002930358695, 3.1138609375883397]
+    @test river_flow_model.variables.storage ≈ [495875.84798393055, 609843.6005807515]
+    @test river_flow_model.variables.q ≈ [139.7837242183345, 159.9486203252721]
+    @test river_flow_model.variables.q_cumulative ≈ [224056.7400718776, 256378.67820075116]
+
+    Wflow.update_floodplain_model!(river_flow_model, domain, dt)
+
+    @test river_flow_model.floodplain.variables.storage ≈
+          [52745.30941770246, 41649.967280840756]
+    @test river_flow_model.floodplain.variables.q ≈ [3.752513987924126, 2.7693140810933157]
+    @test river_flow_model.floodplain.variables.q_cumulative ≈
+          [6014.835102655834, 4438.882199731312]
 end
 
 @testitem "unit: update_directional_flow!" begin
     # Test local inertial overland flow routing in x-direction at edge 2 (i = 2) using 3
     # nodes.
     n = 3
-    overland_flow_model = Wflow.LocalInertialOverlandFlowModel(;
+    overland_flow_model = Wflow.OverlandFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7),
         boundary_conditions = Wflow.LocalInertialOverlandFlowBC(;
             n,
@@ -1072,13 +999,17 @@ end
         ),
         parameters = Wflow.LocalInertialOverlandFlowParameters(;
             n,
-            ywidth = [926.6857061478484, 869.7426481339323, 812.7995901200163],
-            xwidth = [],
-            zx_max = [257.3280029296875, 232.67100524902344, 232.67100524902344],
+            ywidth_at_edge = [926.6857061478484, 869.7426481339323, 812.7995901200163],
+            xwidth_at_edge = [],
+            zx_max_at_edge = [257.3280029296875, 232.67100524902344, 232.67100524902344],
             theta = 1.0,
             h_thresh = 1e-3,
-            zy_max = [],
-            mannings_n_sq = [0.24167056670421605, 0.2883451232664811, 0.3928782408368683],
+            zy_max_at_edge = [],
+            mannings_n_sq_at_edge = [
+                0.24167056670421605,
+                0.2883451232664811,
+                0.3928782408368683,
+            ],
             z = [257.3280029296875, 227.5050048828125, 232.67100524902344],
             froude_limit = true,
         ),
@@ -1095,8 +1026,8 @@ end
             network = Wflow.NetworkLand(;
                 edge_indices = Wflow.EdgeConnectivity(;
                     n = 3,
-                    idx_right = [2, 3, 4],
-                    idx_left = [4, 1, 2],
+                    ind_x_up = [2, 3, 4],
+                    ind_x_down = [4, 1, 2],
                 ),
             ),
             parameters = Wflow.LandParameters(;
@@ -1121,7 +1052,8 @@ end
     # overland and river flow (subgrid channel) routing, with 2 land cells and 1 river cell:
     # | land | land | river |.
     n_land = 3
-    overland_flow_model = Wflow.LocalInertialOverlandFlowModel(;
+    overland_flow_model = Wflow.OverlandFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7),
         variables = Wflow.LocalInertialOverlandFlowVariables(;
             n = n_land,
@@ -1136,42 +1068,39 @@ end
         ),
         parameters = Wflow.LocalInertialOverlandFlowParameters(;
             n = n_land,
-            xwidth = [],
-            ywidth = [],
+            xwidth_at_edge = [],
+            ywidth_at_edge = [],
             theta = 1.0,
             h_thresh = 1e-3,
-            zx_max = [],
-            zy_max = [],
-            mannings_n_sq = [],
+            zx_max_at_edge = [],
+            zy_max_at_edge = [],
+            mannings_n_sq_at_edge = [],
             z = [],
             froude_limit = true,
         ),
     )
     n_river = 1
-    river_flow_model = Wflow.LocalInertialRiverFlowModel(;
+    river_flow_model = Wflow.RiverFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7),
-        boundary_conditions = Wflow.RiverFlowBC(;
+        boundary_conditions = Wflow.RiverFlowBC(; n = n_river, reservoir = nothing),
+        parameters = Wflow.RiverFlowStaggeredParameters(;
             n = n_river,
-            reservoir = nothing,
-            external_inflow = zeros(n_river),
-        ),
-        parameters = Wflow.LocalInertialRiverFlowParameters(;
-            n = n_river,
-            ne = 2,
+            n_edges = 2,
             active_n = [1],
             active_e = [1, 2],
             froude_limit = true,
             h_thresh = 1e-3,
             zb = [],
-            zb_max = [],
+            zb_at_edge = [],
             bankfull_storage = [171137.5821314017],
             bankfull_depth = [1.3683528900146484],
-            mannings_n_sq = [],
-            mannings_n = [],
+            mannings_n_sq_at_edge = [],
+            mannings_n_at_edge = [],
             flow_length_at_edge = [],
             flow_width_at_edge = [],
         ),
-        variables = Wflow.LocalInertialRiverFlowVariables(;
+        variables = Wflow.RiverFlowStaggeredVariables(;
             n_cells = n_river,
             n_edges = 2,
             q = [56.685647296907476, 53.70963118023338],
@@ -1189,8 +1118,8 @@ end
                 river_indices = [0, 0, 1],
                 edge_indices = Wflow.EdgeConnectivity(;
                     n = 3,
-                    idx_left = [4, 1, 2],
-                    idx_down = [4, 4, 4],
+                    ind_x_down = [4, 1, 2],
+                    ind_y_down = [4, 4, 4],
                 ),
             ),
             parameters = Wflow.LandParameters(;
@@ -1212,7 +1141,7 @@ end
     dt = Wflow.stable_timestep(river_flow_model, domain.river.parameters.flow_length)
     @test dt ≈ 201.394687315008
 
-    # Test functions called by function update_river_cell_storage_and_depth!
+    # Test functions called by function update_river_and_land_storage_and_depth!
     storage_change = Wflow.compute_river_storage_change(
         overland_flow_model,
         river_flow_model,
@@ -1231,8 +1160,8 @@ end
     @test Wflow.compute_water_depths(total_storage, 1, 3, river_flow_model, domain) |>
           collect ≈ [river_h_expected, land_h_expected, river_storage_expected]
 
-    # Test update_river_cell_storage_and_depth!
-    Wflow.update_river_cell_storage_and_depth!(
+    # Test update_river_and_land_storage_and_depth!
+    Wflow.update_river_and_land_storage_and_depth!(
         overland_flow_model,
         river_flow_model,
         domain,
@@ -1243,7 +1172,7 @@ end
     @test overland_flow_model.variables.h[3] ≈ land_h_expected
     @test river_flow_model.variables.storage[1] ≈ river_storage_expected
 
-    # Test function called by function update_land_cell_storage_and_depth!
+    # Test function called by function update_land_storage_and_depth!
     @test Wflow.compute_land_storage_change(
         overland_flow_model,
         domain.land.network,
@@ -1251,8 +1180,8 @@ end
         dt,
     ) ≈ 880.1704436259577
 
-    # Test update_land_cell_storage_and_depth!
-    Wflow.update_land_cell_storage_and_depth!(overland_flow_model, domain.land, 2, dt)
+    # Test update_land_storage_and_depth!
+    Wflow.update_land_storage_and_depth!(overland_flow_model, domain.land, 2, dt)
     @test overland_flow_model.variables.storage[2] ≈ 784038.1273051423
     @test overland_flow_model.variables.h[2] ≈ 1.3770166561681556
 end
@@ -1298,23 +1227,24 @@ end
 
     # for each edge the src and dst node is required
     nodes_at_edge = Wflow.adjacent_nodes_at_edge(graph)
-    _ne = ne(graph)
+    n_edges = ne(graph)
 
     # determine z, width, length and manning's n at edges
-    zb_max = fill(0.0, _ne)
-    width_at_edge = fill(0.0, _ne)
-    length_at_edge = fill(0.0, _ne)
-    mannings_n_sq = fill(0.0, _ne)
-    for i in 1:_ne
-        zb_max[i] = max(zb[nodes_at_edge.src[i]], zb[nodes_at_edge.dst[i]])
-        width_at_edge[i] = min(width[nodes_at_edge.dst[i]], width[nodes_at_edge.src[i]])
-        length_at_edge[i] = 0.5 * (dl[nodes_at_edge.dst[i]] + dl[nodes_at_edge.src[i]])
-        mannings_n =
+    zb_at_edge = fill(0.0, n_edges)
+    flow_width_at_edge = fill(0.0, n_edges)
+    flow_length_at_edge = fill(0.0, n_edges)
+    mannings_n_sq_at_edge = fill(0.0, n_edges)
+    for i in 1:n_edges
+        zb_at_edge[i] = max(zb[nodes_at_edge.src[i]], zb[nodes_at_edge.dst[i]])
+        flow_width_at_edge[i] =
+            min(width[nodes_at_edge.dst[i]], width[nodes_at_edge.src[i]])
+        flow_length_at_edge[i] = 0.5 * (dl[nodes_at_edge.dst[i]] + dl[nodes_at_edge.src[i]])
+        mannings_n_at_edge =
             (
                 n_river[nodes_at_edge.dst[i]] * dl[nodes_at_edge.dst[i]] +
                 n_river[nodes_at_edge.src[i]] * dl[nodes_at_edge.src[i]]
             ) / (dl[nodes_at_edge.dst[i]] + dl[nodes_at_edge.src[i]])
-        mannings_n_sq[i] = mannings_n * mannings_n
+        mannings_n_sq_at_edge[i] = mannings_n_at_edge * mannings_n_at_edge
     end
 
     river_network = Wflow.NetworkRiver(;
@@ -1338,30 +1268,29 @@ end
     push!(h_init, h_a[n])
 
     timestepping = Wflow.TimeStepping(; alpha_coefficient = 0.7)
-    parameters = Wflow.LocalInertialRiverFlowParameters(;
+    parameters = Wflow.RiverFlowStaggeredParameters(;
         n,
-        ne = _ne,
+        n_edges,
         active_n = collect(1:(n - 1)),
-        active_e = collect(1:_ne),
+        active_e = collect(1:n_edges),
         h_thresh,
-        zb_max,
-        mannings_n_sq,
-        mannings_n = n_river,
-        flow_width_at_edge = width_at_edge,
-        flow_length_at_edge = length_at_edge,
+        zb_at_edge,
+        mannings_n_sq_at_edge,
+        flow_width_at_edge,
+        flow_length_at_edge,
         bankfull_storage = fill(Wflow.MISSING_VALUE, n),
         bankfull_depth = fill(Wflow.MISSING_VALUE, n),
         zb,
         froude_limit,
     )
 
-    variables =
-        Wflow.LocalInertialRiverFlowVariables(; n_cells = n, n_edges = _ne, h = h_init)
+    variables = Wflow.RiverFlowStaggeredVariables(; n_cells = n, n_edges, h = h_init)
 
     boundary_conditions =
         Wflow.RiverFlowBC(; n, external_inflow = zeros(n), reservoir = nothing)
 
-    sw_river = Wflow.LocalInertialRiverFlowModel(;
+    sw_river = Wflow.RiverFlowModel(;
+        routing_method = Wflow.LocalInertial(),
         timestepping,
         boundary_conditions,
         parameters,
@@ -1377,7 +1306,7 @@ end
         sw_river.boundary_conditions.inwater[1] = 20.0
         h0 = mean(sw_river.variables.h)
         dt = Wflow.stable_timestep(sw_river, flow_length)
-        Wflow.local_inertial_river_update!(sw_river, domain, dt, true)
+        Wflow.staggered_scheme_river_update!(sw_river, domain, dt, true)
         d = abs(h0 - mean(sw_river.variables.h))
         if d <= epsilon
             break
@@ -1397,7 +1326,7 @@ end
     A = 0.04970535373017763
     R = 0.0011733219820725962
     length = 533.453125
-    mannings_n_sq = 0.0008999999597668652
+    mannings_n_sq_at_edge = 0.0008999999597668652
     froude_limit = true
     dt = 89.29563868855615
 
@@ -1409,7 +1338,7 @@ end
         A,
         R,
         length,
-        mannings_n_sq,
+        mannings_n_sq_at_edge,
         froude_limit,
         dt,
     ) ≈ 0.005331926324969742
@@ -1424,7 +1353,7 @@ end
     water_depth_at_edge = 0.00310727852479431
     width = 620.6649135473787
     length = 926.602742473319
-    mannings_n_sq = 0.1773345894316103
+    mannings_n_sq_at_edge = 0.1773345894316103
     froude_limit = true
     dt = 49.774905820268735
 
@@ -1438,7 +1367,7 @@ end
         water_depth_at_edge,
         width,
         length,
-        mannings_n_sq,
+        mannings_n_sq_at_edge,
         froude_limit,
         dt,
     ) ≈ 0.00017992597962222483
@@ -1453,7 +1382,7 @@ end
     A = 0.04970535373017763
     R = 0.0011733219820725962
     length = 533.453125
-    mannings_n_sq = 0.0008999999597668652
+    mannings_n_sq_at_edge = 0.0008999999597668652
     froude_limit = true
     dt = 89.29563868855615
 
@@ -1465,7 +1394,7 @@ end
         A,
         R,
         length,
-        mannings_n_sq,
+        mannings_n_sq_at_edge,
         froude_limit,
         dt,
     ) ≈ 0.005331926324969742
@@ -1480,7 +1409,7 @@ end
     water_depth_at_edge = 0.00310727852479431
     width = 620.6649135473787
     length = 926.602742473319
-    mannings_n_sq = 0.1773345894316103
+    mannings_n_sq_at_edge = 0.1773345894316103
     froude_limit = true
     dt = 49.774905820268735
 
@@ -1494,7 +1423,7 @@ end
         water_depth_at_edge,
         width,
         length,
-        mannings_n_sq,
+        mannings_n_sq_at_edge,
         froude_limit,
         dt,
     ) ≈ 0.00017992597962222483
diff --git a/Wflow/test/run_sbm.jl b/Wflow/test/run_sbm.jl
index 894575398..3f537beb7 100644
--- a/Wflow/test/run_sbm.jl
+++ b/Wflow/test/run_sbm.jl
@@ -357,7 +357,7 @@
             Dict(:to_river_average => 7.373928060045205e-5),
         )
         @test test_means(
-            overland_flow.variables.flow,
+            overland_flow.variables,
             Dict(
                 :qin => 0.0008814607677874327,
                 :storage => 16.837627069604082,
@@ -557,7 +557,7 @@
             Dict(:to_river_average => 0.0006441014921670179),
         )
         @test test_means(
-            overland_flow.variables.flow,
+            overland_flow.variables,
             Dict(
                 :qin => 0.006930522031750806,
                 :storage => 88.26124686793727,
@@ -882,7 +882,7 @@ end
 end
 
 @testitem "Local-inertial option for river flow including 1D floodplain schematization" begin
-    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-local-inertial_config.toml")
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-staggered-scheme_config.toml")
     config = Wflow.Config(tomlpath)
     config.dir_output = mktempdir()
     model = Wflow.Model(config)
@@ -896,8 +896,9 @@ end
 
     @testset "river flow (local inertial) floodplain schematization" begin
         # floodplain geometry checks (index 3)
-        @test profile.storage[:, 3] ≈ [0.0, 8641.0, 19011.0, 31685.0, 51848.0, 80653.0]
-        @test profile.width[:, 3] ≈ [
+        idx = 3
+        @test profile.storage[:, idx] ≈ [0.0, 8641.0, 19011.0, 31685.0, 51848.0, 80653.0]
+        @test profile.width[:, idx] ≈ [
             30.0,
             99.28617594254938,
             119.15260323159785,
@@ -905,7 +906,7 @@ end
             231.6754039497307,
             330.9730700179533,
         ]
-        @test profile.wetted_perimeter[:, 3] ≈ [
+        @test profile.wetted_perimeter[:, idx] ≈ [
             69.28617594254938,
             70.28617594254938,
             91.15260323159785,
@@ -913,7 +914,7 @@ end
             205.6754039497307,
             305.9730700179533,
         ]
-        @test profile.flow_area[:, 3] ≈ [
+        @test profile.flow_area[:, idx] ≈ [
             0.0,
             49.64308797127469,
             109.21938958707361,
@@ -921,116 +922,63 @@ end
             297.8700179533214,
             463.35655296229805,
         ]
-        @test dh .* profile.width[2:end, 3] * flow_length[3] ≈ Δv
-        @test profile.flow_area[:, 3] * flow_length[3] ≈ profile.storage[:, 3]
+        @test dh .* profile.width[2:end, idx] * flow_length[idx] ≈ Δv
+        @test profile.flow_area[:, 3] * flow_length[idx] ≈ profile.storage[:, idx]
         # flood depth from flood storage (8000.0)
         flood_vol = 8000.0
-        river_flow.variables.storage[3] =
-            flood_vol + river_flow.parameters.bankfull_storage[3]
-        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, 3])
+        river_flow.variables.storage[idx] =
+            flood_vol + river_flow.parameters.bankfull_storage[idx]
+        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, idx])
         @test (i1, i2) == (1, 2)
-        flood_depth = Wflow.flood_depth(profile, flood_vol, flow_length[3], 3)
+        flood_depth = Wflow.compute_flood_depth(profile, flood_vol, flow_length[idx], idx)
         @test flood_depth ≈ 0.46290938548779076
-        @test (flood_depth - profile.depth[i1]) * profile.width[i2, 3] * flow_length[3] +
-              profile.storage[i1, 3] ≈ flood_vol
+        @test (flood_depth - profile.depth[i1]) *
+              profile.width[i2, idx] *
+              flow_length[idx] + profile.storage[i1, idx] ≈ flood_vol
         # flood depth from flood storage (12000.0)
         flood_vol = 12000.0
-        river_flow.variables.storage[3] =
-            flood_vol + river_flow.parameters.bankfull_storage[3]
-        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, 3])
-        @test (i1, i2) == (2, 3)
-        flood_depth = Wflow.flood_depth(profile, flood_vol, flow_length[3], 3)
+        river_flow.variables.storage[idx] =
+            flood_vol + river_flow.parameters.bankfull_storage[idx]
+        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, idx])
+        @test (i1, i2) == (2, idx)
+        flood_depth = Wflow.compute_flood_depth(profile, flood_vol, flow_length[idx], idx)
         @test flood_depth ≈ 0.6619575699132112
-        @test (flood_depth - profile.depth[i1]) * profile.width[i2, 3] * flow_length[3] +
-              profile.storage[i1, 3] ≈ flood_vol
+        @test (flood_depth - profile.depth[i1]) *
+              profile.width[i2, idx] *
+              flow_length[idx] + profile.storage[i1, idx] ≈ flood_vol
         # test extrapolation of segment
         flood_vol = 95000.0
-        river_flow.variables.storage[3] =
-            flood_vol + river_flow.parameters.bankfull_storage[3]
-        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, 3])
+        river_flow.variables.storage[idx] =
+            flood_vol + river_flow.parameters.bankfull_storage[idx]
+        i1, i2 = Wflow.interpolation_indices(flood_vol, profile.storage[:, idx])
         @test (i1, i2) == (6, 6)
-        flood_depth = Wflow.flood_depth(profile, flood_vol, flow_length[3], 3)
+        flood_depth = Wflow.compute_flood_depth(profile, flood_vol, flow_length[idx], idx)
         @test flood_depth ≈ 2.749036625585836
-        @test (flood_depth - profile.depth[i1]) * profile.width[i2, 3] * flow_length[3] +
-              profile.storage[i1, 3] ≈ flood_vol
-        river_flow.variables.storage[3] = 0.0 # reset storage
+        @test (flood_depth - profile.depth[i1]) *
+              profile.width[i2, idx] *
+              flow_length[idx] + profile.storage[i1, idx] ≈ flood_vol
+        river_flow.variables.storage[idx] = 0.0 # reset storage
         # flow area and wetted perimeter based on hf
         h = 0.5
         i1, i2 = Wflow.interpolation_indices(h, profile.depth)
-        @test Wflow.flow_area(
-            profile.width[i2, 3],
-            profile.flow_area[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 49.64308797127469
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 70.28617594254938
+        @test Wflow.compute_flood_flow_area(profile, h, idx, i1, i2) ≈ 49.64308797127469
+        @test Wflow.compute_wetted_perimeter(profile, h, idx, i1) ≈ 70.28617594254938
         h = 1.5
         i1, i2 = Wflow.interpolation_indices(h, profile.depth)
-        @test Wflow.flow_area(
-            profile.width[i2, 3],
-            profile.flow_area[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 182.032315978456
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 118.62585278276481
+        @test Wflow.compute_flood_flow_area(profile, h, idx, i1, i2) ≈ 182.032315978456
+        @test Wflow.compute_wetted_perimeter(profile, h, idx, i1) ≈ 118.62585278276481
         h = 1.7
         i1, i2 = Wflow.interpolation_indices(h, profile.depth)
-        @test Wflow.flow_area(
-            profile.width[i2, 3],
-            profile.flow_area[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 228.36739676840216
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 119.02585278276482
+        @test Wflow.compute_flood_flow_area(profile, h, idx, i1, i2) ≈ 228.36739676840216
+        @test Wflow.compute_wetted_perimeter(profile, h, idx, i1) ≈ 119.02585278276482
         h = 3.2
         i1, i2 = Wflow.interpolation_indices(h, profile.depth)
-        @test Wflow.flow_area(
-            profile.width[i2, 3],
-            profile.flow_area[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 695.0377019748654
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 307.3730700179533
+        @test Wflow.compute_flood_flow_area(profile, h, idx, i1, i2) ≈ 695.0377019748654
+        @test Wflow.compute_wetted_perimeter(profile, h, idx, i1) ≈ 307.3730700179533
         h = 4.0
         i1, i2 = Wflow.interpolation_indices(h, profile.depth)
-        @test Wflow.flow_area(
-            profile.width[i2, 3],
-            profile.flow_area[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 959.816157989228
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 3],
-            profile.depth[i1],
-            h,
-        ) ≈ 308.9730700179533
-        @test Wflow.flow_area(
-            profile.width[i2, 4],
-            profile.flow_area[i1, 4],
-            profile.depth[i1],
-            h,
-        ) ≈ 407.6395313908081
-        @test Wflow.wetted_perimeter(
-            profile.wetted_perimeter[i1, 4],
-            profile.depth[i1],
-            h,
-        ) ≈ 90.11775307900271
+        @test Wflow.compute_flood_flow_area(profile, h, idx, i1, i2) ≈ 959.816157989228
+        @test Wflow.compute_wetted_perimeter(profile, h, idx, i1) ≈ 308.9730700179533
     end
 
     Wflow.run_timestep!(model)
@@ -1219,6 +1167,57 @@ end
     end
 end
 
+@testitem "River flow including 1D floodplain schematization using Manning's equation on a staggered grid" begin
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-staggered-scheme_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.model.river_routing = "manning_staggered"
+    config.dir_output = mktempdir()
+    model = Wflow.Model(config)
+    (; river_flow) = model.routing
+
+    Wflow.run_timestep!(model)
+    Wflow.run_timestep!(model)
+
+    (; q_average, h) = river_flow.variables
+    @test sum(q_average) ≈ 2156.7740284069446
+    @test q_average[1622] ≈ 0.000216893766604954
+    @test q_average[43] ≈ 10.362921270030188
+    @test q_average[501] ≈ 0.021301823784704774
+    @test q_average[5808] ≈ 0.0049185166802011605
+    @test h[1622] ≈ 0.001787783986962345
+    @test h[43] ≈ 1.3195684133779013
+    @test h[501] ≈ 0.005358321292557224
+    @test h[5808] ≈ 0.005888613968265871
+    (; q_average, h) = river_flow.floodplain.variables
+    @test maximum(q_average) ≈ 0.8692896468290046
+    @test maximum(h) ≈ 1.1111707025189723
+end
+
+@testitem "Kinematic river flow including 1D floodplain schematization" begin
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-kw_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.dir_output = mktempdir()
+    model = Wflow.Model(config)
+    (; river_flow) = model.routing
+
+    Wflow.run_timestep!(model)
+    Wflow.run_timestep!(model)
+
+    (; q_average, h) = river_flow.variables
+    @test sum(q_average) ≈ 3664.1750799952056
+    @test q_average[1622] ≈ 0.000750262072108146
+    @test q_average[43] ≈ 11.647181706571397
+    @test q_average[501] ≈ 0.1624521220430912
+    @test q_average[5808] ≈ 0.04397694950619904
+    @test h[1622] ≈ 0.0009790317470641876
+    @test h[43] ≈ 0.16302929690891788
+    @test h[501] ≈ 0.1542630086179196
+    @test h[5808] ≈ 0.005801865652235309
+    (; q_average, h) = river_flow.floodplain.variables
+    @test maximum(q_average) ≈ 0.606783057904105
+    @test maximum(h) ≈ 1.5799312217287604
+end
+
 @testitem "run wflow sbm" begin
     tomlpath = joinpath(@__DIR__, "sbm_config.toml")
     config = Wflow.Config(tomlpath)
@@ -1293,7 +1292,7 @@ end
 end
 
 @testitem "water balance river local inertial routing with floodplain" begin
-    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-local-inertial_config.toml")
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-staggered-scheme_config.toml")
     config = Wflow.Config(tomlpath)
     config.dir_output = mktempdir()
     config.model.water_mass_balance__flag = true
@@ -1318,6 +1317,53 @@ end
     Wflow.close_files(model; delete_output = false)
 end
 
+@testitem "water balance river flow with floodplain using Manning's equation on a staggered grid" begin
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-staggered-scheme_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.model.river_routing = "manning_staggered"
+    config.model.water_mass_balance__flag = true
+    config.dir_output = mktempdir()
+    model = Wflow.Model(config)
+    (; river_water_balance) = model.mass_balance.routing
+    Wflow.run_timestep!(model)
+    @testset "water balance first timestep" begin
+        @test all(e -> abs(e) < 1e-9, river_water_balance.error)
+        @test all(re -> abs(re) < 1e-9, river_water_balance.relative_error)
+    end
+    Wflow.run_timestep!(model)
+    @testset "water balance second timestep" begin
+        @test all(e -> abs(e) < 0.00012, river_water_balance.error)
+        @test all(re -> abs(re) < 0.33, river_water_balance.relative_error)
+        inds = findall(x -> x > 2.0e-3, model.routing.river_flow.variables.q_average)
+        @test all(e -> abs(e) < 1e-9, river_water_balance.error[inds])
+        @test all(re -> abs(re) < 1e-9, river_water_balance.relative_error[inds])
+    end
+    Wflow.close_files(model; delete_output = false)
+end
+
+@testitem "water balance river kinematic wave routing with floodplain" begin
+    tomlpath = joinpath(@__DIR__, "sbm_river-floodplain-kw_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.dir_output = mktempdir()
+    config.model.water_mass_balance__flag = true
+    model = Wflow.Model(config)
+    (; river_water_balance) = model.mass_balance.routing
+    Wflow.run_timestep!(model)
+    @testset "water balance first timestep" begin
+        @test all(e -> abs(e) < 1e-9, river_water_balance.error)
+        @test all(re -> abs(re) < 1e-9, river_water_balance.relative_error)
+    end
+    Wflow.run_timestep!(model)
+    @testset "water balance second timestep" begin
+        @test all(e -> abs(e) < 3e-5, river_water_balance.error)
+        @test all(re -> abs(re) < 12.2, river_water_balance.relative_error)
+        inds = findall(x -> x > 1e-3, model.routing.river_flow.variables.q_average)
+        @test all(re -> abs(re) < 1e-9, river_water_balance.error[inds])
+        @test all(re -> abs(re) < 1e-9, river_water_balance.relative_error[inds])
+    end
+    Wflow.close_files(model; delete_output = false)
+end
+
 @testitem "water balance river and land local inertial routing" begin
     tomlpath = joinpath(@__DIR__, "sbm_river-land-local-inertial_config.toml")
     config = Wflow.Config(tomlpath)
diff --git a/Wflow/test/sbm_river-floodplain-kw_config.toml b/Wflow/test/sbm_river-floodplain-kw_config.toml
new file mode 100644
index 000000000..688899fdb
--- /dev/null
+++ b/Wflow/test/sbm_river-floodplain-kw_config.toml
@@ -0,0 +1,232 @@
+# This is a TOML configuration file for Wflow.
+# Relative file paths are interpreted as being relative to this TOML file.
+# Wflow documentation https://deltares.github.io/Wflow.jl/dev/
+# TOML documentation: https://github.com/toml-lang/toml
+
+dir_input = "data/input"
+dir_output = "data/output"
+
+[logging]
+loglevel = "info"
+
+[time]
+calendar = "proleptic_gregorian"
+endtime = "2000-02-01T00:00:00"
+starttime = "2000-01-01T00:00:00"
+time_units = "days since 1900-01-01 00:00:00"
+timestepsecs = 86400
+
+[state]
+path_input = "instates-moselle.nc"
+path_output = "outstates-moselle.nc"
+
+# if listed, the variable must be present in the NetCDF or error
+# if not listed, the variable can get a default value if it has one
+
+[state.variables]
+vegetation_canopy_water__depth = "canopystorage"
+
+soil_water_saturated_zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_layer_water_unsaturated_zone__depth = "ustorelayerdepth"
+
+snowpack_dry_snow__leq_depth = "snow"
+snowpack_liquid_water__depth = "snowwater"
+
+river_water__depth = "h_river"
+river_water__instantaneous_volume_flow_rate = "q_river"
+
+floodplain_water__volume = "storage_floodplain"
+reservoir_water_surface__elevation = "waterlevel_reservoir"
+
+subsurface_water__instantaneous_volume_flow_rate = "ssf"
+
+land_surface_water__instantaneous_volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+
+[input]
+path_forcing = "forcing-moselle.nc"
+path_static = "staticmaps-moselle.nc"
+
+# these are not directly part of the model
+basin__local_drain_direction = "wflow_ldd"
+river_location__mask = "wflow_river"
+reservoir_area__count = "wflow_reservoirareas"
+reservoir_location__count = "wflow_reservoirlocs"
+subbasin_location__count = "wflow_subcatch"
+river_gauge__count = "wflow_gauges_grdc"
+
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.static]
+atmosphere_air__snowfall_temperature_threshold = "TT"
+atmosphere_air__snowfall_temperature_interval = "TTI"
+
+land_water_covered__area_fraction = "WaterFrac"
+
+snowpack__melting_temperature_threshold = "TTM"
+snowpack__degree_day_coefficient = "Cfmax"
+snowpack__liquid_water_holding_capacity = "WHC"
+
+soil_layer_water__brooks_corey_exponent = "c"
+soil_surface_water__infiltration_reduction_parameter = "cf_soil"
+soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "f"
+compacted_soil_surface_water__infiltration_capacity = "InfiltCapPath"
+soil_water__residual_volume_fraction = "thetaR"
+soil_water__saturated_volume_fraction = "thetaS"
+soil_water_saturated_zone_bottom__max_leakage_volume_flux = "MaxLeakage"
+compacted_soil__area_fraction = "PathFrac"
+soil_wet_root__sigmoid_function_shape_parameter = "rootdistpar"
+soil__thickness = "SoilThickness"
+
+vegetation_canopy_water__mean_evaporation_to_mean_precipitation_ratio = "EoverR"
+vegetation_canopy__light_extinction_coefficient = "Kext"
+vegetation__specific_leaf_storage = "Sl"
+vegetation_wood_water__storage_capacity = "Swood"
+vegetation_root__depth = "RootingDepth"
+
+river__length = "wflow_riverlength"
+river_water_flow__manning_n_parameter = "N_River"
+river__slope = "RiverSlope"
+river__width = "wflow_riverwidth"
+river_bank_water__depth = "RiverDepth"
+river_bank_water__elevation = "RiverZ"
+floodplain__slope = "RiverSlope"
+floodplain_water__sum_of_volume_per_depth = "floodplain_volume"
+
+land_surface_water_flow__manning_n_parameter = "N"
+land_surface__slope = "Slope"
+land_surface__elevation = "wflow_dem"
+
+reservoir_surface__area = "reservoir_area"
+reservoir_water_demand__required_downstream_volume_flow_rate = "ResDemand"
+reservoir_water_release_below_spillway__max_volume_flow_rate = "ResMaxRelease"
+reservoir_water__max_volume = "ResMaxVolume"
+reservoir_water__target_full_volume_fraction = "ResTargetFullFrac"
+reservoir_water__target_min_volume_fraction = "ResTargetMinFrac"
+reservoir_water_surface__initial_elevation = "waterlevel_reservoir"
+reservoir_water__rating_curve_type_count = "outflowfunc"
+reservoir_water__storage_curve_type_count = "storfunc"
+
+subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio = "KsatHorFrac"
+
+[input.cyclic]
+vegetation__leaf_area_index = "LAI"
+
+[input.static.soil_surface_water__vertical_saturated_hydraulic_conductivity]
+netcdf_variable_name = "KsatVer"
+scale = 1.0
+offset = 0.0
+
+[model]
+kinematic_wave__adaptive_time_step_flag = true
+floodplain_1d__flag = true
+snow_gravitational_transport__flag = true
+cold_start__flag = true
+reservoir__flag = true
+snow__flag = true
+soil_layer__thickness = [100, 300, 800]
+type = "sbm"
+river_streamorder__min_count = 6
+land_streamorder__min_count = 5
+
+[output.netcdf_grid]
+path = "output_moselle.nc"
+
+[output.netcdf_grid.variables]
+soil_water_saturated_zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_layer_water_unsaturated_zone__depth = "ustorelayerdepth"
+snowpack_dry_snow__leq_depth = "snow"
+snowpack_liquid_water__depth = "snowwater"
+river_water__depth = "h_river"
+river_water__volume_flow_rate = "q_av_river"
+reservoir_water__volume = "storage_reservoir"
+subsurface_water__volume_flow_rate = "ssf_av"
+land_surface_water__volume_flow_rate = "q_av_land"
+land_surface_water__depth = "h_land"
+land.interception.variables.canopy_storage = "canopystorage"
+
+[output.netcdf_scalar]
+path = "output_scalar_moselle.nc"
+
+[[output.netcdf_scalar.variable]]
+name = "Q"
+map = "river_gauge__count"
+parameter = "routing.river_flow.variables.q"
+
+[[output.netcdf_scalar.variable]]
+name = "T"
+map = "river_gauge__count"
+parameter = "atmosphere_air__temperature"
+
+[[output.netcdf_scalar.variable]]
+coordinate.x = 6.255
+coordinate.y = 50.012
+name = "temp_coord"
+location = "temp_bycoord"
+parameter = "atmosphere_air__temperature"
+
+[[output.netcdf_scalar.variable]]
+location = "temp_byindex"
+name = "temp_index"
+index.x = 100
+index.y = 264
+parameter = "atmosphere_air__temperature"
+
+[output.csv]
+path = "output_moselle.csv"
+
+[[output.csv.column]]
+header = "Q"
+parameter = "river_water__volume_flow_rate"
+reducer = "maximum"
+
+[[output.csv.column]]
+header = "storage"
+index = 1
+parameter = "reservoir_water__volume"
+
+[[output.csv.column]]
+coordinate.x = 6.255
+coordinate.y = 50.012
+header = "temp_bycoord"
+parameter = "atmosphere_air__temperature"
+
+[[output.csv.column]]
+coordinate.x = 6.255
+coordinate.y = 50.012
+header = "vwc_layer2_bycoord"
+parameter = "land.soil.variables.volumetric_water_content"
+layer = 2
+
+[[output.csv.column]]
+header = "temp_byindex"
+index.x = 100
+index.y = 264
+parameter = "atmosphere_air__temperature"
+
+[[output.csv.column]]
+header = "Q"
+map = "river_gauge__count"
+parameter = "river_water__volume_flow_rate"
+
+[[output.csv.column]]
+header = "recharge"
+map = "subbasin_location__count"
+parameter = "soil_water_saturated_zone_top__net_recharge_volume_flux"
+reducer = "mean"
+
+[API]
+variables = [
+  "river_water__volume_flow_rate",
+  "soil_water_unsaturated_zone__depth",
+  "soil_water__transpiration_volume_flux",
+  "soil_layer_1_water_unsaturated_zone__depth",
+  "soil_layer_2_water_unsaturated_zone__depth",
+  "soil_layer_3_water_unsaturated_zone__depth",
+  "soil_layer_4_water_unsaturated_zone__depth",
+]
diff --git a/Wflow/test/sbm_river-floodplain-local-inertial_config.toml b/Wflow/test/sbm_river-floodplain-staggered-scheme_config.toml
similarity index 100%
rename from Wflow/test/sbm_river-floodplain-local-inertial_config.toml
rename to Wflow/test/sbm_river-floodplain-staggered-scheme_config.toml
diff --git a/Wflow/test/utils.jl b/Wflow/test/utils.jl
index c164a61e6..e2ed79921 100644
--- a/Wflow/test/utils.jl
+++ b/Wflow/test/utils.jl
@@ -137,7 +137,7 @@ end
 
     for file_name in [
         "sbm_gwf_config.toml",
-        "sbm_river-floodplain-local-inertial_config.toml",
+        "sbm_river-floodplain-staggered-scheme_config.toml",
         "sbm_river-land-local-inertial_config.toml",
         "sbm_gwf_piave_demand_config.toml",
         "sediment_config.toml",
@@ -189,6 +189,9 @@ end
                     "soil_exponential_vertical_saturated_hydraulic_conductivity_profile_below_surface__depth",
                     "soil_layer_water__vertical_saturated_hydraulic_conductivity",
                 ])
+            elseif map_name == "routing"
+                # The lens of this standard name is valid but not part of any test model.
+                Set(["floodplain__slope"])
             else
                 Set{String}()
             end
diff --git a/docs/_quarto.yml b/docs/_quarto.yml
index 9c351eeec..ff6991370 100644
--- a/docs/_quarto.yml
+++ b/docs/_quarto.yml
@@ -109,7 +109,7 @@ website:
         - text: Routing concepts
           contents:
           - file: model_docs/routing/kinwave.qmd
-          - file: model_docs/routing/local-inertial.qmd
+          - file: model_docs/routing/staggered-grid.qmd
           - file: model_docs/routing/reservoirs.qmd
           - file: model_docs/routing/gwf.qmd
           - file: model_docs/routing/sediment_flux.qmd
diff --git a/docs/changelog.qmd b/docs/changelog.qmd
index 6f133ae05..f3d0da994 100644
--- a/docs/changelog.qmd
+++ b/docs/changelog.qmd
@@ -67,6 +67,16 @@ subsurface_water__instantaneous_volume_flow_rate = "ssf"
   domain via `input.subbasin_location__count`, or specify specific subcatchments by linking
   `input.subbasin_location__count` to a subcatchment map, and indicate which subbasins
   should be active by providing their IDs in `input.subbasin_active_location__count`.
+- An optional 1D floodplain schematization to kinematic wave routing of river flow. This is
+  the same 1D floodplain schematization used by the local inertial river flow model. River
+  channel-floodplain water exchange is assumed to be instantaneous. Routing is done
+  separately for the river channel and floodplain. Floodplain water is routed down slope
+  using the `accucapacityflux` function with a transport capacity based on Manning's
+  equation.
+- Support the use of Manning's equation (assuming that the topography controls water flow
+  mostly) for river and floodplain routing on a staggered grid. This is an additional
+  routing scheme to the local inertial river routing scheme, both make us of the staggered
+  grid method.
 - Added an optional `wflow_version` field to the config .toml file. This setting is used by
   the wflow model builder `hydromt_wflow` for incremental Wflow.jl version upgrading.
 
diff --git a/docs/developments/model_struct.qmd b/docs/developments/model_struct.qmd
index 6a277c60b..ba976d82d 100644
--- a/docs/developments/model_struct.qmd
+++ b/docs/developments/model_struct.qmd
@@ -8,7 +8,7 @@ configuration and input and output.
 
 ```julia
 """
-    Model{R <: Routing, L <: AbstractLandModel, M <: AbstractMassBalance, T <: AbstractModelType} <: AbstractModel{T}
+    Model{R <: Routing, L <: AbstractLandModel, M <: AbstractMassBalance, W <: Writer, T <: AbstractModelType} <: AbstractModel{T}
 
 Composite type that represents all different aspects of a Wflow Model, such as the network,
 parameters, clock, configuration and input and output.
@@ -17,6 +17,7 @@ struct Model{
     R <: Routing,
     L <: AbstractLandModel,
     M <: AbstractMassBalance,
+    W <: Writer,
     T <: AbstractModelType,
 } <: AbstractModel{T}
     config::Config                  # all configuration options
@@ -26,7 +27,7 @@ struct Model{
     mass_balance::M                 # mass balance error
     clock::Clock                    # to keep track of simulation time
     reader::NCReader                # provides the model with dynamic input
-    writer::Writer                  # writes model output
+    writer::W                       # writes model output
     type::T                         # model type
 end
 ```
diff --git a/docs/model_docs/model_configurations.qmd b/docs/model_docs/model_configurations.qmd
index af4b54d2f..b3cf7ba57 100644
--- a/docs/model_docs/model_configurations.qmd
+++ b/docs/model_docs/model_configurations.qmd
@@ -19,10 +19,14 @@ supported in wflow:
 
     and with the following surface routing schemes:
     - kinematic wave for overland and river flow
-    - kinematic wave for overland flow + local inertial river (+ optional floodplain)
+    - kinematic wave for overland flow + local inertial river
+    - kinematic wave for overland flow + river flow on a staggered grid using Manning's equation
     - local inertial river (1D) and land (2D)
  - wflow\_sediment as post processing of wflow\_sbm output
 
+The river flow routing schemes support an optional 1D floodplain schematization (routing is
+done separately for the river channel and floodplain).
+
 Below, some explanation will be given on how to prepare a basic wflow\_sbm model. Example
 data for other model configurations is provided in the section with [sample
 data](../getting_started/download_example_models.qmd).
@@ -64,16 +68,21 @@ four directions (adjacent cells) is used. This is described in more detail in th
 type = "sbm_gwf"
 ```
 
-### Local inertial river
-By default the model types `sbm` and `sbm_gwf` use the kinematic wave approach for river flow.
-There is also the option to use the local inertial model for river flow with an optional 1D
-floodplain schematization (routing is done separately for the river channel and floodplain), by
-providing the following in the TOML file:
+### River flow on a staggered grid
+By default the model types `sbm` and `sbm_gwf` use the kinematic wave approach for river
+flow, including an optional 1D floodplain schematization. There is also the option to use
+river flow routing on a staggered grid, with the following two routing schemes:
+
+- Local inertial routing
+- Mannings flow
+
+Both support an optional 1D floodplain schematization.
+The river routing scheme can be set as follows in the TOML file:
 
 ```toml
 [model]
-river_routing = "local_inertial"    # optional, default is "kinematic_wave"
-floodplain_1d__flag = true                # optional, default is false
+river_routing = "local_inertial"    # one of ["kinematic_wave", "local_inertial", "manning_staggered"], default is "kinematic_wave"
+floodplain_1d__flag = true          # optional, default is false
 ```
 
 ### Local inertial river (1D) and land (2D)
diff --git a/docs/model_docs/model_settings.qmd b/docs/model_docs/model_settings.qmd
index 54ad547c4..bc5c2fb94 100644
--- a/docs/model_docs/model_settings.qmd
+++ b/docs/model_docs/model_settings.qmd
@@ -55,6 +55,7 @@ The Tables in this subsection show the available models settings for the model t
 | `subsurface_kinematic_wave__alpha_coefficient` | Stability coefficient for internal model time step subsurface flow routing (CFL condition) | - | 1.0 |
 | `drain__flag` | Include drainage boundary condition  | - | `false` |
 | `river_subsurface_exchange_head_based__flag` | Include head-dependent river flux boundary | - | `false` |
+| `floodplain_1d__flag` | Include 1D floodplain schematization for kinematic wave river routing | - | `false` |
 : {.striped .hover}
 
 ### Local inertial routing
@@ -63,13 +64,21 @@ The Tables in this subsection show the available models settings for the model t
 | `river_local_inertial_flow__alpha_coefficient` | Stability coefficient for internal model time step river flow routing (CFL condition) | - | 0.7 |
 | `land_local_inertial_flow__alpha_coefficient` | Stability coefficient for internal model time step overland flow routing (CFL condition) | - | 0.7 |
 | `land_local_inertial_flow__theta_coefficient` | Weighting factor that adjusts the amount of artificial numerical diffusion for local inertial overland flow | - | 1.0 |
-| `river_water_flow_threshold__depth` | Water depth threshold for calculating river flow between cells | m | 0.001 |
-| `land_surface_water_flow_threshold__depth` | Water depth threshold for calculating overland flow between cells | m | 0.001 |
+| `river_water_flow_threshold__depth` | Water depth threshold for calculating river flow between nodes | m | 0.001 |
+| `land_surface_water_flow_threshold__depth` | Water depth threshold for calculating overland flow between nodes | m | 0.001 |
 | `river_water_flow__froude_limit_flag` | Limit river flow to subcritical-critical according to Froude number | - | `true` |
 | `land_surface_water_flow__froude_limit_flag` | Limit overland flow to subcritical-critical according to Froude number | - | `true` |
 | `floodplain_1d__flag` | Include 1D floodplain schematization for local inertial river routing | - | `false` |
 : {.striped .hover}
 
+### River routing using Manning's equation on a staggered grid
+| Setting | Description | Unit | Default |
+| --- | ----- | - | - |
+| `river_staggered_manning_flow__alpha_coefficient` | Stability coefficient for internal model time step river flow routing (CFL condition) | - | 0.7 |
+| `river_water_flow_threshold__depth` | Water depth threshold for calculating river flow between nodes | m | 0.001 |
+| `floodplain_1d__flag` | Include 1D floodplain schematization | - | `false` |
+: {.striped .hover}
+
 ### Groundwater flow
 | Setting | Description | Unit | Default |
 | --- | ----- | - | - |
diff --git a/docs/model_docs/parameters_routing.qmd b/docs/model_docs/parameters_routing.qmd
index cb2d4cd93..d86a2b9a2 100644
--- a/docs/model_docs/parameters_routing.qmd
+++ b/docs/model_docs/parameters_routing.qmd
@@ -70,6 +70,25 @@ generate_table(Wflow.Routing, :kinematic_wave_river_state; relative_widths = [3,
 generate_table(Wflow.Routing, :kinematic_wave_river_output; relative_widths = [3, 3, 1, 1, 2])
 ```
 
+### 1D floodplain flow
+#### Input
+```{julia}
+# | code-fold: true
+generate_table(Wflow.Routing, :kinematic_wave_floodplain_1D_flow_input; relative_widths = [6, 4, 2, 2, 2, 2])
+```
+
+#### States
+```{julia}
+# | code-fold: true
+generate_table(Wflow.Routing, :kinematic_wave_floodplain_1D_flow_state; relative_widths = [3, 3, 1, 1, 2])
+```
+
+#### Output
+```{julia}
+# | code-fold: true
+generate_table(Wflow.Routing, :kinematic_wave_floodplain_1D_flow_output; relative_widths = [3, 3, 1, 1, 2])
+```
+
 ### Overland flow
 #### Input
 ```{julia}
@@ -108,52 +127,55 @@ generate_table(Wflow.Routing, :kinematic_lateral_subsurface_state; relative_widt
 generate_table(Wflow.Routing, :kinematic_lateral_subsurface_output; relative_widths = [3, 3, 1, 1, 2])
 ```
 
-## Local inertial
+## Staggered grid
+River and floodplain flow routing on a staggered grid make use of the local inertial or
+Manning's equation.
+
 ### River flow
 #### Input
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_river_input; relative_widths = [6, 4, 2, 2, 2, 2])
+generate_table(Wflow.Routing, :staggered_grid_river_input; relative_widths = [6, 4, 2, 2, 2, 2])
 ```
 
 #### Static or cyclic/forcing input
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_river_static_cyclic_forcing_input; relative_widths = [6, 4, 2, 2, 2, 2])
+generate_table(Wflow.Routing, :staggered_grid_river_static_cyclic_forcing_input; relative_widths = [6, 4, 2, 2, 2, 2])
 ```
 
 #### States
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_river_state; relative_widths = [3, 3, 1, 1, 2])
+generate_table(Wflow.Routing, :staggered_grid_river_state; relative_widths = [3, 3, 1, 1, 2])
 ```
 
 #### Output
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_river_output; relative_widths = [3, 3, 1, 1, 2])
+generate_table(Wflow.Routing, :staggered_grid_river_output; relative_widths = [3, 3, 1, 1, 2])
 ```
 
 ### 1D floodplain flow
 #### Input
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_floodplain_1D_flow_input; relative_widths = [6, 4, 2, 2, 2, 2])
+generate_table(Wflow.Routing, :staggered_grid_floodplain_1D_flow_input; relative_widths = [6, 4, 2, 2, 2, 2])
 ```
 
 #### States
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_floodplain_1D_flow_state; relative_widths = [3, 3, 1, 1, 2])
+generate_table(Wflow.Routing, :staggered_grid_floodplain_1D_flow_state; relative_widths = [3, 3, 1, 1, 2])
 ```
 
 #### Output
 ```{julia}
 # | code-fold: true
-generate_table(Wflow.Routing, :local_inertial_floodplain_1D_flow_output; relative_widths = [3, 3, 1, 1, 2])
+generate_table(Wflow.Routing, :staggered_grid_floodplain_1D_flow_output; relative_widths = [3, 3, 1, 1, 2])
 ```
 
-### Overland flow
+### Local inertial overland flow (2D)
 #### Input
 ```{julia}
 # | code-fold: true
diff --git a/docs/model_docs/routing/kinwave.qmd b/docs/model_docs/routing/kinwave.qmd
index e28bd5d9e..147bb3736 100644
--- a/docs/model_docs/routing/kinwave.qmd
+++ b/docs/model_docs/routing/kinwave.qmd
@@ -76,6 +76,32 @@ kinematic_wave__adaptive_time_step_flag = true    # optional, default is false
 [Reservoir](./reservoirs.qmd) models representing reservoirs, (natural) lakes or other water
 storage features can be included as part of the river kinematic wave network.
 
+## Floodplain routing
+It is possible to include a [floodplain
+schematization](./staggered-grid.qmd#floodplain-schematization) as part of the kinematic
+wave river flow routing, using the following setting in the TOML file:
+
+```toml
+[model]
+floodplain_1d__flag = true  # include 1D floodplain schematization (default = false)
+```
+
+At the start of each time step, channel-floodplain water exchange, assumed to be
+instantaneous, is computed based on total water volume and the channel storage capacity. The
+exchange is added to the lateral inflow term of the kinematic wave equation. Routing is done
+separately for the river channel and floodplain. Floodplain water is routed down slope using
+the `accucapacityflux` function with a transport capacity based on Manning's equation:
+
+$$
+Q_{t+\Delta t} = n^{-1} R_t^{2/3} \sqrt{c_\mathrm{floodplain\ slope}} A_t
+$$
+
+where $\SIb{Q_{t+\Delta t}}{m^3\ s^{-1}}$ is the floodplain flow at time step $t+\Delta t$,
+$\SIb{n}{m^{-\frac{1}{3}}\ s}$ is the Manning's roughness coefficient, $\SIb{R_t}{m}$ is the
+hydraulic radius at the previous time step,  $c_\mathrm{floodplain\ slope}$ is the
+floodplain bed slope and $\SIb{A_t}{m^2}$ is the cross sectional flow area at the previous
+time step.
+
 ## External inflows
 External inflows $\SIb{}{m^3\ s^{-1}}$ for example water supply or abstractions can be added
 to the kinematic wave for river flow routing, as a static or cyclic parameter or as part of
diff --git a/docs/model_docs/routing/local-inertial.qmd b/docs/model_docs/routing/staggered-grid.qmd
similarity index 69%
rename from docs/model_docs/routing/local-inertial.qmd
rename to docs/model_docs/routing/staggered-grid.qmd
index 7a69231df..ea04f399a 100644
--- a/docs/model_docs/routing/local-inertial.qmd
+++ b/docs/model_docs/routing/staggered-grid.qmd
@@ -1,8 +1,11 @@
 ---
-title: Local inertial
+title: Staggered grid
 ---
 
+This section provides an overview of available flow routing schemes on a staggered grid.
+
 ## River and floodplain routing
+### Local inertial approximation
 The local inertial approximation of shallow water flow neglects only the convective
 acceleration term in the Saint-Venant momentum conservation equation. The numerical solution of
 the local inertial approximation on a staggered grid is as follows (Bates et al., 2010):
@@ -18,17 +21,17 @@ time step, $\SIb{Q_t}{m^3\ s^{-1}}$ is the river flow at the previous time step,
 the water surface slope at the previous time step and $\SIb{n}{m^{-\frac{1}{3}}\ s}$ is the
 Manning's roughness coefficient.
 
-The momentum equation is applied to each link between two river grid cells, while the
-continuity equation over $\Delta t$ is applied to each river cell:
+The momentum equation is applied to each edge between two river nodes, while the continuity
+equation over $\Delta t$ is applied to each river node:
 
 $$
 h^{t+\Delta t} = h^t + \Delta t \frac{\subtext{Q^{t+\Delta t}}{src} - \subtext{Q^{t+\Delta t}}{dst}}{A}
 $$
 
-where $\SIb{h^{t+\Delta t}}{m}$ is the water depthat time step $t+\Delta t$, $\SIb{h^t}{m}$
+where $\SIb{h^{t+\Delta t}}{m}$ is the water depth at time step $t+\Delta t$, $\SIb{h^t}{m}$
 is the water depth at the previous time step, $\SIb{A}{m^2}$ is the river area and
 $\SIb{\subtext{Q}{src}}{m^3\ s^{-1}}$ and $\SIb{\subtext{Q}{dst}}{m^3\ s^{-1}}$ represent
-river flow at the upstream and downstream link of the river cell, respectively.
+river flow at the upstream and downstream edge of the river node, respectively.
 
 The model time step $\Delta t$ for the local inertial model is estimated based on the
 Courant-Friedrichs-Lewy condition (Bates et al., 2010):
@@ -37,8 +40,8 @@ $$
 \Delta t = \alpha \min\left(\frac{\Delta x_i}{\sqrt{gh_i}}\right)
 $$
 
-where $\sqrt{gh_i}$ is the wave celerity for river cell $i$, $\SIb{\Delta x_i}{m}$ is the
-river length for river cell $i$ and $\alpha$ is a coefficient (typically between $0.2$ and
+where $\sqrt{gh_i}$ is the wave celerity for river node $i$, $\SIb{\Delta x_i}{m}$ is the
+river length for river node $i$ and $\alpha$ is a coefficient (typically between $0.2$ and
 $0.7$) to enhance the stability of the simulation.
 
 In the TOML file the following properties related to the local inertial model can be provided
@@ -66,10 +69,7 @@ model_boundary_condition_river_bank_water__depth = "riverdepth_bc"  # optional r
 These boundary conditions and river bed elevation and Manning's roughness coefficient are
 copied to a ghost node (downstream of the river outlet node) in the code.
 
-The optional 1D floodplain schematization is based on provided flood volumes as a function of
-flood depth (per flood depth interval) for each river cell. Wflow calculates from these flood
-volumes a rectangular floodplain profile for each flood depth interval. Routing is done
-separately for the river channel and floodplain.
+It is possible to include a [floodplain schematization] for this routing type.
 
 The momentum equation is most stable for low slope environments, and to keep the simulation
 stable for (partly) steep environments the `river_water_flow__froude_limit_flag` option is
@@ -87,15 +87,68 @@ the water velocity from the local inertial model is causing the Froude number to
 than $1.0$ , the water velocity (and flow) is reduced in order to maintain a Froude number of
 $1.0$.
 
+### Using Manning's equation
+The use of Manning's equation (assuming that the topography controls water flow mostly) for
+river and floodplain routing on a staggered grid is suitable for steep terrain where
+backwater effects are negligible:
+
+
+$$
+Q_{t+\Delta t} = n^{-1} R_t^{2/3} \sqrt{c_\mathrm{river\ slope}} A_t
+$$
+
+where $\SIb{Q_{t+\Delta t}}{m^3\ s^{-1}}$ is the river flow at time step $t+\Delta t$,
+$\SIb{n}{m^{-\frac{1}{3}}\ s}$ is the Manning's roughness coefficient, $\SIb{R_t}{m}$ is the
+hydraulic radius at the previous time step,  $c_\mathrm{river\ slope}$ is the river bed
+slope between nodes and $\SIb{A_t}{m^2}$ is the cross sectional flow area at the previous
+time step.
+
+As for the [local inertial model](#local-inertial-approximation), the momentum equation is
+applied to each edge between two river nodes, while the continuity equation over $\Delta t$
+is applied to each river node.
+
+The model time step $\Delta t$ for Manning's equation on a staggered grid is estimated as
+follows assuming a wide rectangular river channel:
+
+$$
+\begin{gathered}
+c_i = \frac{5}{3} n^{-1} R^{2/3} \sqrt{c_\mathrm{river\ slope}},\\
+\Delta t = \alpha \min\left(\frac{\Delta x_i}{c_i}\right)
+\end{gathered}
+$$
+
+where $c_i$ is the wave celerity for river edge $i$, $\SIb{\Delta x_i}{m}$ is the river
+length for river edge $i$ and $\alpha$ is a coefficient (typically between $0.2$ and $0.7$)
+to enhance the stability of the simulation.
+
+In the TOML file the following properties related to Manning's flow on a staggered grid can
+be provided for the `sbm` and `sbm_gwf` model types:
+
+```toml
+[model]
+river_routing = "manning_staggered"                   # default is "kinematic_wave"
+river_staggered_manning_flow__alpha_coefficient = 0.5 # alpha coefficient for model stability (default = 0.7)
+floodplain_1d__flag = true                            # include 1D floodplain schematization (default = false)
+```
+
+It is possible to include a [floodplain schematization] for this routing type.
+
+## Floodplain schematization
+The optional 1D floodplain schematization is based on provided flood volumes as a function
+of flood depth (per flood depth interval) for each river node. Wflow calculates from these
+flood volumes a rectangular floodplain profile for each flood depth interval. Routing is
+done separately for the river channel and floodplain.
+
+## Reservoirs
 [Reservoir](./reservoirs.qmd) models representing reservoirs, (natural) lakes or other water
-storage features can be included as part of the local inertial model for river flow (1D) and
-river and overland flow combined (see next section). Reservoir models are included as a
-boundary point with zero water depth for both river and overland flow. For river flow the
-reservoir model replaces the local inertial model at the reservoir location, and river flow
-$Q$ is set by the outflow from the reservoir. Overland flow at a reservoir location is not
-allowed to or from the downstream river grid cell.
-
-## Overland flow (2D)
+storage features can be included as part of 1D river flow and river and overland flow
+combined (see next section) on a staggered grid. Reservoir models are included as a boundary
+point with zero water depth for both local inertial river and overland flow routing. For
+river flow the reservoir model replaces the river flow model at the reservoir location, and
+river flow $Q$ is set by the outflow from the reservoir. Overland flow at a reservoir
+location is not allowed to or from the downstream river node.
+
+## Local inertial overland flow (2D)
 For the simulation of 2D overland flow on a staggered grid the numerical scheme proposed by de
 Almeida et al. (2012) is adopted. The explicit solution for the estimation of water discharge
 between two cells in the x-direction is of the following form (following the notation of
@@ -147,18 +200,18 @@ routing](#river-and-floodplain-routing) section of the local inertial model.
 
 ## External inflows
 External inflows $\SIb{}{m^3\ s^{-1}}$ for example water supply or abstractions can to added
-the local inertial model for river flow (1D) and river and overland flow combined (1D-2D),
-as a static or cyclic parameter or as part of forcing. For example, cyclic inflow can be
-provided in the TOML file as follows:
+the 1D river flow model on a staggered grid and local inertial river and overland flow
+combined (1D-2D), as a static or cyclic parameter or as part of forcing. For example, cyclic
+inflow can be provided in the TOML file as follows:
 
 ```toml
 [input.cyclic]
 river_water__external_inflow_volume_flow_rate = "river_inflow"
 ```
-These inflows are added or subtracted as part of the continuity equation of the local
-inertial model. Abstractions are limited to a maximum of $\SI{80}{\%}$ of the river storage.
-Negative external inflows are satisfied first, and the remaining river storage is available
-for internal abstractions as part of water demand and allocation computations.
+These inflows are added or subtracted as part of the continuity equation of river flow model
+on a staggered grid. Abstractions are limited to a maximum of $\SI{80}{\%}$ of the river
+storage. Negative external inflows are satisfied first, and the remaining river storage is
+available for internal abstractions as part of water demand and allocation computations.
 
 ## Abstractions
 Internal abstractions $\SIb{}{m^3\ s^{-1}}$ from the river are possible when water demand and
diff --git a/docs/user_guide/faq.qmd b/docs/user_guide/faq.qmd
index 73b80ed2f..245d71904 100644
--- a/docs/user_guide/faq.qmd
+++ b/docs/user_guide/faq.qmd
@@ -34,7 +34,14 @@ output](./toml_file.qmd#output-netcdf-section).
 
 The choice of a specific river routing option can vary depending on the model and use case.
 However, the numerical properties of the routing schemes provide an indication of their
-advantages and disadvantages.
+advantages and disadvantages. The following river routing schemes are available:
+
+- Kinematic wave
+- Mannings flow on a staggered grid
+- Local inertial routing on a staggered grid
+
+These routing schemes support an optional 1D floodplain schematization (routing is done
+separately for the river channel and floodplain).
 
 The kinematic wave method is driven by the river slope within each cell, assuming that the
 water surface slope is parallel to the bed slope. This results in sharper discharge peaks
@@ -43,6 +50,10 @@ approach is especially suitable for steep (upstream) areas where flow propagatio
 dominated by topography. The kinematic wave approach does not include any form of backwater
 effects.
 
+As the kinematic wave method, the use of Manning's equation on a staggered grid is also
+driven by the river slope and thus most suitable for steep terrain where backwater effects
+are negligible.
+
 On the other hand, the local inertial method incorporates the slope of the water surface
 into the momentum equation. This results in 'damping effects' on flow propagating through
 the cells and yields better results in scenarios where the water surface slope differs from
@@ -51,19 +62,16 @@ slope is limited or during the propagation of larger flood waves where the water
 slope is greater than the river slope itself. The local inertial method does incorporate
 backwater effects, although not as comprehensively as the full dynamic wave equation.
 
-Generally, the kinematic wave method is computationally faster than the local inertial
-method for river (and land) routing. For the multi-threading execution of the kinematic
-wave, solved with a nonlinear scheme using Newton's method, the order of execution
-(sub-basins) is important (from upstream to downstream sub-basin). For the momentum equation
-of the local inertial approach an explicit numerical method is used and can be solved
-independently for each cell during multi-threading execution. As a consequence, speedups for
-local inertial routing are larger than for kinematic wave routing using multi-threading
-(compared to a single thread run) and the difference between local inertial and kinematic
-wave routing run times gets smaller as the number of threads increase.
-
-Currently, only the local inertial approach supports the inclusion of floodplains in river
-routing, see also [this
-question](./faq.qmd#what-is-the-difference-between-1d-2d-and-no-floodplains).
+Generally, the kinematic wave method and the use of Manning's equation on a staggered grid
+is computationally faster than the local inertial method for river (and land) routing. For
+the multi-threading execution of the kinematic wave, solved with a nonlinear scheme using
+Newton's method, the order of execution (sub-basins) is important (from upstream to
+downstream sub-basin). For the momentum equation of river flow on a staggered grid an
+explicit numerical method is used and can be solved independently for each cell during
+multi-threading execution. As a consequence, speedups for river flow on a staggered grid are
+larger than for kinematic wave routing using multi-threading (compared to a single thread
+run) and the difference between river flow on a staggered grid and kinematic wave routing
+run times gets smaller as the number of threads increase.
 
 ### Which land routing option should I choose?
 
@@ -85,11 +93,11 @@ combination of kinematic wave for land routing and local inertial for river rout
 
 Effects of floodplain flow can be included in several ways. A one-dimensional sub-grid
 approximation (hence the name 1D floodplains) can be included when kinematic wave land
-routing is combined with local inertial river routing. The water surface elevations of the
-river channel and the floodplain are the same and it is assumed that water is exchanged
-instantaneously between the river channel and the floodplain. When the water depth in the
-river channel rises above the river bankfull depth, water is exchanged between the river
-channel and the floodplain resulting in an attenuation of the discharge peaks.
+routing is combined with river routing. The water surface elevations of the river channel
+and the floodplain are the same and it is assumed that water is exchanged instantaneously
+between the river channel and the floodplain. When the water depth in the river channel
+rises above the river bankfull depth, water is exchanged between the river channel and the
+floodplain resulting in an attenuation of the discharge peaks.
 
 ![Schematization of the 1D floodplain concept (side view)](../images/floodplains_1d.png)
 
@@ -110,10 +118,9 @@ the 1D floodplain method can be used instead.
 
 ![Schematization of the 2D floodplain concept (top view)](../images/floodplains_2d.png)
 
-It is also possible to run the model without 1D floodplains. This means that the model is
-using kinematic wave for land routing and either kinematic wave or local inertial for river
-routing. This can lead to unrealistically large water depths when the bankfull capacity is
-exceeded as there is no exchange of water between the river channel and floodplains.
+It is also possible to run the model without 1D floodplains. This can lead to
+unrealistically large water depths when the bankfull capacity is exceeded as there is no
+exchange of water between the river channel and floodplains.
 
 ### How do I run only a part of my model?
 You can easily run only a subsection of your model, by modifying the
diff --git a/docs/user_guide/model_config.qmd b/docs/user_guide/model_config.qmd
index 157eaff0b..6f6180312 100644
--- a/docs/user_guide/model_config.qmd
+++ b/docs/user_guide/model_config.qmd
@@ -59,13 +59,12 @@ reservoir_water_surface__elevation = "waterlevel_reservoir"
 ```
 
 ## Enabling floodplain routing
-As part of the local inertial model for river flow.
 
-### 1D floodplains
+### 1D floodplains on a staggered grid
 
 ```toml
 [model]
-river_routing = "local_inertial"
+river_routing = "local_inertial"    # one of ("local_inertial" or "manning_staggered")
 floodplain_1d__flag = true
 
 [input.static]
@@ -77,7 +76,22 @@ floodplain_water__instantaneous_volume_flow_rate = "q_floodplain"
 floodplain_water__depth = "h_floodplain"
 ```
 
-### 2D floodplains
+### 1D floodplains (kinematic wave routing)
+
+```toml
+[model]
+river_routing = "kinematic_wave"
+floodplain_1d__flag = true
+
+[input.static]
+floodplain_water__sum_of_volume_per_depth = "floodplain_volume"
+floodplain_water_flow__manning_n_parameter = "floodplain_n"
+
+[state.variables]
+floodplain_water__volume = "storage_floodplain"
+```
+
+### 2D floodplains (local inertial routing)
 
 ```toml
 [model]
diff --git a/docs/user_guide/toml_file.qmd b/docs/user_guide/toml_file.qmd
index 0d0ed35d9..07f3b20ad 100644
--- a/docs/user_guide/toml_file.qmd
+++ b/docs/user_guide/toml_file.qmd
@@ -75,7 +75,7 @@ Model-specific settings can be included in the model section of the TOML file.
 ```toml
 [model]
 type = "sbm"                                    # one of ("sbm" or "sbm_gwf")
-snow_gravitational_transport__flag = true         # include lateral snow transport in the model, default is false
+snow_gravitational_transport__flag = true       # include lateral snow transport in the model, default is false
 snow__flag = true                               # include snow modelling, default is false
 cold_start__flag = true                         # cold (cold_start__flag = true) or warm state (cold_start__flag = false), default is true
 reservoir__flag = true                          # include reservoir modelling, default is false
diff --git a/utils/export_parameter_data/routing_metadata.json b/utils/export_parameter_data/routing_metadata.json
index 0632ba86e..16b8e105a 100644
--- a/utils/export_parameter_data/routing_metadata.json
+++ b/utils/export_parameter_data/routing_metadata.json
@@ -1,4 +1,13 @@
 {
+  "floodplain__slope": {
+    "unit": "m m-1",
+    "default": null,
+    "fill": null,
+    "type": "Float64",
+    "description": "Floodplain slope",
+    "allow_missing": false,
+    "allow_dynamic_input": false
+  },
   "floodplain_water__depth": {
     "unit": "m",
     "default": null,