Support kinematic wave floodplain routing (#870)

* Add `floodplain` to kinematic river flow
Initialized as `nothing`.

* Move `floodplain` to separate file

* Add `bankfull_storage` to kinematic wave river

* Separate init kinwave and local inertial floodplain

* Refactor surface flow routing
Introduce one struct for riverflow routing `RiverFlowModel` and one struct for overlandflow routing `OverlandFlowModel`. The routing method (subtype of `AbstractRoutingMethod`) is used to differentiate between different routing approaches for river and overland flow. Extending to other routing methods is now easier.

* Change abstract surface flow routing types
Remove `AbstractRoutingMethod` from these abstract types and add the routing method as first field to `RiverFlowModel` and `OverlandFlowModel`. This simplifies dispatching a bit and looks a bit cleaner.

* Update docstrings

* Renaming of init functions surface flow routing
Also add `slope` parameter that can be used by kinematic wave routing (as alternative to local inertial routing) using a staggered scheme.

* Refactor init local inertial river flow parameters
Use a separate function for reading river parameters that can also be used by kinematic wave routing with a staggered scheme.

* Update info log message

* Add init kinwave staggered river flow

* Refactor init river flow staggered a bit

* Use split up of local_inertial_river_update! function
As implemented in PR#835.

* Init parameter `h_thresh` staggered kinematic wave

* Fix typo

* Fix slope computation

* Add update manning's river flow staggered grid
Using the term kinematic wave is not entirely correct as the Manning's equation is used as an alternative for the local inertial approach on a staggered grid.

* Refactor init staggered routing

* Add `stable_timestep` function `ManningStaggered`
For manning flow on a staggered grid.

* Some fixes
To run manning river and floodplain flow on a staggered grid.

* Revert `river_routing` TOML setting change

* Add TOML setting `alpha` staggered manning flow

* Rename file `surface_local_inertial.jl`

* Add test river and floodplain Manning's flow on staggered grid

* Add missing toml file

* Some ordering and renaming
Edge/node variables and parameters.

* Add floodplain model for kinematic wave

* Simplify dispatching on floodplain model

* Init manning flow parameters
Remove update of `alpha` at each model timestep as dynamic input of manning roughness is not allowed.

* Use separate structs river and overland flow variables
For kinematic wave routing.

* Add kinematic wave flood routing
For routing type `KinematicWave` using a compound channel (not separate routing for floodplain) with a compound `alpha` value.

* Use function for wetted perimeter channel

* Refactor computation of floodplain variables

* Use function for active floodplain cells

* Move update kinematic wave floodplain to function

* Separate routing for floodplain
Making use of a simple routing scheme ( Manning flow and `accucapacityflux`) as kinematic wave routing gets quite slow for very small flows.

* Support mass balance computation `ManningStaggered`

* Consistent use of `at_edge`
For parameter and variable names stored at edges that normally are stored at cell centres/nodes.

* Explain `FloodPlainModel` `domain`

* Address review comment
Make use of const `LDD_PIT`.

* Address review comment

* Remove `Float64` `hydraulic_radius_pow`

* Address review comment
Remove `pow` from `celerity` computation.

* Use constant `BETA_KINWAVE` for `celerity` computation

* Compute stable timestep using variables at edges
For `RiverFlowModel{<:ManningStaggered}`.

* Remove magic number
In `stable_timestep` functions.

* Address review comment

* Improve docstring for `FloodPlainModel`

* Add docstrings related to `FloodPlainModel`

* Update and add docstrings
For routing on a staggered grid. Also removed `cell` from functions, for a staggered grid typically `node` and `edge` terms are used.

* Some more docstrings

* Update `stable_timestep` staggered scheme

* Add small comment

* Add docstring and renaming function

* Update docs routing on staggered grid

* Fix links staggered grid docs

* Add docs floodplain routing kinematic wave

* Update docs and introduce standard name `floodplain__slope`
To support different slope values for river channel and floodplain for `kinematic_wave` routing.

* Update routing_metadata json output

* Update model_config docs

* Small formatting

* Update faq docs

* Remove flow_area from test description

* Use `floodplain__slope` for kinematic wave floodplain

* Add unit test river flow with floodplain
For river flow on a staggered grid using Manning's equation.

* Small change comment in unit test

* Add unit test kinematic wave routing with floodplain

* Update model settings docs

* Rename inds `EdgeConnectivity`
Keeping `x` and `y` dir in the name, bit more generic than `left`/`right`.

* Update struct model docs

* Update changelog

* Update docs model parameters staggered grid

* Fix some doc tags model parameters

* Formatting and typo docs routing parameters

* Fix missing space in changelog

* Rename routing method kinematic wave floodplain
From `KinematicWave` to `Manning`, a better description as it uses the Manning equation and not the full kinematic wave equation.

* Fix `NaN`/`Inf` values floodplain hydraulic radius
Removed flow area and hydraulic radius fields of local inertial structs as these are used only internally (per function) and are not part of model output. Added check for hydraulic radius floodplain.

* Remove check active floodplain index
Already part of release/v1.0 branch: this check was based on `water_depth_at_edge` of river nodes and should be based on `water_depth_at_edge` of floodplain nodes. However, testing showed there is no performance benefit from this check and looping over active floodplain indices.

* Fix large `q` fluctuations `manning_staggered` routing
By using the "upwind"option (instead of "max" as used for local inertial routing) for water level and depth estimation at edge.

* Address review comment

* Small change "upwind" option `manning_staggered`
For zero flows the maximum waterlevel between two adjacent cells was computed. For slope based mannings flow (one direction) this is not required.

* Address review comment
Fix floodplain states check.

* Fix states kinematic wave flood routing
The floodplain water depth was incorrect as it was based on the floodplain storage excluding river channel storage using the floodplain profile, while river channel storage should be included. This is now resolved by using the floodplain storage as state instead of the water depth. Additonally, for floodplain routing a check if the downstream floodplain is active (based on flow area) has been added.

* Update docs Model structure

* Update docs state kinematic wave floodplain

* Update tests

* Remove `hydraulic_radius` check floodplain

* Fix `dt` input `update_floodplain_model!`

* Fix reading `alpha_coefficient` `manning_staggered`

* Revert small change used for testing

---------

Co-authored-by: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com>

Author: vers-w

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