diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl
index 90e033229..d80660861 100644
--- a/core/src/allocation_init.jl
+++ b/core/src/allocation_init.jl
@@ -95,7 +95,7 @@ function get_subnetwork_capacity(
     return capacity
 end
 
-const allocation_source_nodetypes =
+const boundary_source_nodetypes =
     Set{NodeType.T}([NodeType.LevelBoundary, NodeType.FlowBoundary])
 
 """
@@ -138,10 +138,6 @@ function get_capacity(
     capacity = get_subnetwork_capacity(p, subnetwork_id)
     add_subnetwork_connections!(capacity, p, subnetwork_id)
 
-    if !valid_sources(p, capacity, subnetwork_id)
-        error("Errors in sources in allocation network.")
-    end
-
     return capacity
 end
 
@@ -272,36 +268,55 @@ function add_constraints_user_source!(
 end
 
 """
-Add the source constraints to the allocation problem.
+Add the boundary source constraints to the allocation problem.
 The actual threshold values will be set before each allocation solve.
 The constraint indices are (edge_source_id, edge_dst_id).
 
 Constraint:
-flow over source edge <= source flow in subnetwork
+flow over source edge <= source flow in physical layer
 """
-function add_constraints_source!(
+function add_constraints_boundary_source!(
     problem::JuMP.Model,
     p::Parameters,
     subnetwork_id::Int32,
 )::Nothing
-    (; graph) = p
-    edges_source = Tuple{NodeID, NodeID}[]
+    # Source edges (without the basins)
+    edges_source =
+        [edge for edge in source_edges_subnetwork(p, subnetwork_id) if edge[1] != edge[2]]
     F = problem[:F]
 
-    # Find the edges in the whole model which are a source for
-    # this subnetwork
-    for edge_metadata in values(graph.edge_data)
-        (; edge) = edge_metadata
-        if graph[edge...].subnetwork_id_source == subnetwork_id
-            push!(edges_source, edge)
-        end
-    end
+    problem[:source_boundary] = JuMP.@constraint(
+        problem,
+        [edge_id = edges_source],
+        F[edge_id] <= 0.0,
+        base_name = "source_boundary"
+    )
+    return nothing
+end
 
-    problem[:source] = JuMP.@constraint(
+"""
+Add main network source constraints to the allocation problem.
+The actual threshold values will be set before each allocation solve.
+The constraint indices are (edge_source_id, edge_dst_id).
+
+Constraint:
+flow over main network to subnetwork connection edge <= either 0 or allocated amount from the main network
+"""
+function add_constraints_main_network_source!(
+    problem::JuMP.Model,
+    p::Parameters,
+    subnetwork_id::Int32,
+)::Nothing
+    F = problem[:F]
+    (; main_network_connections, subnetwork_ids) = p.allocation
+    subnetwork_id = searchsortedfirst(subnetwork_ids, subnetwork_id)
+    edges_source = main_network_connections[subnetwork_id]
+
+    problem[:source_main_network] = JuMP.@constraint(
         problem,
         [edge_id = edges_source],
         F[edge_id] <= 0.0,
-        base_name = "source"
+        base_name = "source_main_network"
     )
     return nothing
 end
@@ -451,9 +466,9 @@ function allocation_problem(
 
     # Add constraints to problem
     add_constraints_conservation_node!(problem, p, subnetwork_id)
-
     add_constraints_capacity!(problem, capacity, p, subnetwork_id)
-    add_constraints_source!(problem, p, subnetwork_id)
+    add_constraints_boundary_source!(problem, p, subnetwork_id)
+    add_constraints_main_network_source!(problem, p, subnetwork_id)
     add_constraints_user_source!(problem, p, subnetwork_id)
     add_constraints_basin_flow!(problem)
     add_constraints_buffer!(problem)
@@ -484,12 +499,12 @@ function get_sources_in_order(
         sources[edge] = AllocationSource(; edge, type = AllocationSourceType.user_return)
     end
 
-    # Source edges (within subnetwork)
-    for edge in
-        sort(only(problem[:source].axes); by = edge -> (edge[1].value, edge[2].value))
-        if graph[edge[1]].subnetwork_id == graph[edge[2]].subnetwork_id
-            sources[edge] = AllocationSource(; edge, type = AllocationSourceType.edge)
-        end
+    # Boundary node sources
+    for edge in sort(
+        only(problem[:source_boundary].axes);
+        by = edge -> (edge[1].value, edge[2].value),
+    )
+        sources[edge] = AllocationSource(; edge, type = AllocationSourceType.boundary_node)
     end
 
     # Basins with level demand
diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl
index 9a3f0d4e0..6e2404862 100644
--- a/core/src/allocation_optim.jl
+++ b/core/src/allocation_optim.jl
@@ -127,8 +127,7 @@ function assign_allocations!(
         # 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 optimization_type == OptimizationType.collect_demands
-            if graph[edge...].subnetwork_id_source == subnetwork_id &&
-               edge ∈ main_network_source_edges
+            if edge in main_network_source_edges
                 allocated = flow[edge]
                 subnetwork_demands[edge][priority_idx] += allocated
             end
@@ -203,22 +202,13 @@ function set_initial_capacities_source!(
     p::Parameters,
 )::Nothing
     (; sources) = allocation_model
-    (; graph, allocation) = p
-    (; mean_input_flows) = allocation
     (; subnetwork_id) = allocation_model
-    main_network_source_edges = get_main_network_connections(p, subnetwork_id)
 
-    for edge_metadata in values(graph.edge_data)
-        (; edge) = edge_metadata
-        if graph[edge...].subnetwork_id_source == subnetwork_id
-            # If it is a source edge for this allocation problem
-            if edge ∉ main_network_source_edges
-                # Reset the source to the averaged flow over the last allocation period
-                source = sources[edge]
-                @assert source.type == AllocationSourceType.edge
-                source.capacity = mean_input_flows[edge][]
-            end
-        end
+    mean_input_flows_subnetwork_ = mean_input_flows_subnetwork(p, subnetwork_id)
+
+    for edge in keys(mean_input_flows_subnetwork_)
+        source = sources[edge]
+        source.capacity = mean_input_flows_subnetwork_[edge]
     end
     return nothing
 end
@@ -231,7 +221,7 @@ function reduce_source_capacity!(problem::JuMP.Model, source::AllocationSource):
 
     used_capacity =
         if source.type in (
-            AllocationSourceType.edge,
+            AllocationSourceType.boundary_node,
             AllocationSourceType.main_to_sub,
             AllocationSourceType.user_return,
         )
@@ -343,11 +333,10 @@ function get_basin_data(
     u::ComponentVector,
     node_id::NodeID,
 )
-    (; graph, allocation, basin) = p
-    (; Δt_allocation) = allocation_model
-    (; mean_input_flows) = allocation
+    (; graph, basin) = p
+    (; Δt_allocation, subnetwork_id) = allocation_model
     @assert node_id.type == NodeType.Basin
-    influx = mean_input_flows[(node_id, node_id)][]
+    influx = mean_input_flows_subnetwork(p, subnetwork_id)[(node_id, node_id)]
     storage_basin = basin.current_properties.current_storage[parent(u)][node_id.idx]
     control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control)
     if isempty(control_inneighbors)
@@ -841,9 +830,10 @@ function set_source_capacity!(
     optimization_type::OptimizationType.T,
 )::Nothing
     (; problem, sources) = allocation_model
-    constraints_source_edge = problem[:source]
-    constraints_source_basin = problem[:basin_outflow]
+    constraints_source_boundary = problem[:source_boundary]
     constraints_source_user_out = problem[:source_user]
+    constraints_source_main_network = problem[:source_main_network]
+    constraints_source_basin = problem[:basin_outflow]
     constraints_source_buffer = problem[:flow_buffer_outflow]
 
     for source in values(sources)
@@ -860,16 +850,17 @@ function set_source_capacity!(
             0.0
         end
 
-        constraint =
-            if source.type in (AllocationSourceType.edge, AllocationSourceType.main_to_sub)
-                constraints_source_edge[edge]
-            elseif source.type == AllocationSourceType.basin
-                constraints_source_basin[edge[1]]
-            elseif source.type == AllocationSourceType.user_return
-                constraints_source_user_out[edge[1]]
-            elseif source.type == AllocationSourceType.buffer
-                constraints_source_buffer[edge[1]]
-            end
+        constraint = if source.type == AllocationSourceType.boundary_node
+            constraints_source_boundary[edge]
+        elseif source.type == AllocationSourceType.main_to_sub
+            constraints_source_main_network[edge]
+        elseif source.type == AllocationSourceType.basin
+            constraints_source_basin[edge[1]]
+        elseif source.type == AllocationSourceType.user_return
+            constraints_source_user_out[edge[1]]
+        elseif source.type == AllocationSourceType.buffer
+            constraints_source_buffer[edge[1]]
+        end
 
         JuMP.set_normalized_rhs(constraint, capacity_effective)
     end
@@ -1062,7 +1053,8 @@ function empty_sources!(allocation_model::AllocationModel, allocation::Allocatio
     (; problem) = allocation_model
     (; subnetwork_demands) = allocation
 
-    for constraint_set_name in [:source, :source_user, :basin_outflow, :flow_buffer_outflow]
+    for constraint_set_name in
+        [:source_boundary, :source_user, :basin_outflow, :flow_buffer_outflow]
         constraint_set = problem[constraint_set_name]
         for key in only(constraint_set.axes)
             # Do not set the capacity to 0.0 if the edge
diff --git a/core/src/callback.jl b/core/src/callback.jl
index 213e97284..467e71060 100644
--- a/core/src/callback.jl
+++ b/core/src/callback.jl
@@ -153,8 +153,11 @@ function update_cumulative_flows!(u, t, integrator)::Nothing
     end
 
     # Update realized flows for allocation input
-    for edge in keys(allocation.mean_input_flows)
-        allocation.mean_input_flows[edge] += flow_update_on_edge(integrator, edge)
+    for subnetwork_id in allocation.subnetwork_ids
+        mean_input_flows_subnetwork_ = mean_input_flows_subnetwork(p, subnetwork_id)
+        for edge in keys(mean_input_flows_subnetwork_)
+            mean_input_flows_subnetwork_[edge] += flow_update_on_edge(integrator, edge)
+        end
     end
 
     # Update realized flows for allocation output
@@ -773,8 +776,10 @@ function update_allocation!(integrator)::Nothing
 
     # Divide by the allocation Δt to get the mean input flows from the cumulative flows
     (; Δt_allocation) = allocation_models[1]
-    for edge in keys(mean_input_flows)
-        mean_input_flows[edge] /= Δt_allocation
+    for mean_input_flows_subnetwork in values(mean_input_flows)
+        for edge in keys(mean_input_flows_subnetwork)
+            mean_input_flows_subnetwork[edge] /= Δt_allocation
+        end
     end
 
     # Divide by the allocation Δt to get the mean realized flows from the cumulative flows
@@ -797,11 +802,14 @@ function update_allocation!(integrator)::Nothing
     end
 
     # Reset the mean flows
-    for mean_flows in (mean_input_flows, mean_realized_flows)
+    for mean_flows in mean_input_flows
         for edge in keys(mean_flows)
             mean_flows[edge] = 0.0
         end
     end
+    for edge in keys(mean_realized_flows)
+        mean_realized_flows[edge] = 0.0
+    end
 end
 
 "Load updates from 'TabulatedRatingCurve / time' into the parameters"
diff --git a/core/src/graph.jl b/core/src/graph.jl
index d7e32d6a8..e2a3273e0 100644
--- a/core/src/graph.jl
+++ b/core/src/graph.jl
@@ -19,8 +19,7 @@ function create_graph(db::DB, config::Config)::MetaGraph
             FromNode.node_type AS from_node_type,
             ToNode.node_id AS to_node_id,
             ToNode.node_type AS to_node_type,
-            Edge.edge_type,
-            Edge.subnetwork_id
+            Edge.edge_type
         FROM Edge
         LEFT JOIN Node AS FromNode ON FromNode.node_id = Edge.from_node_id
         LEFT JOIN Node AS ToNode ON ToNode.node_id = Edge.to_node_id
@@ -28,8 +27,7 @@ function create_graph(db::DB, config::Config)::MetaGraph
     )
     # Node IDs per subnetwork
     node_ids = Dict{Int32, Set{NodeID}}()
-    # Source edges per subnetwork
-    edges_source = Dict{Int32, Set{EdgeMetadata}}()
+
     # The metadata of the flow edges in the order in which they are in the input
     # and will be in the output
     flow_edges = EdgeMetadata[]
@@ -57,15 +55,8 @@ function create_graph(db::DB, config::Config)::MetaGraph
     end
 
     errors = false
-    for (;
-        edge_id,
-        from_node_type,
-        from_node_id,
-        to_node_type,
-        to_node_id,
-        edge_type,
-        subnetwork_id,
-    ) in edge_rows
+    for (; edge_id, from_node_type, from_node_id, to_node_type, to_node_id, edge_type) in
+        edge_rows
         try
             # hasfield does not work
             edge_type = getfield(EdgeType, Symbol(edge_type))
@@ -74,15 +65,8 @@ function create_graph(db::DB, config::Config)::MetaGraph
         end
         id_src = NodeID(from_node_type, from_node_id, db)
         id_dst = NodeID(to_node_type, to_node_id, db)
-        if ismissing(subnetwork_id)
-            subnetwork_id = 0
-        end
-        edge_metadata = EdgeMetadata(;
-            id = edge_id,
-            type = edge_type,
-            subnetwork_id_source = subnetwork_id,
-            edge = (id_src, id_dst),
-        )
+        edge_metadata =
+            EdgeMetadata(; id = edge_id, type = edge_type, edge = (id_src, id_dst))
         if edge_type == EdgeType.flow
             push!(flow_edges, edge_metadata)
         end
@@ -91,12 +75,6 @@ function create_graph(db::DB, config::Config)::MetaGraph
             @error "Duplicate edge" id_src id_dst
         end
         graph[id_src, id_dst] = edge_metadata
-        if subnetwork_id != 0
-            if !haskey(edges_source, subnetwork_id)
-                edges_source[subnetwork_id] = Set{EdgeMetadata}()
-            end
-            push!(edges_source[subnetwork_id], edge_metadata)
-        end
     end
     if errors
         error("Invalid edges found")
@@ -106,7 +84,7 @@ function create_graph(db::DB, config::Config)::MetaGraph
         error("Incomplete connectivity in subnetwork")
     end
 
-    graph_data = (; node_ids, edges_source, flow_edges, config.solver.saveat)
+    graph_data = (; node_ids, flow_edges, config.solver.saveat)
     @reset graph.graph_data = graph_data
 
     return graph
diff --git a/core/src/model.jl b/core/src/model.jl
index e0047ca0d..2197d72b7 100644
--- a/core/src/model.jl
+++ b/core/src/model.jl
@@ -47,6 +47,7 @@ function Model(config::Config)::Model
     # so we can directly close it again.
     db = SQLite.DB(db_path)
 
+    database_warning(db)
     if !valid_nodes(db)
         error("Invalid nodes found.")
     end
diff --git a/core/src/parameter.jl b/core/src/parameter.jl
index 4c56611be..03b909799 100644
--- a/core/src/parameter.jl
+++ b/core/src/parameter.jl
@@ -128,7 +128,7 @@ The caches are always initialized with zeros
 """
 cache(len::Int)::Cache = LazyBufferCache(Returns(len); initializer! = set_zero!)
 
-@enumx AllocationSourceType edge basin main_to_sub user_return buffer
+@enumx AllocationSourceType boundary_node basin main_to_sub user_return buffer
 
 """
 Data structure for a single source within an allocation subnetwork.
@@ -180,20 +180,21 @@ main_network_connections: (from_id, to_id) from the main network to the subnetwo
 priorities: All used priority values.
 subnetwork_demands: The demand of an edge from the main network to a subnetwork
 subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork
-mean_input_flows: Flows averaged over Δt_allocation over edges that are allocation sources
+mean_input_flows: Per subnetwork, flows averaged over Δt_allocation over edges that are allocation sources
 mean_realized_flows: Flows averaged over Δt_allocation over edges that realize a demand
 record_demand: A record of demands and allocated flows for nodes that have these
 record_flow: A record of all flows computed by allocation optimization, eventually saved to
     output file
 """
 @kwdef struct Allocation
-    subnetwork_ids::Vector{Int32} = []
-    allocation_models::Vector{AllocationModel} = []
-    main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} = []
+    subnetwork_ids::Vector{Int32} = Int32[]
+    allocation_models::Vector{AllocationModel} = AllocationModel[]
+    main_network_connections::Vector{Vector{Tuple{NodeID, NodeID}}} =
+        Vector{Tuple{NodeID, NodeID}}[]
     priorities::Vector{Int32}
     subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} = Dict()
     subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}} = Dict()
-    mean_input_flows::Dict{Tuple{NodeID, NodeID}, Float64}
+    mean_input_flows::Vector{Dict{Tuple{NodeID, NodeID}, Float64}}
     mean_realized_flows::Dict{Tuple{NodeID, NodeID}, Float64}
     record_demand::@NamedTuple{
         time::Vector{Float64},
@@ -255,14 +256,11 @@ end
 Type for storing metadata of edges in the graph:
 id: ID of the edge (only used for labeling flow output)
 type: type of the edge
-subnetwork_id_source: ID of subnetwork where this edge is a source
-  (0 if not a source)
 edge: (from node ID, to node ID)
 """
 @kwdef struct EdgeMetadata
     id::Int32
     type::EdgeType.T
-    subnetwork_id_source::Int32
     edge::Tuple{NodeID, NodeID}
 end
 
@@ -902,7 +900,6 @@ const ModelGraph = MetaGraph{
     EdgeMetadata,
     @NamedTuple{
         node_ids::Dict{Int32, Set{NodeID}},
-        edges_source::Dict{Int32, Set{EdgeMetadata}},
         flow_edges::Vector{EdgeMetadata},
         saveat::Float64,
     },
diff --git a/core/src/read.jl b/core/src/read.jl
index edcce108c..e99d207e0 100644
--- a/core/src/read.jl
+++ b/core/src/read.jl
@@ -1268,40 +1268,56 @@ function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid
 end
 
 function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation
-    mean_input_flows = Dict{Tuple{NodeID, NodeID}, Float64}()
+    mean_input_flows = Dict{Tuple{NodeID, NodeID}, Float64}[]
+    mean_realized_flows = Dict{Tuple{NodeID, NodeID}, Float64}()
+    subnetwork_ids = sort(collect(keys(graph[].node_ids)))
 
-    # Find edges which serve as sources in allocation
-    for edge_metadata in values(graph.edge_data)
-        (; subnetwork_id_source, edge) = edge_metadata
-        if subnetwork_id_source != 0
-            mean_input_flows[edge] = 0.0
+    if config.allocation.use_allocation
+        for subnetwork_id in subnetwork_ids
+            push!(mean_input_flows, Dict{Tuple{NodeID, NodeID}, Float64}())
         end
-    end
 
-    # Find basins with a level demand
-    for node_id in values(graph.vertex_labels)
-        if has_external_demand(graph, node_id, :level_demand)[1]
-            mean_input_flows[(node_id, node_id)] = 0.0
+        # Find edges which serve as sources in allocation
+        for edge_metadata in values(graph.edge_data)
+            (; edge) = edge_metadata
+            id_source, _ = edge
+            if id_source.type in boundary_source_nodetypes
+                (; subnetwork_id) = graph[id_source]
+                # Check whether the source node is part of a subnetwork
+                if subnetwork_id ≠ 0
+                    subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id)
+                    mean_input_flows[subnetwork_idx][edge] = 0.0
+                end
+            end
         end
-    end
 
-    mean_realized_flows = Dict{Tuple{NodeID, NodeID}, Float64}()
+        # Find basins with a level demand
+        for node_id in values(graph.vertex_labels)
+            if has_external_demand(graph, node_id, :level_demand)[1]
+                subnetwork_id = graph[node_id].subnetwork_id
+                subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id)
+                mean_input_flows[subnetwork_idx][(node_id, node_id)] = 0.0
+            end
+        end
 
-    # Find edges that realize a demand
-    for edge_metadata in values(graph.edge_data)
-        (; type, edge) = edge_metadata
-
-        src_id, dst_id = edge
-        user_demand_inflow = (type == EdgeType.flow) && (dst_id.type == NodeType.UserDemand)
-        level_demand_inflow =
-            (type == EdgeType.control) && (src_id.type == NodeType.LevelDemand)
-        flow_demand_inflow =
-            (type == EdgeType.flow) && has_external_demand(graph, dst_id, :flow_demand)[1]
-
-        if user_demand_inflow || flow_demand_inflow
-            mean_realized_flows[edge] = 0.0
-        elseif level_demand_inflow
-            mean_realized_flows[(dst_id, dst_id)] = 0.0
+        # Find edges that realize a demand
+        for edge_metadata in values(graph.edge_data)
+            (; type, edge) = edge_metadata
+
+            src_id, dst_id = edge
+            user_demand_inflow =
+                (type == EdgeType.flow) && (dst_id.type == NodeType.UserDemand)
+            level_demand_inflow =
+                (type == EdgeType.control) && (src_id.type == NodeType.LevelDemand)
+            flow_demand_inflow =
+                (type == EdgeType.flow) &&
+                has_external_demand(graph, dst_id, :flow_demand)[1]
+
+            if user_demand_inflow || flow_demand_inflow
+                mean_realized_flows[edge] = 0.0
+            elseif level_demand_inflow
+                mean_realized_flows[(dst_id, dst_id)] = 0.0
+            end
         end
     end
 
diff --git a/core/src/util.jl b/core/src/util.jl
index 379f2abd5..2eac3f267 100644
--- a/core/src/util.jl
+++ b/core/src/util.jl
@@ -587,15 +587,17 @@ function set_initial_allocation_mean_flows!(integrator)::Nothing
     du = get_du(integrator)
     water_balance!(du, u, p, t)
 
-    for edge in keys(mean_input_flows)
-        if edge[1] == edge[2]
-            q = get_influx(du, edge[1], p)
-        else
-            q = get_flow(du, p, t, edge)
+    for mean_input_flows_subnetwork in values(mean_input_flows)
+        for edge in keys(mean_input_flows_subnetwork)
+            if edge[1] == edge[2]
+                q = get_influx(du, edge[1], p)
+            else
+                q = get_flow(du, p, t, edge)
+            end
+            # Multiply by Δt_allocation as averaging divides by this factor
+            # in update_allocation!
+            mean_input_flows_subnetwork[edge] = q * Δt_allocation
         end
-        # Multiply by Δt_allocation as averaging divides by this factor
-        # in update_allocation!
-        mean_input_flows[edge] = q * Δt_allocation
     end
 
     # Mean realized demands for basins are calculated as Δstorage/Δt
@@ -1004,12 +1006,8 @@ function set_state_flow_edges(p::Parameters, u0::ComponentVector)::Parameters
     state_inflow_edges = Vector{EdgeMetadata}[]
     state_outflow_edges = Vector{EdgeMetadata}[]
 
-    placeholder_edge = EdgeMetadata(
-        0,
-        EdgeType.flow,
-        0,
-        (NodeID(:Terminal, 0, 0), NodeID(:Terminal, 0, 0)),
-    )
+    placeholder_edge =
+        EdgeMetadata(0, EdgeType.flow, (NodeID(:Terminal, 0, 0), NodeID(:Terminal, 0, 0)))
 
     for node_name in keys(u0)
         if hasfield(Parameters, node_name)
@@ -1172,3 +1170,12 @@ function min_low_user_demand_level_factor(
         one(T)
     end
 end
+
+function mean_input_flows_subnetwork(p::Parameters, subnetwork_id::Int32)
+    (; mean_input_flows, subnetwork_ids) = p.allocation
+    subnetwork_idx = searchsortedfirst(subnetwork_ids, subnetwork_id)
+    return mean_input_flows[subnetwork_idx]
+end
+
+source_edges_subnetwork(p::Parameters, subnetwork_id::Int32) =
+    keys(mean_input_flows_subnetwork(p, subnetwork_id))
diff --git a/core/src/validation.jl b/core/src/validation.jl
index cfed36c6c..81a372421 100644
--- a/core/src/validation.jl
+++ b/core/src/validation.jl
@@ -182,6 +182,14 @@ function valid_nodes(db::DB)::Bool
     return !errors
 end
 
+function database_warning(db::DB)::Nothing
+    cols = SQLite.columns(db, "Edge")
+    if "subnetwork_id" in cols.name
+        @warn "The 'subnetwork_id' column in the 'Edge' table is deprecated since ribasim v2024.12."
+    end
+    return nothing
+end
+
 """
 Test for each node given its node type whether the nodes that
 # are downstream ('down-edge') of this node are of an allowed type
@@ -626,53 +634,6 @@ function valid_discrete_control(p::Parameters, config::Config)::Bool
     return !errors
 end
 
-"""
-An allocation source edge is valid if either:
-    - The edge connects the main network to a subnetwork
-    - The edge comes from a source node
-"""
-function valid_sources(
-    p::Parameters,
-    capacity::JuMP.Containers.SparseAxisArray{Float64, 2, Tuple{NodeID, NodeID}},
-    subnetwork_id::Int32,
-)::Bool
-    (; graph) = p
-
-    errors = false
-
-    # Loop over edges that were assigned a capacity
-    for edge in keys(capacity.data)
-
-        # For an edge (id_a, id_b) in the physical model
-        # the reverse (id_b, id_a) can exist in the allocation subnetwork
-        if !haskey(graph, edge...)
-            edge = reverse(edge)
-        end
-
-        (id_source, id_dst) = edge
-
-        # Whether the current edge is a source for the current subnetwork
-        if graph[edge...].subnetwork_id_source == subnetwork_id
-            from_source_node = id_source.type in allocation_source_nodetypes
-
-            if is_main_network(subnetwork_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].subnetwork_id)
-
-                if !from_source_node && !from_main_network
-                    errors = true
-                    @error "The source node of source edge $edge for subnetwork $subnetwork_id is neither a source node nor is it coming from the main network."
-                end
-            end
-        end
-    end
-    return !errors
-end
-
 function valid_priorities(priorities::Vector{Int32}, use_allocation::Bool)::Bool
     if use_allocation && any(iszero, priorities)
         return false
diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl
index 6da2fec36..3dd6c7ba0 100644
--- a/core/test/allocation_test.jl
+++ b/core/test/allocation_test.jl
@@ -14,7 +14,8 @@
 
     (; graph, allocation) = p
 
-    allocation.mean_input_flows[(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] = 4.5
+    allocation.mean_input_flows[1][(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] =
+        4.5
     allocation_model = p.allocation.allocation_models[1]
     (; flow) = allocation_model
     u = ComponentVector(; storage = zeros(length(p.basin.node_id)))
@@ -105,11 +106,11 @@ end
 
     # In each subnetwork, the connection from the main network to the subnetwork is
     # interpreted as a source
-    @test Ribasim.get_allocation_model(p, Int32(3)).problem[:source].axes[1] ==
+    @test Ribasim.get_allocation_model(p, Int32(3)).problem[:source_main_network].axes[1] ==
           [(NodeID(:Basin, 2, p), NodeID(:Pump, 11, p))]
-    @test Ribasim.get_allocation_model(p, Int32(5)).problem[:source].axes[1] ==
+    @test Ribasim.get_allocation_model(p, Int32(5)).problem[:source_main_network].axes[1] ==
           [(NodeID(:Basin, 6, p), NodeID(:Pump, 24, p))]
-    @test Ribasim.get_allocation_model(p, Int32(7)).problem[:source].axes[1] ==
+    @test Ribasim.get_allocation_model(p, Int32(7)).problem[:source_main_network].axes[1] ==
           [(NodeID(:Basin, 10, p), NodeID(:Pump, 38, p))]
 end
 
@@ -174,7 +175,7 @@ end
 
     # Running full allocation algorithm
     (; Δt_allocation) = allocation_models[1]
-    mean_input_flows[(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] =
+    mean_input_flows[1][(NodeID(:FlowBoundary, 1, p), NodeID(:Basin, 2, p))] =
         4.5 * Δt_allocation
     Ribasim.update_allocation!(model.integrator)
 
@@ -224,9 +225,9 @@ end
         allocation
     t = 0.0
 
-    # Set flows of sources in
-    mean_input_flows[(NodeID(:FlowBoundary, 58, p), NodeID(:Basin, 16, p))] = 1.0
-    mean_input_flows[(NodeID(:FlowBoundary, 59, p), NodeID(:Basin, 44, p))] = 1e-3
+    # Set flows of sources in subnetworks
+    mean_input_flows[2][(NodeID(:FlowBoundary, 58, p), NodeID(:Basin, 16, p))] = 1.0
+    mean_input_flows[4][(NodeID(:FlowBoundary, 59, p), NodeID(:Basin, 44, p))] = 1e-3
 
     # Collecting demands
     u = ComponentVector(; storage = zeros(length(basin.node_id)))
@@ -251,24 +252,24 @@ end
     du = get_du(model.integrator)
     Ribasim.formulate_storages!(current_storage, du, u, p, t)
 
-    current_storage ≈ Float32[
-        1.0346e6,
-        1.0346e6,
-        1.0346e6,
-        1.0346e6,
-        1.0346e6,
-        13.83,
-        40.10,
-        6.029e5,
-        4641,
-        2402,
-        6.03995,
-        928.8,
-        8.017,
-        10417,
-        5.619,
-        10417,
-        4.057,
+    @test current_storage ≈ Float32[
+        1.0346908f6,
+        1.03469f6,
+        1.0346894f6,
+        1.034689f6,
+        1.0346888f6,
+        13.833241,
+        40.109993,
+        187761.73,
+        4641.365,
+        2402.6687,
+        6.039952,
+        928.84283,
+        8.0175905,
+        10419.247,
+        5.619053,
+        10419.156,
+        4.057502,
     ]
 end
 
@@ -287,7 +288,7 @@ end
     (; user_demand, graph, allocation, basin, level_demand) = p
 
     # Initial "integrated" vertical flux
-    @test allocation.mean_input_flows[(NodeID(:Basin, 2, p), NodeID(:Basin, 2, p))] ≈ 1e2
+    @test allocation.mean_input_flows[1][(NodeID(:Basin, 2, p), NodeID(:Basin, 2, p))] ≈ 1e2
 
     Ribasim.solve!(model)
 
@@ -614,5 +615,5 @@ end
     # Given a max_level of Inf, the basin capacity is 0.0 because it is not possible for the basin level to be > Inf
     @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[1]) == 0.0
     @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[2]) == 0.0
-    @test Ribasim.get_basin_capacity(allocation_models[1], u, p, t, basin.node_id[3]) == 0.0
+    @test Ribasim.get_basin_capacity(allocation_models[2], u, p, t, basin.node_id[3]) == 0.0
 end
diff --git a/core/test/validation_test.jl b/core/test/validation_test.jl
index 547a11ff4..e46c0ef74 100644
--- a/core/test/validation_test.jl
+++ b/core/test/validation_test.jl
@@ -86,12 +86,7 @@ end
     graph[NodeID(:Pump, 6, 1)] = NodeMetadata(:pump, 9)
 
     function set_edge_metadata!(id_1, id_2, edge_type)
-        graph[id_1, id_2] = EdgeMetadata(;
-            id = 0,
-            type = edge_type,
-            subnetwork_id_source = 0,
-            edge = (id_1, id_2),
-        )
+        graph[id_1, id_2] = EdgeMetadata(; id = 0, type = edge_type, edge = (id_1, id_2))
         return nothing
     end
 
@@ -148,12 +143,7 @@ end
     graph[NodeID(:Basin, 7, 1)] = NodeMetadata(:basin, 0)
 
     function set_edge_metadata!(id_1, id_2, edge_type)
-        graph[id_1, id_2] = EdgeMetadata(;
-            id = 0,
-            type = edge_type,
-            subnetwork_id_source = 0,
-            edge = (id_1, id_2),
-        )
+        graph[id_1, id_2] = EdgeMetadata(; id = 0, type = edge_type, edge = (id_1, id_2))
         return nothing
     end
 
diff --git a/docs/changelog.qmd b/docs/changelog.qmd
index d9c04350e..373ab6cc8 100644
--- a/docs/changelog.qmd
+++ b/docs/changelog.qmd
@@ -5,6 +5,9 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 
 ## [Unreleased]
 
+### Changed
+- The Edge table no longer supports `subnetwork_id`; this is automatically inferred. [#1956](https://github.com/Deltares/Ribasim/pull/1956)
+
 ## [v2024.11.0] - 2024-10-08
 
 This major new release contains many improvements.
diff --git a/docs/concept/allocation.qmd b/docs/concept/allocation.qmd
index f90a83cec..91318fd84 100644
--- a/docs/concept/allocation.qmd
+++ b/docs/concept/allocation.qmd
@@ -42,9 +42,9 @@ The allocation problem is solved per subnetwork, which is given by a subset $S \
 
 Sources are indicated by a set of edges in the subnetwork
 $$
-E_S^\text{source} \subset E.
+E_S^\text{source} \subset E,
 $$
-That is, if $(i,j) \in E_S^\text{source}$, then the average over the last allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge
+which are automatically inferred as all edges that point out of LevelBoundary or FlowBoundary nodes. That is, if $(i,j) \in E_S^\text{source}$, then the average over the last allocation interval $\Delta t_{\text{alloc}}$ of the of the flow over this edge
 $$
     \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^tQ_{ij}(t') dt'
 $$
diff --git a/docs/dev/callstacks.qmd b/docs/dev/callstacks.qmd
index 5bd4fe372..a36c75cdf 100644
--- a/docs/dev/callstacks.qmd
+++ b/docs/dev/callstacks.qmd
@@ -74,10 +74,12 @@ toml_path = normpath(
     @__DIR__,
     "../../generated_testmodels/main_network_with_subnetworks/ribasim.toml",
 )
-config = Ribasim.Config(toml_path; allocation_use_allocation = false)
+config = Ribasim.Config(toml_path)
 db_path = Ribasim.database_path(config)
 db = SQLite.DB(db_path)
 p = Ribasim.Parameters(db, config)
+empty!(p.allocation.subnetwork_ids)
+empty!(p.allocation.main_network_connections)
 graph, verts = tracecall((Ribasim,), Ribasim.initialize_allocation!, (p, config))
 plot_graph(graph)
 ```
diff --git a/docs/guide/examples.ipynb b/docs/guide/examples.ipynb
index 31e899dd6..3e09e5986 100644
--- a/docs/guide/examples.ipynb
+++ b/docs/guide/examples.ipynb
@@ -1203,7 +1203,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1)\n",
+    "model.edge.add(model.flow_boundary[1], model.basin[2])\n",
     "model.edge.add(model.basin[2], model.user_demand[3])\n",
     "model.edge.add(model.basin[2], model.linear_resistance[4])\n",
     "model.edge.add(model.linear_resistance[4], model.basin[5])\n",
@@ -1499,7 +1499,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)\n",
+    "model.edge.add(model.flow_boundary[1], model.basin[2])\n",
     "model.edge.add(model.basin[2], model.user_demand[3])\n",
     "model.edge.add(model.level_demand[4], model.basin[2])\n",
     "model.edge.add(model.user_demand[3], model.basin[5])\n",
diff --git a/docs/reference/usage.qmd b/docs/reference/usage.qmd
index 63eeedd99..afe2f1529 100644
--- a/docs/reference/usage.qmd
+++ b/docs/reference/usage.qmd
@@ -200,7 +200,6 @@ to_node_id    | Int32                         | -
 edge_type     | String                        | must be "flow" or "control"
 geom          | LineString or MultiLineString | (optional)
 name          | String                        | (optional, does not have to be unique)
-subnetwork_id | Int32                         | (optional, denotes source in allocation network)
 
 Similarly to the node table, you can use a geometry to visualize the connections between the
 nodes in QGIS. For instance, you can draw a line connecting the two node coordinates.
diff --git a/python/ribasim/ribasim/geometry/edge.py b/python/ribasim/ribasim/geometry/edge.py
index dc8f6e202..bc473458f 100644
--- a/python/ribasim/ribasim/geometry/edge.py
+++ b/python/ribasim/ribasim/geometry/edge.py
@@ -47,7 +47,6 @@ class EdgeSchema(_GeoBaseSchema):
     from_node_id: Series[Int32] = pa.Field(default=0)
     to_node_id: Series[Int32] = pa.Field(default=0)
     edge_type: Series[str] = pa.Field(default="flow")
-    subnetwork_id: Series[pd.Int32Dtype] = pa.Field(default=pd.NA, nullable=True)
     geometry: GeoSeries[LineString] = pa.Field(default=None, nullable=True)
 
     @classmethod
@@ -73,7 +72,6 @@ def add(
         to_node: NodeData,
         geometry: LineString | MultiLineString | None = None,
         name: str = "",
-        subnetwork_id: int | None = None,
         edge_id: Optional[NonNegativeInt] = None,
         **kwargs,
     ):
@@ -92,9 +90,8 @@ def add(
             The geometry of a line. If not supplied, it creates a straight line between the nodes.
         name : str
             An optional name for the edge.
-        subnetwork_id : int | None
-            An optional subnetwork id for the edge. This edge indicates a source for
-            the allocation algorithm, and should thus not be set for every edge in a subnetwork.
+        id : int
+            An optional non-negative edge ID. If not supplied, it will be automatically generated.
         **kwargs : Dict
         """
         if not can_connect(from_node.node_type, to_node.node_type):
@@ -125,7 +122,6 @@ def add(
                 "to_node_id": [to_node.node_id],
                 "edge_type": [edge_type],
                 "name": [name],
-                "subnetwork_id": [subnetwork_id],
                 **kwargs,
             },
             geometry=geometry_to_append,
diff --git a/python/ribasim/ribasim/migrations.py b/python/ribasim/ribasim/migrations.py
index fc51eff52..31a326e83 100644
--- a/python/ribasim/ribasim/migrations.py
+++ b/python/ribasim/ribasim/migrations.py
@@ -27,6 +27,9 @@ def edgeschema_migration(gdf: GeoDataFrame, schema_version: int) -> GeoDataFrame
         warnings.warn("Migrating outdated Edge table.", UserWarning)
         assert gdf["edge_id"].is_unique, "Edge IDs have to be unique."
         gdf.set_index("edge_id", inplace=True)
+    if "subnetwork_id" in gdf.columns:
+        warnings.warn("Migrating outdated Edge table.", UserWarning)
+        gdf.drop(columns="subnetwork_id", inplace=True, errors="ignore")
 
     return gdf
 
diff --git a/python/ribasim/tests/test_model.py b/python/ribasim/tests/test_model.py
index ba703d5fa..f5e341561 100644
--- a/python/ribasim/tests/test_model.py
+++ b/python/ribasim/tests/test_model.py
@@ -145,7 +145,6 @@ def test_edge_table(basic):
     df = model.edge.df
     assert df.geometry.is_unique
     assert df.from_node_id.dtype == np.int32
-    assert df.subnetwork_id.dtype == pd.Int32Dtype()
     assert df.crs == CRS.from_epsg(28992)
 
 
diff --git a/python/ribasim_testmodels/ribasim_testmodels/allocation.py b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
index dce660924..dce319344 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/allocation.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/allocation.py
@@ -153,7 +153,7 @@ def subnetwork_model() -> Model:
     )
     model.outlet.add(Node(13, Point(2, 4), subnetwork_id=2), outlet_data)
 
-    model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)
+    model.edge.add(model.flow_boundary[1], model.basin[2])
     model.edge.add(model.basin[2], model.outlet[3])
     model.edge.add(model.outlet[3], model.terminal[4])
     model.edge.add(model.basin[2], model.user_demand[10])
@@ -268,7 +268,7 @@ def looped_subnetwork_model() -> Model:
         ],
     )
 
-    model.edge.add(model.flow_boundary[5], model.basin[2], subnetwork_id=2)
+    model.edge.add(model.flow_boundary[5], model.basin[2])
     model.edge.add(model.basin[2], model.outlet[3])
     model.edge.add(model.outlet[3], model.terminal[4])
     model.edge.add(model.basin[2], model.user_demand[1])
@@ -346,7 +346,7 @@ def minimal_subnetwork_model() -> Model:
         ],
     )
 
-    model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)
+    model.edge.add(model.flow_boundary[1], model.basin[2])
     model.edge.add(model.basin[2], model.pump[3])
     model.edge.add(model.pump[3], model.basin[4])
     model.edge.add(model.basin[4], model.user_demand[5])
@@ -406,7 +406,7 @@ def allocation_example_model() -> Model:
     )
     model.terminal.add(Node(8, Point(5, 0), subnetwork_id=2))
 
-    model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)
+    model.edge.add(model.flow_boundary[1], model.basin[2])
     model.edge.add(model.basin[2], model.user_demand[3])
     model.edge.add(model.basin[2], model.linear_resistance[4])
     model.edge.add(model.linear_resistance[4], model.basin[5])
@@ -628,7 +628,7 @@ def main_network_with_subnetworks_model() -> Model:
         [user_demand.Static(return_factor=[0.9], priority=2, min_level=0.0)],
     )
 
-    model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=1)
+    model.edge.add(model.flow_boundary[1], model.basin[2])
     model.edge.add(model.basin[2], model.linear_resistance[3])
     model.edge.add(model.linear_resistance[3], model.basin[4])
     model.edge.add(model.basin[4], model.linear_resistance[5])
@@ -686,9 +686,9 @@ def main_network_with_subnetworks_model() -> Model:
     model.edge.add(model.user_demand[53], model.basin[50])
     model.edge.add(model.basin[44], model.user_demand[57])
     model.edge.add(model.user_demand[57], model.basin[44])
-    model.edge.add(model.basin[2], model.pump[11], subnetwork_id=3)
-    model.edge.add(model.basin[6], model.pump[24], subnetwork_id=5)
-    model.edge.add(model.basin[10], model.pump[38], subnetwork_id=7)
+    model.edge.add(model.basin[2], model.pump[11])
+    model.edge.add(model.basin[6], model.pump[24])
+    model.edge.add(model.basin[10], model.pump[38])
     model.edge.add(model.basin[8], model.user_demand[60])
     model.edge.add(model.user_demand[60], model.basin[8])
 
@@ -708,8 +708,8 @@ def subnetworks_with_sources_model() -> Model:
         [flow_boundary.Static(flow_rate=[0.003])],
     )
 
-    model.edge.add(model.flow_boundary[58], model.basin[16], subnetwork_id=3)
-    model.edge.add(model.flow_boundary[59], model.basin[44], subnetwork_id=7)
+    model.edge.add(model.flow_boundary[58], model.basin[16])
+    model.edge.add(model.flow_boundary[59], model.basin[44])
 
     return model
 
@@ -763,7 +763,7 @@ def level_demand_model() -> Model:
         [basin.Profile(area=1000.0, level=[0.0, 1.0]), basin.State(level=[2.0])],
     )
 
-    model.edge.add(model.flow_boundary[1], model.basin[2], subnetwork_id=2)
+    model.edge.add(model.flow_boundary[1], model.basin[2])
     model.edge.add(model.basin[2], model.user_demand[3])
     model.edge.add(model.level_demand[4], model.basin[2])
     model.edge.add(model.user_demand[3], model.basin[5])
@@ -835,7 +835,6 @@ def flow_demand_model() -> Model:
     model.edge.add(
         model.level_boundary[1],
         model.tabulated_rating_curve[2],
-        subnetwork_id=2,
     )
     model.edge.add(model.tabulated_rating_curve[2], model.basin[3])
     model.edge.add(model.basin[3], model.user_demand[4])
@@ -877,7 +876,7 @@ def linear_resistance_demand_model():
         [flow_demand.Static(priority=[1], demand=2.0)],
     )
 
-    model.edge.add(model.basin[1], model.linear_resistance[2], subnetwork_id=1)
+    model.edge.add(model.basin[1], model.linear_resistance[2])
     model.edge.add(model.linear_resistance[2], model.basin[3])
     model.edge.add(model.flow_demand[4], model.linear_resistance[2])
 
@@ -966,7 +965,7 @@ def fair_distribution_model():
         ],
     )
 
-    model.edge.add(model.level_boundary[1], model.pump[2], subnetwork_id=1)
+    model.edge.add(model.level_boundary[1], model.pump[2])
     model.edge.add(model.pump[2], model.basin[3])
     model.edge.add(model.basin[3], model.linear_resistance[4])
     model.edge.add(model.linear_resistance[4], model.basin[5])
diff --git a/python/ribasim_testmodels/ribasim_testmodels/invalid.py b/python/ribasim_testmodels/ribasim_testmodels/invalid.py
index 049562e6d..1e0657f6e 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/invalid.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/invalid.py
@@ -214,7 +214,6 @@ def invalid_priorities_model() -> Model:
     model.edge.add(
         model.level_boundary[1],
         model.tabulated_rating_curve[2],
-        subnetwork_id=2,
     )
     model.edge.add(model.tabulated_rating_curve[2], model.basin[3])
     model.edge.add(model.basin[3], model.user_demand[4])
diff --git a/ribasim_qgis/core/nodes.py b/ribasim_qgis/core/nodes.py
index 0f6eb5e7e..35d8f3970 100644
--- a/ribasim_qgis/core/nodes.py
+++ b/ribasim_qgis/core/nodes.py
@@ -254,7 +254,6 @@ def attributes(cls) -> list[QgsField]:
             QgsField("from_node_id", QVariant.Int),
             QgsField("to_node_id", QVariant.Int),
             QgsField("edge_type", QVariant.String),
-            QgsField("subnetwork_id", QVariant.Int),
         ]
 
     @classmethod