diff --git a/Manifest.toml b/Manifest.toml
index 4b4538c67..42cbf3109 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.0"
 manifest_format = "2.0"
-project_hash = "7e18cea33bf679231c5c94c6b99a1b635f9a1b61"
+project_hash = "742e16f5de5af91dfff30d91944981a5c7673d44"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245"
diff --git a/core/src/allocation.jl b/core/src/allocation.jl
index ac6a0745f..f8248bbc3 100644
--- a/core/src/allocation.jl
+++ b/core/src/allocation.jl
@@ -1,19 +1,45 @@
+"""Find the edges from the main network to a subnetwork."""
+function find_subnetwork_connections!(p::Parameters)::Nothing
+    (; allocation, graph, user) = p
+    n_priorities = length(user.demand[1])
+    (; subnetwork_demands, subnetwork_allocateds) = allocation
+    for node_id in graph[].node_ids[1]
+        for outflow_id in outflow_ids(graph, node_id)
+            if graph[outflow_id].allocation_network_id != 1
+                main_network_source_edges =
+                    get_main_network_connections(p, graph[outflow_id].allocation_network_id)
+                edge = (node_id, outflow_id)
+                push!(main_network_source_edges, edge)
+                subnetwork_demands[edge] = zeros(n_priorities)
+                subnetwork_allocateds[edge] = zeros(n_priorities)
+            end
+        end
+    end
+    return nothing
+end
+
 """
 Find all nodes in the subnetwork which will be used in the allocation network.
 Some nodes are skipped to optimize allocation optimization.
 """
 function allocation_graph_used_nodes!(p::Parameters, allocation_network_id::Int)::Nothing
-    (; graph, basin, fractional_flow) = p
+    (; graph, basin, fractional_flow, allocation) = p
+    (; main_network_connections) = allocation
 
     node_ids = graph[].node_ids[allocation_network_id]
     used_nodes = Set{NodeID}()
     for node_id in node_ids
+        use_node = false
         has_fractional_flow_outneighbors =
             get_fractional_flow_connected_basins(node_id, basin, fractional_flow, graph)[3]
         node_type = graph[node_id].type
         if node_type in [:user, :basin, :terminal]
-            push!(used_nodes, node_id)
+            use_node = true
         elseif has_fractional_flow_outneighbors
+            use_node = true
+        end
+
+        if use_node
             push!(used_nodes, node_id)
         end
     end
@@ -21,13 +47,24 @@ function allocation_graph_used_nodes!(p::Parameters, allocation_network_id::Int)
     # Add nodes in the allocation graph for nodes connected to the source edges
     # One of these nodes can be outside the subnetwork, as long as the edge
     # connects to the subnetwork
-    for edge_metadata in graph[].edges_source[allocation_network_id]
+    edges_source = graph[].edges_source
+    for edge_metadata in get(edges_source, allocation_network_id, Set{EdgeMetadata}())
         (; from_id, to_id) = edge_metadata
         push!(used_nodes, from_id)
         push!(used_nodes, to_id)
     end
 
     filter!(in(used_nodes), node_ids)
+
+    # For the main network, include nodes that connect the main network to a subnetwork
+    # (also includes nodes not in the main network in the input)
+    if is_main_network(allocation_network_id)
+        for connections_subnetwork in main_network_connections
+            for connection in connections_subnetwork
+                union!(node_ids, connection)
+            end
+        end
+    end
     return nothing
 end
 
@@ -241,6 +278,8 @@ function process_allocation_graph_edges!(
     return capacity
 end
 
+const allocation_source_nodetypes = Set{Symbol}([:level_boundary, :flow_boundary])
+
 """
 The source nodes must only have one allocation outneighbor and no allocation inneighbors.
 """
@@ -251,24 +290,23 @@ function valid_sources(p::Parameters, allocation_network_id::Int)::Bool
 
     errors = false
 
-    for (id_source, id_dst) in edge_ids
+    for edge in edge_ids
+        (id_source, id_dst) = edge
         if graph[id_source, id_dst].allocation_network_id_source == allocation_network_id
-            ids_allocation_in = [
-                label for label in inneighbor_labels(graph, id_source) if
-                graph[label, id_source].allocation_flow
-            ]
-            if length(ids_allocation_in) !== 0
-                errors = true
-                @error "Source edge ($id_source, $id_dst) is not an entry point of subnetwork $allocation_network_id"
-            end
+            from_source_node = graph[id_source].type in allocation_source_nodetypes
 
-            ids_allocation_out = [
-                label for label in outneighbor_labels(graph, id_source) if
-                graph[id_source, label].allocation_flow
-            ]
-            if length(ids_allocation_out) !== 1
-                errors = true
-                @error "Source edge ($id_source, $id_dst) is not the only allocation edge coming from $id_source"
+            if is_main_network(allocation_network_id)
+                if !from_source_node
+                    errors = true
+                    @error "The source node of source edge $edge in the main network must be one of $allocation_source_nodetypes."
+                end
+            else
+                from_main_network = is_main_network(graph[id_source].allocation_network_id)
+
+                if !from_source_node && !from_main_network
+                    errors = true
+                    @error "The source node of source edge $edge for subnetwork $allocation_network_id is neither a source node nor is it coming from the main network."
+                end
             end
         end
     end
@@ -297,6 +335,25 @@ function avoid_using_own_returnflow!(p::Parameters, allocation_network_id::Int):
     return nothing
 end
 
+"""
+Add the edges connecting the main network work to a subnetwork to both the main network
+and subnetwork allocation graph.
+"""
+function add_subnetwork_connections!(p::Parameters, allocation_network_id::Int)::Nothing
+    (; graph, allocation) = p
+    (; main_network_connections) = allocation
+    edge_ids = graph[].edge_ids[allocation_network_id]
+
+    if is_main_network(allocation_network_id)
+        for connections in main_network_connections
+            union!(edge_ids, connections)
+        end
+    else
+        union!(edge_ids, get_main_network_connections(p, allocation_network_id))
+    end
+    return nothing
+end
+
 """
 Build the graph used for the allocation problem.
 """
@@ -312,6 +369,7 @@ function allocation_graph(
 
     # Process the edges in the allocation graph
     process_allocation_graph_edges!(capacity, edges_composite, p, allocation_network_id)
+    add_subnetwork_connections!(p, allocation_network_id)
 
     if !valid_sources(p, allocation_network_id)
         error("Errors in sources in allocation graph.")
@@ -351,10 +409,21 @@ function add_variables_absolute_value!(
     allocation_network_id::Int,
     config::Config,
 )::Nothing
-    (; graph) = p
-    node_ids = graph[].node_ids[allocation_network_id]
-    node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user]
+    (; graph, allocation) = p
+    (; main_network_connections) = allocation
     if startswith(config.allocation.objective_type, "linear")
+        node_ids = graph[].node_ids[allocation_network_id]
+        node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user]
+
+        # For the main network, connections to subnetworks are treated as users
+        if is_main_network(allocation_network_id)
+            for connections_subnetwork in main_network_connections
+                for connection in connections_subnetwork
+                    push!(node_ids_user, connection[2])
+                end
+            end
+        end
+
         problem[:F_abs] = JuMP.@variable(problem, F_abs[node_id = node_ids_user])
     end
     return nothing
@@ -375,11 +444,12 @@ function add_constraints_capacity!(
     allocation_network_id::Int,
 )::Nothing
     (; graph) = p
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
     F = problem[:F]
     edge_ids = graph[].edge_ids[allocation_network_id]
     edge_ids_finite_capacity = Tuple{NodeID, NodeID}[]
     for edge in edge_ids
-        if !isinf(capacity[edge...])
+        if !isinf(capacity[edge...]) && edge ∉ main_network_source_edges
             push!(edge_ids_finite_capacity, edge)
         end
     end
@@ -461,13 +531,19 @@ function add_constraints_flow_conservation!(
     p::Parameters,
     allocation_network_id::Int,
 )::Nothing
-    (; graph) = p
+    (; graph, allocation) = p
     F = problem[:F]
     node_ids = graph[].node_ids[allocation_network_id]
-    node_ids_basin = [node_id for node_id in node_ids if graph[node_id].type == :basin]
+    node_ids_conservation =
+        [node_id for node_id in node_ids if graph[node_id].type == :basin]
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
+    for edge in main_network_source_edges
+        push!(node_ids_conservation, edge[2])
+    end
+    unique!(node_ids_conservation)
     problem[:flow_conservation] = JuMP.@constraint(
         problem,
-        [node_id = node_ids_basin],
+        [node_id = node_ids_conservation],
         sum([
             F[(node_id, outneighbor_id)] for
             outneighbor_id in outflow_ids_allocation(graph, node_id)
@@ -523,12 +599,23 @@ function add_constraints_absolute_value!(
     allocation_network_id::Int,
     config::Config,
 )::Nothing
-    (; graph) = p
-    node_ids = graph[].node_ids[allocation_network_id]
+    (; graph, allocation) = p
+    (; main_network_connections) = allocation
 
     objective_type = config.allocation.objective_type
     if startswith(objective_type, "linear")
+        node_ids = graph[].node_ids[allocation_network_id]
         node_ids_user = [node_id for node_id in node_ids if graph[node_id].type == :user]
+
+        # For the main network, connections to subnetworks are treated as users
+        if is_main_network(allocation_network_id)
+            for connections_subnetwork in main_network_connections
+                for connection in connections_subnetwork
+                    push!(node_ids_user, connection[2])
+                end
+            end
+        end
+
         node_ids_user_inflow = Dict(
             node_id_user => only(inflow_ids_allocation(graph, node_id_user)) for
             node_id_user in node_ids_user
@@ -689,6 +776,60 @@ function AllocationModel(
     )
 end
 
+"""
+Add a term to the expression of the objective function corresponding to
+the demand of a user.
+"""
+function add_user_term!(
+    ex::Union{JuMP.QuadExpr, JuMP.AffExpr},
+    edge::Tuple{NodeID, NodeID},
+    objective_type::Symbol,
+    demand::Float64,
+    model::AllocationModel,
+)::Nothing
+    (; problem) = model
+    F = problem[:F]
+    F_edge = F[edge]
+    node_id_user = edge[2]
+
+    if objective_type == :quadratic_absolute
+        # Objective function ∑ (F - d)^2
+        JuMP.add_to_expression!(ex, 1, F_edge, F_edge)
+        JuMP.add_to_expression!(ex, -2 * demand, F_edge)
+        JuMP.add_to_expression!(ex, demand^2)
+
+    elseif objective_type == :quadratic_relative
+        # Objective function ∑ (1 - F/d)^2
+        if demand ≈ 0
+            return nothing
+        end
+        JuMP.add_to_expression!(ex, 1.0 / demand^2, F_edge, F_edge)
+        JuMP.add_to_expression!(ex, -2.0 / demand, F_edge)
+        JuMP.add_to_expression!(ex, 1.0)
+
+    elseif objective_type == :linear_absolute
+        # Objective function ∑ |F - d|
+        JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -demand)
+        JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], demand)
+
+    elseif objective_type == :linear_relative
+        # Objective function ∑ |1 - F/d|
+        JuMP.set_normalized_coefficient(
+            problem[:abs_positive][node_id_user],
+            F_edge,
+            iszero(demand) ? 0 : 1 / demand,
+        )
+        JuMP.set_normalized_coefficient(
+            problem[:abs_negative][node_id_user],
+            F_edge,
+            iszero(demand) ? 0 : -1 / demand,
+        )
+    else
+        error("Invalid allocation objective type $objective_type.")
+    end
+    return nothing
+end
+
 """
 Set the objective for the given priority.
 For an objective with absolute values this also involves adjusting constraints.
@@ -700,8 +841,9 @@ function set_objective_priority!(
     priority_idx::Int,
 )::Nothing
     (; objective_type, problem, allocation_network_id) = allocation_model
-    (; graph, user) = p
+    (; graph, user, allocation) = p
     (; demand, node_id) = user
+    (; main_network_connections, subnetwork_demands) = allocation
     edge_ids = graph[].edge_ids[allocation_network_id]
 
     F = problem[:F]
@@ -713,6 +855,18 @@ function set_objective_priority!(
 
     demand_max = 0.0
 
+    # Terms for subnetworks as users
+    if is_main_network(allocation_network_id)
+        for connections_subnetwork in main_network_connections
+            for connection in connections_subnetwork
+                d = subnetwork_demands[connection][priority_idx]
+                demand_max = max(demand_max, d)
+                add_user_term!(ex, connection, objective_type, d, allocation_model)
+            end
+        end
+    end
+
+    # Terms for user nodes
     for edge_id in edge_ids
         node_id_user = edge_id[2]
         if graph[node_id_user].type != :user
@@ -722,43 +876,7 @@ function set_objective_priority!(
         user_idx = findsorted(node_id, node_id_user)
         d = demand[user_idx][priority_idx](t)
         demand_max = max(demand_max, d)
-        F_edge = F[edge_id]
-
-        if objective_type == :quadratic_absolute
-            # Objective function ∑ (F - d)^2
-            JuMP.add_to_expression!(ex, 1, F_edge, F_edge)
-            JuMP.add_to_expression!(ex, -2 * d, F_edge)
-            JuMP.add_to_expression!(ex, d^2)
-
-        elseif objective_type == :quadratic_relative
-            # Objective function ∑ (1 - F/d)^2S
-            if d ≈ 0
-                continue
-            end
-            JuMP.add_to_expression!(ex, 1.0 / d^2, F_edge, F_edge)
-            JuMP.add_to_expression!(ex, -2.0 / d, F_edge)
-            JuMP.add_to_expression!(ex, 1.0)
-
-        elseif objective_type == :linear_absolute
-            # Objective function ∑ |F - d|
-            JuMP.set_normalized_rhs(problem[:abs_positive][node_id_user], -d)
-            JuMP.set_normalized_rhs(problem[:abs_negative][node_id_user], d)
-
-        elseif objective_type == :linear_relative
-            # Objective function ∑ |1 - F/d|
-            JuMP.set_normalized_coefficient(
-                problem[:abs_positive][node_id_user],
-                F_edge,
-                iszero(d) ? 0 : 1 / d,
-            )
-            JuMP.set_normalized_coefficient(
-                problem[:abs_negative][node_id_user],
-                F_edge,
-                iszero(d) ? 0 : -1 / d,
-            )
-        else
-            error("Invalid allocation objective type $objective_type.")
-        end
+        add_user_term!(ex, edge_id, objective_type, d, allocation_model)
     end
 
     # Add flow cost
@@ -788,35 +906,66 @@ function assign_allocations!(
     allocation_model::AllocationModel,
     p::Parameters,
     t::Float64,
-    priority_idx::Int,
+    priority_idx::Int;
+    collect_demands::Bool = false,
 )::Nothing
     (; problem, allocation_network_id) = allocation_model
-    (; graph, user) = p
+    (; graph, user, allocation) = p
+    (;
+        subnetwork_demands,
+        subnetwork_allocateds,
+        allocation_network_ids,
+        main_network_connections,
+    ) = allocation
     (; record) = user
     edge_ids = graph[].edge_ids[allocation_network_id]
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
     F = problem[:F]
     for edge_id in edge_ids
+        # If this edge is a source edge from the main network to a subnetwork,
+        # and demands are being collected, add its flow to the demand of this edge
+        if collect_demands &&
+           graph[edge_id...].allocation_network_id_source == allocation_network_id &&
+           edge_id ∈ main_network_source_edges
+            allocated = JuMP.value(F[edge_id])
+            subnetwork_demands[edge_id][priority_idx] += allocated
+        end
+
         user_node_id = edge_id[2]
-        if graph[user_node_id].type != :user
-            continue
+
+        if graph[user_node_id].type == :user
+            allocated = JuMP.value(F[edge_id])
+            user_idx = findsorted(user.node_id, user_node_id)
+            user.allocated[user_idx][priority_idx] = allocated
+
+            # Save allocations to record
+            push!(record.time, t)
+            push!(record.allocation_network_id, allocation_model.allocation_network_id)
+            push!(record.user_node_id, Int(user_node_id))
+            push!(record.priority, user.priorities[priority_idx])
+            push!(record.demand, user.demand[user_idx][priority_idx](t))
+            push!(record.allocated, allocated)
+            # TODO: This is now the last abstraction before the allocation update,
+            # should be the average abstraction since the last allocation solve
+            push!(
+                record.abstracted,
+                get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0),
+            )
+        end
+    end
+
+    # Write the flows to the subnetworks as allocated flows
+    # in the allocation object
+    if is_main_network(allocation_network_id)
+        for (allocation_network_id, main_network_source_edges) in
+            zip(allocation_network_ids, main_network_connections)
+            if is_main_network(allocation_network_id)
+                continue
+            end
+            for edge_id in main_network_source_edges
+                subnetwork_allocateds[edge_id][priority_idx] = JuMP.value(F[edge_id])
+            end
         end
-        user_idx = findsorted(user.node_id, user_node_id)
-        allocated = JuMP.value(F[edge_id])
-        user.allocated[user_idx][priority_idx] = allocated
-
-        # Save allocations to record
-        push!(record.time, t)
-        push!(record.allocation_network_id, allocation_model.allocation_network_id)
-        push!(record.user_node_id, Int(user_node_id))
-        push!(record.priority, user.priorities[priority_idx])
-        push!(record.demand, user.demand[user_idx][priority_idx](t))
-        push!(record.allocated, allocated)
-        # TODO: This is now the last abstraction before the allocation update,
-        # should be the average abstraction since the last allocation solve
-        push!(
-            record.abstracted,
-            get_flow(graph, inflow_id(graph, user_node_id), user_node_id, 0),
-        )
     end
     return nothing
 end
@@ -824,27 +973,43 @@ end
 """
 Adjust the source flows.
 """
-function adjust_source_flows!(
+function adjust_source_capacities!(
     allocation_model::AllocationModel,
     p::Parameters,
-    priority_idx::Int,
+    priority_idx::Int;
+    collect_demands::Bool = false,
 )::Nothing
     (; problem) = allocation_model
-    (; graph) = p
+    (; graph, allocation) = p
     (; allocation_network_id) = allocation_model
+    (; subnetwork_allocateds) = allocation
     edge_ids = graph[].edge_ids[allocation_network_id]
     source_constraints = problem[:source]
     F = problem[:F]
 
-    # It is assumed that the allocation procedure does not have to be differentiated.
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
+
     for edge_id in edge_ids
-        # If it is a source edge.
         if graph[edge_id...].allocation_network_id_source == allocation_network_id
+            # If it is a source edge for this allocation problem
             if priority_idx == 1
-                # Reset the source to the current flow.
+                # If the optimization was just started, i.e. sources have to be reset
+                if edge_id in main_network_source_edges
+                    if collect_demands
+                        # Set the source capacity to effectively unlimited if subnetwork demands are being collected
+                        source_capacity = Inf
+                    else
+                        # Set the source capacity to the value allocated to the subnetwork over this edge
+                        source_capacity = subnetwork_allocateds[edge_id][priority_idx]
+                    end
+                else
+                    # Reset the source to the current flow from the physical layer.
+                    source_capacity = get_flow(graph, edge_id..., 0)
+                end
                 JuMP.set_normalized_rhs(
                     source_constraints[edge_id],
-                    get_flow(graph, edge_id..., 0),
+                    # It is assumed that the allocation procedure does not have to be differentiated.
+                    source_capacity,
                 )
             else
                 # Subtract the allocated flow from the source.
@@ -876,11 +1041,15 @@ function adjust_edge_capacities!(
     constraints_capacity = problem[:capacity]
     F = problem[:F]
 
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
+
     for edge_id in edge_ids
         c = capacity[edge_id...]
 
-        # Edges with infinite capacity have no capacity constraints
-        if isinf(c)
+        # These edges have no capacity constraints:
+        # - With infinite capacity
+        # - Being a source from the main network to a subnetwork
+        if isinf(c) || edge_id ∈ main_network_source_edges
             continue
         end
 
@@ -902,18 +1071,34 @@ end
 Update the allocation optimization problem for the given subnetwork with the problem state
 and flows, solve the allocation problem and assign the results to the users.
 """
-function allocate!(p::Parameters, allocation_model::AllocationModel, t::Float64)::Nothing
-    (; user) = p
-    (; problem) = allocation_model
+function allocate!(
+    p::Parameters,
+    allocation_model::AllocationModel,
+    t::Float64;
+    collect_demands::Bool = false,
+)::Nothing
+    (; user, allocation) = p
+    (; problem, allocation_network_id) = allocation_model
     (; priorities) = user
+    (; subnetwork_demands) = allocation
 
     # TODO: Compute basin flow from vertical fluxes and basin volume.
     # Set as basin demand if the net flow is negative, set as source
     # in the flow_conservation constraints if the net flow is positive.
     # Solve this as a separate problem before the priorities below
 
+    main_network_source_edges = get_main_network_connections(p, allocation_network_id)
+
+    if collect_demands
+        for main_network_connection in keys(subnetwork_demands)
+            if main_network_connection in main_network_source_edges
+                subnetwork_demands[main_network_connection] .= 0.0
+            end
+        end
+    end
+
     for priority_idx in eachindex(priorities)
-        adjust_source_flows!(allocation_model, p, priority_idx)
+        adjust_source_capacities!(allocation_model, p, priority_idx; collect_demands)
 
         # Subtract the flows used by the allocation of the previous priority from the capacities of the edges
         # or set edge capacities if priority_idx = 1
@@ -930,10 +1115,14 @@ function allocate!(p::Parameters, allocation_model::AllocationModel, t::Float64)
         JuMP.optimize!(problem)
         @debug JuMP.solution_summary(problem)
         if JuMP.termination_status(problem) !== JuMP.OPTIMAL
-            error("Allocation coudn't find optimal solution.")
+            (; allocation_network_id) = allocation_model
+            priority = priorities[priority_index]
+            error(
+                "Allocation of subnetwork $allocation_network_id, priority $priority coudn't find optimal solution.",
+            )
         end
 
         # Assign the allocations to the users for this priority
-        assign_allocations!(allocation_model, p, t, priority_idx)
+        assign_allocations!(allocation_model, p, t, priority_idx; collect_demands)
     end
 end
diff --git a/core/src/bmi.jl b/core/src/bmi.jl
index 24c3b763a..c014dad0d 100644
--- a/core/src/bmi.jl
+++ b/core/src/bmi.jl
@@ -476,16 +476,30 @@ function discrete_control_affect!(
     return control_state_change
 end
 
-function get_allocation_model(
+function get_allocation_model(p::Parameters, allocation_network_id::Int)::AllocationModel
+    (; allocation) = p
+    (; allocation_network_ids, allocation_models) = allocation
+    idx = findsorted(allocation_network_ids, allocation_network_id)
+    if isnothing(idx)
+        error("Invalid allocation network ID $allocation_network_id.")
+    else
+        return allocation_models[idx]
+    end
+end
+
+function get_main_network_connections(
     p::Parameters,
     allocation_network_id::Int,
-)::Union{AllocationModel, Nothing}
-    for allocation_model in p.allocation_models
-        if allocation_model.allocation_network_id == allocation_network_id
-            return allocation_model
-        end
+)::Vector{Tuple{NodeID, NodeID}}
+    (; allocation) = p
+    (; allocation_network_ids, main_network_connections) = allocation
+    idx = findsorted(allocation_network_ids, allocation_network_id)
+    if isnothing(idx)
+        error("Invalid allocation network ID $allocation_network_id.")
+    else
+        return main_network_connections[idx]
     end
-    return nothing
+    return
 end
 
 """
@@ -589,7 +603,20 @@ end
 "Solve the allocation problem for all users and assign allocated abstractions to user nodes."
 function update_allocation!(integrator)::Nothing
     (; p, t) = integrator
-    for allocation_model in integrator.p.allocation_models
+    (; allocation) = p
+    (; allocation_models) = allocation
+
+    # If a main network is present, collect demands of subnetworks
+    if has_main_network(allocation)
+        for allocation_model in Iterators.drop(allocation_models, 1)
+            allocate!(p, allocation_model, t; collect_demands = true)
+        end
+    end
+
+    # Solve the allocation problems
+    # If a main network is present this is solved first,
+    # which provides allocation to the subnetworks
+    for allocation_model in allocation_models
         allocate!(p, allocation_model, t)
     end
 end
diff --git a/core/src/create.jl b/core/src/create.jl
index 331a944f3..7c5494117 100644
--- a/core/src/create.jl
+++ b/core/src/create.jl
@@ -196,16 +196,26 @@ end
 const nonconservative_nodetypes =
     Set{String}(["Basin", "LevelBoundary", "FlowBoundary", "Terminal", "User"])
 
-function generate_allocation_models!(p::Parameters, config::Config)::Nothing
-    (; graph, allocation_models) = p
+function initialize_allocation!(p::Parameters, config::Config)::Nothing
+    (; graph, allocation) = p
+    (; allocation_network_ids, allocation_models, main_network_connections) = allocation
+    allocation_network_ids_ = sort(collect(keys(graph[].node_ids)))
 
     errors = non_positive_allocation_network_id(graph)
-
     if errors
         error("Allocation network initialization failed.")
     end
 
-    for allocation_network_id in keys(graph[].node_ids)
+    for allocation_network_id in allocation_network_ids_
+        push!(allocation_network_ids, allocation_network_id)
+        push!(main_network_connections, Tuple{NodeID, NodeID}[])
+    end
+
+    if first(allocation_network_ids_) == 1
+        find_subnetwork_connections!(p)
+    end
+
+    for allocation_network_id in allocation_network_ids_
         push!(
             allocation_models,
             AllocationModel(config, allocation_network_id, p, config.allocation.timestep),
@@ -634,7 +644,7 @@ function User(db::DB, config::Config)::User
         error("Problems encountered when parsing User static and time node IDs.")
     end
 
-    # The highest priority number given, which corresponds to the least important demands
+    # All provided priorities
     priorities = sort(unique(union(static.priority, time.priority)))
 
     active = BitVector()
@@ -798,7 +808,13 @@ function Parameters(db::DB, config::Config)::Parameters
     n_states = length(get_ids(db, "Basin")) + length(get_ids(db, "PidControl"))
     chunk_sizes = get_chunk_sizes(config, n_states)
     graph = create_graph(db, config, chunk_sizes)
-    allocation_models = Vector{AllocationModel}()
+    allocation = Allocation(
+        Int[],
+        AllocationModel[],
+        Vector{Tuple{NodeID, NodeID}}[],
+        Dict{Tuple{NodeID, NodeID}, Float64}(),
+        Dict{Tuple{NodeID, NodeID}, Float64}(),
+    )
 
     if !valid_edges(graph)
         error("Invalid edge(s) found.")
@@ -835,7 +851,7 @@ function Parameters(db::DB, config::Config)::Parameters
     p = Parameters(
         config.starttime,
         graph,
-        allocation_models,
+        allocation,
         basin,
         linear_resistance,
         manning_resistance,
@@ -859,7 +875,7 @@ function Parameters(db::DB, config::Config)::Parameters
 
     # Allocation data structures
     if config.allocation.use_allocation
-        generate_allocation_models!(p, config)
+        initialize_allocation!(p, config)
     end
     return p
 end
diff --git a/core/src/solve.jl b/core/src/solve.jl
index 2087f2f35..d766a86d5 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -21,6 +21,22 @@ struct AllocationModel
     Δt_allocation::Float64
 end
 
+"""
+Object for all information about allocation
+allocation_network_ids: The unique sorted allocation network IDs
+allocation models: The allocation models for the main network and subnetworks corresponding to
+    allocation_network_ids
+main_network_connections: (from_id, to_id) from the main network to the subnetwork per subnetwork
+subnetwork_demands: The demand of an edge from the main network to a subnetwork
+"""
+struct Allocation
+    allocation_network_ids::Vector{Int}
+    allocation_models::Vector{AllocationModel}
+    main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}}
+    subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
+    subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
+end
+
 @enumx EdgeType flow control none
 
 """
@@ -471,7 +487,7 @@ struct Parameters{T, C1, C2}
         MetaGraphsNext.var"#11#13",
         Float64,
     }
-    allocation_models::Vector{AllocationModel}
+    allocation::Allocation
     basin::Basin{T, C1}
     linear_resistance::LinearResistance
     manning_resistance::ManningResistance
diff --git a/core/src/utils.jl b/core/src/utils.jl
index 27632f41b..f2f743da4 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -42,12 +42,19 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra
         db,
         "SELECT fid, from_node_id, to_node_id, edge_type, allocation_network_id FROM Edge ORDER BY fid",
     )
+    # Node IDs per subnetwork
     node_ids = Dict{Int, Set{NodeID}}()
+    # Allocation edges per subnetwork
     edge_ids = Dict{Int, Set{Tuple{NodeID, NodeID}}}()
+    # Source edges per subnetwork
     edges_source = Dict{Int, Set{EdgeMetadata}}()
+    # The number of flow edges
     flow_counter = 0
+    # Dictionary from flow edge to index in flow vector
     flow_dict = Dict{Tuple{NodeID, NodeID}, Int}()
+    # The number of nodes with vertical flow (interaction with outside of model)
     flow_vertical_counter = 0
+    # Dictionary from node ID to index in vertical flow vector
     flow_vertical_dict = Dict{NodeID, Int}()
     graph = MetaGraph(
         DiGraph();
@@ -1260,3 +1267,11 @@ function allocation_path_exists_in_graph(
     end
     return false
 end
+
+function has_main_network(allocation::Allocation)::Bool
+    return first(allocation.allocation_network_ids) == 1
+end
+
+function is_main_network(allocation_network_id::Int)::Bool
+    return allocation_network_id == 1
+end
diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl
index ff48730e6..5e1eaf0de 100644
--- a/core/test/allocation_test.jl
+++ b/core/test/allocation_test.jl
@@ -1,5 +1,4 @@
 @testitem "Allocation solve" begin
-    using PreallocationTools: get_tmp
     using Ribasim: NodeID
     import SQLite
     import JuMP
@@ -15,7 +14,7 @@
 
     graph = p.graph
     Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5) # Source flow
-    allocation_model = p.allocation_models[1]
+    allocation_model = p.allocation.allocation_models[1]
     Ribasim.allocate!(p, allocation_model, 0.0)
 
     F = allocation_model.problem[:F]
@@ -45,7 +44,7 @@ end
     config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_absolute")
     model = Ribasim.run(config)
     @test successful_retcode(model)
-    problem = model.integrator.p.allocation_models[1].problem
+    problem = model.integrator.p.allocation.allocation_models[1].problem
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.QuadExpr # Quadratic expression
     F = problem[:F]
@@ -61,7 +60,7 @@ end
     config = Ribasim.Config(toml_path; allocation_objective_type = "quadratic_relative")
     model = Ribasim.run(config)
     @test successful_retcode(model)
-    problem = model.integrator.p.allocation_models[1].problem
+    problem = model.integrator.p.allocation.allocation_models[1].problem
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.QuadExpr # Quadratic expression
     @test objective.aff.constant == 2.0
@@ -78,7 +77,7 @@ end
     config = Ribasim.Config(toml_path; allocation_objective_type = "linear_absolute")
     model = Ribasim.run(config)
     @test successful_retcode(model)
-    problem = model.integrator.p.allocation_models[1].problem
+    problem = model.integrator.p.allocation.allocation_models[1].problem
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.AffExpr # Affine expression
     @test :F_abs in keys(problem.obj_dict)
@@ -95,7 +94,7 @@ end
     config = Ribasim.Config(toml_path; allocation_objective_type = "linear_relative")
     model = Ribasim.run(config)
     @test successful_retcode(model)
-    problem = model.integrator.p.allocation_models[1].problem
+    problem = model.integrator.p.allocation.allocation_models[1].problem
     objective = JuMP.objective_function(problem)
     @test objective isa JuMP.AffExpr # Affine expression
     @test :F_abs in keys(problem.obj_dict)
@@ -121,7 +120,7 @@ end
         "../../generated_testmodels/fractional_flow_subnetwork/ribasim.toml",
     )
     model = Ribasim.BMI.initialize(Ribasim.Model, toml_path)
-    problem = model.integrator.p.allocation_models[1].problem
+    problem = model.integrator.p.allocation.allocation_models[1].problem
     F = problem[:F]
     @test JuMP.normalized_coefficient(
         problem[:fractional_flow][(NodeID(3), NodeID(5))],
@@ -155,7 +154,7 @@ end
     @test record_control.control_state == ["A", "B"]
 
     fractional_flow_constraints =
-        model.integrator.p.allocation_models[1].problem[:fractional_flow]
+        model.integrator.p.allocation.allocation_models[1].problem[:fractional_flow]
     @test JuMP.normalized_coefficient(
         problem[:fractional_flow][(NodeID(3), NodeID(5))],
         F[(NodeID(2), NodeID(3))],
@@ -165,3 +164,105 @@ end
         F[(NodeID(2), NodeID(3))],
     ) ≈ -0.25
 end
+
+@testitem "main allocation network initialization" begin
+    using SQLite
+    using Ribasim: NodeID
+
+    toml_path = normpath(
+        @__DIR__,
+        "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml",
+    )
+    @test ispath(toml_path)
+    cfg = Ribasim.Config(toml_path)
+    db_path = Ribasim.input_path(cfg, cfg.database)
+    db = SQLite.DB(db_path)
+    p = Ribasim.Parameters(db, cfg)
+    close(db)
+    (; allocation, graph) = p
+    (; main_network_connections, allocation_network_ids) = allocation
+    @test Ribasim.has_main_network(allocation)
+    @test Ribasim.is_main_network(first(allocation_network_ids))
+
+    # Connections from main network to subnetworks
+    @test isempty(main_network_connections[1])
+    @test only(main_network_connections[2]) == (NodeID(2), NodeID(11))
+    @test only(main_network_connections[3]) == (NodeID(6), NodeID(24))
+    @test only(main_network_connections[4]) == (NodeID(10), NodeID(38))
+
+    # main-sub connections are part of main network allocation graph
+    allocation_edges_main_network = graph[].edge_ids[1]
+    @test Tuple{NodeID, NodeID}[(2, 11), (6, 24), (10, 38)] ⊆ allocation_edges_main_network
+
+    # Subnetworks interpreted as users require variables and constraints to
+    # support absolute value expressions in the objective function
+    allocation_model_main_network = Ribasim.get_allocation_model(p, 1)
+    problem = allocation_model_main_network.problem
+    @test problem[:F_abs].axes[1] == NodeID[11, 24, 38]
+    @test problem[:abs_positive].axes[1] == NodeID[11, 24, 38]
+    @test problem[:abs_negative].axes[1] == NodeID[11, 24, 38]
+
+    # In each subnetwork, the connection from the main network to the subnetwork is
+    # interpreted as a source
+    @test Ribasim.get_allocation_model(p, 3).problem[:source].axes[1] ==
+          Tuple{NodeID, NodeID}[(2, 11)]
+    @test Ribasim.get_allocation_model(p, 5).problem[:source].axes[1] ==
+          Tuple{NodeID, NodeID}[(6, 24)]
+    @test Ribasim.get_allocation_model(p, 7).problem[:source].axes[1] ==
+          Tuple{NodeID, NodeID}[(10, 38)]
+end
+
+@testitem "allocation with main network optimization problem" begin
+    using SQLite
+    using Ribasim: NodeID
+    using JuMP
+
+    toml_path = normpath(
+        @__DIR__,
+        "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml",
+    )
+    @test ispath(toml_path)
+    cfg = Ribasim.Config(toml_path)
+    db_path = Ribasim.input_path(cfg, cfg.database)
+    db = SQLite.DB(db_path)
+    p = Ribasim.Parameters(db, cfg)
+    close(db)
+
+    (; allocation, user, graph) = p
+    (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation
+    t = 0.0
+
+    # Collecting demands
+    for allocation_model in allocation_models[2:end]
+        Ribasim.allocate!(p, allocation_model, t; collect_demands = true)
+    end
+
+    @test subnetwork_demands[(NodeID(2), NodeID(11))] ≈ [4.0, 4.0, 0.0]
+    @test subnetwork_demands[(NodeID(6), NodeID(24))] ≈ [0.001333333333, 0.0, 0.0]
+    @test subnetwork_demands[(NodeID(10), NodeID(38))] ≈ [0.001, 0.002, 0.002]
+
+    # Solving for the main network,
+    # containing subnetworks as users
+    allocation_model = allocation_models[1]
+    (; problem) = allocation_model
+    Ribasim.allocate!(p, allocation_model, t)
+
+    # Main network objective function
+    objective = JuMP.objective_function(problem)
+    objective_variables = keys(objective.terms)
+    F_abs = problem[:F_abs]
+    @test F_abs[NodeID(11)] ∈ objective_variables
+    @test F_abs[NodeID(24)] ∈ objective_variables
+    @test F_abs[NodeID(38)] ∈ objective_variables
+
+    # Running full allocation algorithm
+    Ribasim.set_flow!(graph, NodeID(1), NodeID(2), 4.5)
+    Ribasim.update_allocation!((; p, t))
+
+    @test subnetwork_allocateds[NodeID(2), NodeID(11)] ≈ [4.0, 0.49766666, 0.0]
+    @test subnetwork_allocateds[NodeID(6), NodeID(24)] ≈ [0.00133333333, 0.0, 0.0]
+    @test subnetwork_allocateds[NodeID(10), NodeID(38)] ≈ [0.001, 0.0, 0.0]
+
+    @test user.allocated[2] ≈ [4.0, 0.0, 0.0]
+    @test user.allocated[7] ≈ [0.001, 0.0, 0.0]
+end
diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd
index 14b16167e..e09d11765 100644
--- a/docs/core/allocation.qmd
+++ b/docs/core/allocation.qmd
@@ -238,15 +238,15 @@ using SQLite
 toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml")
 p = Ribasim.Model(toml_path).integrator.p
 
-allocation_model = p.allocation_models[1]
+allocation_model = p.allocation.allocation_models[1]
 t = 0.0
 priority_idx = 1
 
 Ribasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0)
 
-Ribasim.adjust_source_flows!(allocation_model, p, priority_idx)
+Ribasim.adjust_source_capacities!(allocation_model, p, priority_idx)
 Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx)
 Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx)
 
-println(p.allocation_models[1].problem)
+println(p.allocation.allocation_models[1].problem)
 ```
diff --git a/docs/python/test-models.qmd b/docs/python/test-models.qmd
index 6a60503cf..81891cd1d 100644
--- a/docs/python/test-models.qmd
+++ b/docs/python/test-models.qmd
@@ -2,7 +2,7 @@
 title: "Test models"
 ---
 
-Ribasim developers use the following models in its testbench and in order to test new features.
+Ribasim developers use the following models in their testbench and in order to test new features.
 
 ```{python}
 # | code-fold: true
diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py
index a4f9f552e..f10a2bc22 100644
--- a/python/ribasim/ribasim/geometry/node.py
+++ b/python/ribasim/ribasim/geometry/node.py
@@ -1,11 +1,13 @@
 from collections.abc import Sequence
 from typing import Any
 
+import geopandas as gpd
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import pandera as pa
 import shapely
+from matplotlib.patches import Patch
 from numpy.typing import NDArray
 from pandera.typing import Series
 from pandera.typing.geopandas import GeoSeries
@@ -130,6 +132,47 @@ def connectivity_from_geometry(
         to_id = node_index[edge_node_id[:, 1]].to_numpy()
         return from_id, to_id
 
+    def plot_allocation_networks(self, ax=None, zorder=None) -> Any:
+        if ax is None:
+            _, ax = plt.subplots()
+            ax.axis("off")
+
+        COLOR_SUBNETWORK = "black"
+        COLOR_MAIN_NETWORK = "blue"
+        ALPHA = 0.25
+
+        contains_main_network = False
+        contains_subnetworks = False
+
+        for allocation_subnetwork_id, df_subnetwork in self.df.groupby(
+            "allocation_network_id"
+        ):
+            if allocation_subnetwork_id is None:
+                continue
+            elif allocation_subnetwork_id == 1:
+                contains_main_network = True
+                color = COLOR_MAIN_NETWORK
+            else:
+                contains_subnetworks = True
+                color = COLOR_SUBNETWORK
+
+            hull = gpd.GeoDataFrame(
+                geometry=[df_subnetwork.geometry.unary_union.convex_hull]
+            )
+            hull.plot(ax=ax, color=color, alpha=ALPHA, zorder=zorder)
+
+        handles = []
+        labels = []
+
+        if contains_main_network:
+            handles.append(Patch(facecolor=COLOR_MAIN_NETWORK, alpha=ALPHA))
+            labels.append("Main network")
+        if contains_subnetworks:
+            handles.append(Patch(facecolor=COLOR_SUBNETWORK, alpha=ALPHA))
+            labels.append("Subnetwork")
+
+        return handles, labels
+
     def plot(self, ax=None, zorder=None) -> Any:
         """
         Plot the nodes. Each node type is given a separate marker.
diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py
index 5278b8efd..13b4509df 100644
--- a/python/ribasim/ribasim/model.py
+++ b/python/ribasim/ribasim/model.py
@@ -427,9 +427,9 @@ def plot_control_listen(self, ax):
                 label="Listen edge" if i == 0 else None,
             )
 
-    def plot(self, ax=None) -> Any:
+    def plot(self, ax=None, indicate_subnetworks: bool = True) -> Any:
         """
-        Plot the nodes and edges of the model.
+        Plot the nodes, edges and allocation networks of the model.
 
         Parameters
         ----------
@@ -443,11 +443,22 @@ def plot(self, ax=None) -> Any:
         if ax is None:
             _, ax = plt.subplots()
             ax.axis("off")
+
         self.network.edge.plot(ax=ax, zorder=2)
         self.plot_control_listen(ax)
         self.network.node.plot(ax=ax, zorder=3)
 
-        ax.legend(loc="lower left", bbox_to_anchor=(1, 0.5))
+        handles, labels = ax.get_legend_handles_labels()
+
+        if indicate_subnetworks:
+            (
+                handles_subnetworks,
+                labels_subnetworks,
+            ) = self.network.node.plot_allocation_networks(ax=ax, zorder=1)
+            handles += handles_subnetworks
+            labels += labels_subnetworks
+
+        ax.legend(handles, labels, loc="lower left", bbox_to_anchor=(1, 0.5))
 
         return ax
 
diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
index a73777f05..4616d2746 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py
@@ -9,6 +9,7 @@
     allocation_example_model,
     fractional_flow_subnetwork_model,
     looped_subnetwork_model,
+    main_network_with_subnetworks_model,
     minimal_subnetwork_model,
     subnetwork_model,
     user_model,
@@ -86,6 +87,7 @@
     "trivial_model",
     "two_basin_model",
     "user_model",
+    "main_network_with_subnetworks_model",
 ]
 
 # provide a mapping from model name to its constructor, so we can iterate over all models
diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
index 3be17b397..0e556823e 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
@@ -123,7 +123,9 @@ def user_model():
 
 
 def subnetwork_model():
-    """Create a user testmodel representing a subnetwork."""
+    """Create a user testmodel representing a subnetwork.
+    This model is merged into main_network_with_subnetworks_model.
+    """
 
     # Setup the nodes:
     xy = np.array(
@@ -164,7 +166,7 @@ def subnetwork_model():
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"type": node_type, "allocation_network_id": 1},
+            data={"type": node_type, "allocation_network_id": 2},
             index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -177,7 +179,7 @@ def subnetwork_model():
     )
     to_id = np.array([2, 3, 4, 10, 5, 6, 7, 8, 11, 12, 13, 9, 2, 6, 8], dtype=np.int64)
     allocation_network_id = len(from_id) * [None]
-    allocation_network_id[0] = 1
+    allocation_network_id[0] = 2
     lines = node.geometry_from_connectivity(from_id, to_id)
     edge = ribasim.Edge(
         df=gpd.GeoDataFrame(
@@ -283,7 +285,9 @@ def subnetwork_model():
 
 
 def looped_subnetwork_model():
-    """Create a user testmodel representing a subnetwork containing a loop in the topology."""
+    """Create a user testmodel representing a subnetwork containing a loop in the topology.
+    This model is merged into main_network_with_subnetworks_model.
+    """
     # Setup the nodes:
     xy = np.array(
         [
@@ -345,7 +349,7 @@ def looped_subnetwork_model():
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"type": node_type, "allocation_network_id": 1},
+            data={"type": node_type, "allocation_network_id": 2},
             index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -423,7 +427,7 @@ def looped_subnetwork_model():
     )
     lines = node.geometry_from_connectivity(from_id, to_id)
     allocation_network_id = len(from_id) * [None]
-    allocation_network_id[0] = 1
+    allocation_network_id[0] = 2
     edge = ribasim.Edge(
         df=gpd.GeoDataFrame(
             data={
@@ -556,7 +560,7 @@ def minimal_subnetwork_model():
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"type": node_type, "allocation_network_id": 1},
+            data={"type": node_type, "allocation_network_id": 2},
             index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -573,7 +577,7 @@ def minimal_subnetwork_model():
         dtype=np.int64,
     )
     allocation_network_id = len(from_id) * [None]
-    allocation_network_id[0] = 1
+    allocation_network_id[0] = 2
     lines = node.geometry_from_connectivity(from_id, to_id)
     edge = ribasim.Edge(
         df=gpd.GeoDataFrame(
@@ -677,7 +681,9 @@ def minimal_subnetwork_model():
 
 
 def fractional_flow_subnetwork_model():
-    """Create a small subnetwork that contains fractional flow nodes."""
+    """Create a small subnetwork that contains fractional flow nodes.
+    This model is merged into main_network_with_subnetworks_model.
+    """
 
     xy = np.array(
         [
@@ -711,7 +717,7 @@ def fractional_flow_subnetwork_model():
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"type": node_type, "allocation_network_id": 1},
+            data={"type": node_type, "allocation_network_id": 2},
             index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -728,7 +734,7 @@ def fractional_flow_subnetwork_model():
         dtype=np.int64,
     )
     allocation_network_id = len(from_id) * [None]
-    allocation_network_id[0] = 1
+    allocation_network_id[0] = 2
     lines = node.geometry_from_connectivity(from_id, to_id)
     edge = ribasim.Edge(
         df=gpd.GeoDataFrame(
@@ -903,10 +909,10 @@ def allocation_example_model():
         "User",
     ]
 
-    # All nodes belong to allocation network id 1
+    # All nodes belong to allocation network id 2
     node = ribasim.Node(
         df=gpd.GeoDataFrame(
-            data={"type": node_type, "allocation_network_id": 1},
+            data={"type": node_type, "allocation_network_id": 2},
             index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -924,7 +930,7 @@ def allocation_example_model():
     # Denote the first edge, 1 => 2, as a source edge for
     # allocation network 1
     allocation_network_id = len(from_id) * [None]
-    allocation_network_id[0] = 1
+    allocation_network_id[0] = 2
     lines = node.geometry_from_connectivity(from_id, to_id)
     edge = ribasim.Edge(
         df=gpd.GeoDataFrame(
@@ -1072,3 +1078,660 @@ def allocation_example_model():
     )
 
     return model
+
+
+def main_network_with_subnetworks_model():
+    # Set ip the nodes:
+    xy = np.array(
+        [
+            (0.0, -1.0),
+            (3.0, 1.0),
+            (6.0, -1.0),
+            (9.0, 1.0),
+            (12.0, -1.0),
+            (15.0, 1.0),
+            (18.0, -1.0),
+            (21.0, 1.0),
+            (24.0, -1.0),
+            (27.0, 1.0),
+            (3.0, 4.0),
+            (2.0, 4.0),
+            (1.0, 4.0),
+            (0.0, 4.0),
+            (2.0, 5.0),
+            (2.0, 6.0),
+            (1.0, 6.0),
+            (0.0, 6.0),
+            (2.0, 8.0),
+            (2.0, 3.0),
+            (3.0, 6.0),
+            (0.0, 7.0),
+            (2.0, 7.0),
+            (14.0, 3.0),
+            (14.0, 4.0),
+            (14.0, 5.0),
+            (13.0, 6.0),
+            (12.0, 7.0),
+            (11.0, 8.0),
+            (15.0, 6.0),
+            (16.0, 7.0),
+            (17.0, 8.0),
+            (13.0, 5.0),
+            (26.0, 3.0),
+            (26.0, 4.0),
+            (25.0, 4.0),
+            (24.0, 4.0),
+            (28.0, 4.0),
+            (26.0, 5.0),
+            (24.0, 6.0),
+            (25.0, 6.0),
+            (26.0, 6.0),
+            (27.0, 6.0),
+            (28.0, 6.0),
+            (24.0, 7.0),
+            (26.0, 7.0),
+            (28.0, 7.0),
+            (26.0, 8.0),
+            (27.0, 8.0),
+            (28.0, 8.0),
+            (25.0, 9.0),
+            (26.0, 9.0),
+            (28.0, 9.0),
+            (26.0, 10.0),
+            (26.0, 11.0),
+            (26.0, 12.0),
+            (29.0, 6.0),
+        ]
+    )
+    node_xy = gpd.points_from_xy(x=xy[:, 0], y=xy[:, 1])
+
+    node_type = [
+        "FlowBoundary",
+        "Basin",
+        "LinearResistance",
+        "Basin",
+        "LinearResistance",
+        "Basin",
+        "LinearResistance",
+        "Basin",
+        "LinearResistance",
+        "Basin",
+        "Pump",
+        "Basin",
+        "Outlet",
+        "Terminal",
+        "Pump",
+        "Basin",
+        "Outlet",
+        "Basin",
+        "Terminal",
+        "User",
+        "User",
+        "User",
+        "Outlet",
+        "Pump",
+        "Basin",
+        "TabulatedRatingCurve",
+        "FractionalFlow",
+        "Basin",
+        "User",
+        "FractionalFlow",
+        "Basin",
+        "User",
+        "DiscreteControl",
+        "User",
+        "Basin",
+        "Outlet",
+        "Terminal",
+        "Pump",
+        "Pump",
+        "Basin",
+        "Outlet",
+        "Basin",
+        "Outlet",
+        "Basin",
+        "User",
+        "TabulatedRatingCurve",
+        "TabulatedRatingCurve",
+        "Basin",
+        "Pump",
+        "Basin",
+        "User",
+        "TabulatedRatingCurve",
+        "User",
+        "Basin",
+        "Outlet",
+        "Terminal",
+        "User",
+    ]
+
+    allocation_network_id = np.ones(57, dtype=int)
+    allocation_network_id[10:23] = 3
+    allocation_network_id[23:33] = 5
+    allocation_network_id[33:] = 7
+
+    # Make sure the feature id starts at 1: explicitly give an index.
+    node = ribasim.Node(
+        df=gpd.GeoDataFrame(
+            data={
+                "type": node_type,
+                "allocation_network_id": allocation_network_id,
+            },
+            index=pd.Index(np.arange(len(xy)) + 1, name="fid"),
+            geometry=node_xy,
+            crs="EPSG:28992",
+        )
+    )
+
+    # Setup the edges:
+    from_id = np.array(
+        [
+            1,
+            2,
+            3,
+            4,
+            5,
+            6,
+            7,
+            8,
+            9,
+            11,
+            12,
+            13,
+            12,
+            12,
+            15,
+            16,
+            17,
+            16,
+            18,
+            16,
+            23,
+            20,
+            21,
+            22,
+            24,
+            25,
+            26,
+            27,
+            28,
+            29,
+            26,
+            30,
+            31,
+            32,
+            33,
+            33,
+            38,
+            35,
+            36,
+            35,
+            35,
+            39,
+            42,
+            41,
+            40,
+            42,
+            46,
+            48,
+            49,
+            50,
+            48,
+            52,
+            48,
+            51,
+            54,
+            55,
+            42,
+            43,
+            44,
+            47,
+            34,
+            45,
+            53,
+            44,
+            57,
+            2,
+            6,
+            10,
+        ],
+        dtype=np.int64,
+    )
+    to_id = np.array(
+        [
+            2,
+            3,
+            4,
+            5,
+            6,
+            7,
+            8,
+            9,
+            10,
+            12,
+            13,
+            14,
+            20,
+            15,
+            16,
+            17,
+            18,
+            21,
+            22,
+            23,
+            19,
+            12,
+            16,
+            18,
+            25,
+            26,
+            27,
+            28,
+            29,
+            28,
+            30,
+            31,
+            32,
+            31,
+            27,
+            30,
+            35,
+            36,
+            37,
+            34,
+            39,
+            42,
+            41,
+            40,
+            45,
+            46,
+            48,
+            49,
+            50,
+            53,
+            52,
+            54,
+            51,
+            54,
+            55,
+            56,
+            43,
+            44,
+            47,
+            50,
+            35,
+            40,
+            50,
+            57,
+            44,
+            11,
+            24,
+            38,
+        ],
+        dtype=np.int64,
+    )
+
+    edge_type = 68 * ["flow"]
+    edge_type[34] = "control"
+    edge_type[35] = "control"
+    allocation_network_id = 68 * [None]
+    allocation_network_id[0] = 1
+    allocation_network_id[65] = 3
+    allocation_network_id[66] = 5
+    allocation_network_id[67] = 7
+
+    lines = node.geometry_from_connectivity(from_id.tolist(), to_id.tolist())
+    edge = ribasim.Edge(
+        df=gpd.GeoDataFrame(
+            data={
+                "from_node_id": from_id,
+                "to_node_id": to_id,
+                "edge_type": edge_type,
+                "allocation_network_id": allocation_network_id,
+            },
+            geometry=lines,
+            crs="EPSG:28992",
+        )
+    )
+
+    # Setup the basins:
+    profile = pd.DataFrame(
+        data={
+            "node_id": [
+                2,
+                2,
+                4,
+                4,
+                6,
+                6,
+                8,
+                8,
+                10,
+                10,
+                12,
+                12,
+                16,
+                16,
+                18,
+                18,
+                25,
+                25,
+                28,
+                28,
+                31,
+                31,
+                35,
+                35,
+                40,
+                40,
+                42,
+                42,
+                44,
+                44,
+                48,
+                48,
+                50,
+                50,
+                54,
+                54,
+            ],
+            "area": [
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                100000.0,
+                100000.0,
+                100000.0,
+                100000.0,
+                100000.0,
+                100000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+                1000.0,
+            ],
+            "level": [
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+                0.0,
+                1.0,
+            ],
+        }
+    )
+
+    static = pd.DataFrame(
+        data={
+            "node_id": [
+                2,
+                4,
+                6,
+                8,
+                10,
+                12,
+                16,
+                18,
+                25,
+                28,
+                31,
+                35,
+                40,
+                44,
+                42,
+                48,
+                50,
+                54,
+            ],
+            "drainage": 0.0,
+            "potential_evaporation": 0.0,
+            "infiltration": 0.0,
+            "precipitation": 0.0,
+            "urban_runoff": 0.0,
+        }
+    )
+
+    state = pd.DataFrame(
+        data={
+            "node_id": [
+                2,
+                4,
+                6,
+                8,
+                10,
+                12,
+                16,
+                18,
+                25,
+                28,
+                31,
+                35,
+                40,
+                42,
+                44,
+                48,
+                50,
+                54,
+            ],
+            "level": [
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                10.0,
+                10.0,
+                10.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+                1.0,
+            ],
+        }
+    )
+
+    basin = ribasim.Basin(
+        profile=profile,
+        static=static,
+        state=state,
+    )
+
+    # Setup the discrete control:
+    condition = pd.DataFrame(
+        data={
+            "node_id": [33],
+            "listen_feature_id": [25],
+            "variable": ["level"],
+            "greater_than": [0.003],
+        }
+    )
+
+    logic = pd.DataFrame(
+        data={
+            "node_id": [33, 33],
+            "truth_state": ["F", "T"],
+            "control_state": ["A", "B"],
+        }
+    )
+
+    discrete_control = ribasim.DiscreteControl(condition=condition, logic=logic)
+
+    # Setup flow boundary
+    flow_boundary = ribasim.FlowBoundary(
+        static=pd.DataFrame(data={"node_id": [1], "flow_rate": [1.0]})
+    )
+
+    # Setup fractional flow
+    fractional_flow = ribasim.FractionalFlow(
+        static=pd.DataFrame(
+            data={
+                "node_id": [27, 30, 27, 30],
+                "fraction": [0.25, 0.75, 0.75, 0.25],
+                "control_state": ["A", "A", "B", "B"],
+            }
+        )
+    )
+
+    # Setup linear resistance
+    linear_resistance = ribasim.LinearResistance(
+        static=pd.DataFrame(data={"node_id": [3, 5, 7, 9], "resistance": 0.001})
+    )
+
+    # Setup outlet
+    outlet = ribasim.Outlet(
+        static=pd.DataFrame(
+            data={
+                "node_id": [13, 17, 23, 36, 41, 43, 55],
+                "flow_rate": [3.0, 3.0, 3.0, 0.003, 0.003, 0.003, 0.003],
+                "max_flow_rate": 3.0,
+            }
+        )
+    )
+
+    # Setup pump
+    pump = ribasim.Pump(
+        static=pd.DataFrame(
+            data={
+                "node_id": [15, 39, 49, 11, 24, 38],
+                "flow_rate": [4.0e00, 4.0e-03, 4.0e-03, 1.0e-03, 1.0e-03, 1.0e-03],
+                "max_flow_rate": [4.0, 0.004, 0.004, 1.0, 1.0, 1.0],
+            }
+        )
+    )
+
+    # Setup tabulated rating curve
+    rating_curve = ribasim.TabulatedRatingCurve(
+        static=pd.DataFrame(
+            data={
+                "node_id": [26, 26, 46, 46, 47, 47, 52, 52],
+                "level": [0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0],
+                "flow_rate": [
+                    0.0e00,
+                    1.0e-04,
+                    0.0e00,
+                    2.0e00,
+                    0.0e00,
+                    2.0e00,
+                    0.0e00,
+                    2.0e00,
+                ],
+            }
+        )
+    )
+
+    # Setup terminal node
+    terminal = ribasim.Terminal(static=pd.DataFrame(data={"node_id": [14, 19, 37, 56]}))
+
+    # Setup the user
+    user = ribasim.User(
+        static=pd.DataFrame(
+            data={
+                "node_id": [20, 21, 22, 29, 34, 45, 51, 53, 57],
+                "demand": [
+                    4.0e00,
+                    5.0e00,
+                    3.0e00,
+                    1.0e-03,
+                    1.0e-03,
+                    1.0e-03,
+                    1.0e-03,
+                    1.0e-03,
+                    1.0e-03,
+                ],
+                "return_factor": 0.9,
+                "min_level": 0.9,
+                "priority": [2, 1, 2, 1, 2, 1, 3, 3, 2],
+            }
+        ),
+        time=pd.DataFrame(
+            data={
+                "node_id": [32, 32],
+                "time": ["2020-01-01 00:00:00", "2021-01-01 00:00:00"],
+                "demand": [0.001, 0.002],
+                "return_factor": 0.9,
+                "min_level": 0.9,
+                "priority": 1,
+            }
+        ),
+    )
+
+    # Setup allocation:
+    allocation = ribasim.Allocation(use_allocation=True, timestep=86400)
+
+    model = ribasim.Model(
+        network=ribasim.Network(node=node, edge=edge),
+        basin=basin,
+        discrete_control=discrete_control,
+        flow_boundary=flow_boundary,
+        fractional_flow=fractional_flow,
+        linear_resistance=linear_resistance,
+        outlet=outlet,
+        pump=pump,
+        terminal=terminal,
+        user=user,
+        tabulated_rating_curve=rating_curve,
+        allocation=allocation,
+        starttime="2020-01-01 00:00:00",
+        endtime="2021-01-01 00:00:00",
+    )
+
+    return model