Support kinematic wave floodplain routing

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")

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"

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

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]

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,20 +288,24 @@ 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,
 )
     (; 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) =

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

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
 

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