diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl
index 9bf3ddb30..6ac9f2be0 100644
--- a/core/src/allocation_init.jl
+++ b/core/src/allocation_init.jl
@@ -153,7 +153,8 @@ function add_variables_basin!(
     (; graph) = p
     node_ids_basin = [
         node_id for
-        node_id in graph[].node_ids[subnetwork_id] if graph[node_id].type == :basin
+        node_id in graph[].node_ids[subnetwork_id] if graph[node_id].type == :basin &&
+        has_external_demand(graph, node_id, :level_demand)[1]
     ]
     problem[:F_basin_in] =
         JuMP.@variable(problem, F_basin_in[node_id = node_ids_basin,] >= 0.0)
@@ -211,7 +212,7 @@ function add_variables_absolute_value!(
         type = node_id.type
         if type == NodeType.UserDemand
             push!(node_ids_user_demand, node_id)
-        elseif type == NodeType.Basin
+        elseif has_external_demand(graph, node_id, :level_demand)[1]
             push!(node_ids_level_demand, node_id)
         elseif has_external_demand(graph, node_id, :flow_demand)[1]
             push!(node_ids_flow_demand, node_id)
@@ -392,8 +393,8 @@ function add_constraints_conservation_node!(
             end
         end
 
-        # If the node is a Basin, add basin in- and outflow
-        if node_id.type == NodeType.Basin
+        # If the node is a Basin with a level demand, add basin in- and outflow
+        if has_external_demand(graph, node_id, :level_demand)[1]
             push!(inflows_node, F_basin_out[node_id])
             push!(outflows_node, F_basin_in[node_id])
         end
diff --git a/core/src/allocation_optim.jl b/core/src/allocation_optim.jl
index d3c94f262..01370c228 100644
--- a/core/src/allocation_optim.jl
+++ b/core/src/allocation_optim.jl
@@ -254,7 +254,8 @@ function set_initial_capacities_source!(
     p::Parameters,
 )::Nothing
     (; problem) = allocation_model
-    (; graph) = p
+    (; graph, allocation) = p
+    (; mean_flows) = allocation
     (; subnetwork_id) = allocation_model
     source_constraints = problem[:source]
     main_network_source_edges = get_main_network_connections(p, subnetwork_id)
@@ -264,8 +265,8 @@ function set_initial_capacities_source!(
         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 current flow from the physical layer.
-                source_capacity = get_flow(graph, edge..., 0)
+                # Reset the source to the averaged flow over the last allocation period
+                source_capacity = mean_flows[edge][]
                 JuMP.set_normalized_rhs(
                     source_constraints[edge],
                     # It is assumed that the allocation procedure does not have to be differentiated.
@@ -357,14 +358,14 @@ function get_basin_data(
     u::ComponentVector,
     node_id::NodeID,
 )
-    (; graph, basin, level_demand) = p
+    (; graph, basin, level_demand, allocation) = p
     (; vertical_flux) = basin
     (; Δt_allocation) = allocation_model
+    (; mean_flows) = allocation
     @assert node_id.type == NodeType.Basin
     vertical_flux = get_tmp(vertical_flux, 0)
     _, basin_idx = id_index(basin.node_id, node_id)
-    # NOTE: Instantaneous
-    influx = get_influx(basin, node_id)
+    influx = mean_flows[(node_id, node_id)][]
     _, basin_idx = id_index(basin.node_id, node_id)
     storage_basin = u.storage[basin_idx]
     control_inneighbors = inneighbor_labels_type(graph, node_id, EdgeType.control)
@@ -511,12 +512,15 @@ function set_initial_demands_level!(
     p::Parameters,
     t::Float64,
 )::Nothing
-    (; subnetwork_id) = allocation_model
+    (; subnetwork_id, problem) = allocation_model
     (; graph, basin) = p
     (; node_id, demand) = basin
 
-    for (i, id) in enumerate(node_id)
+    node_ids_level_demand = only(problem[:basin_outflow].axes)
+
+    for id in node_ids_level_demand
         if graph[id].subnetwork_id == subnetwork_id
+            _, i = id_index(node_id, id)
             demand[i] = get_basin_demand(allocation_model, u, p, t, id)
         end
     end
@@ -606,8 +610,9 @@ function adjust_demands_level!(allocation_model::AllocationModel, p::Parameters)
     F_basin_in = problem[:F_basin_in]
 
     # Reduce the demand by what was allocated
-    for (i, id) in enumerate(node_id)
+    for id in only(F_basin_in.axes)
         if graph[id].subnetwork_id == subnetwork_id
+            _, i = id_index(basin.node_id, id)
             demand[i] -= JuMP.value(F_basin_in[id])
         end
     end
@@ -763,8 +768,7 @@ function save_demands_and_allocations!(
             #NOTE: instantaneous
             realized = get_flow(graph, inflow_id(graph, node_id), node_id, 0)
 
-        elseif node_id.type == NodeType.Basin &&
-               has_external_demand(graph, node_id, :level_demand)[1]
+        elseif has_external_demand(graph, node_id, :level_demand)[1]
             basin_priority_idx = get_external_priority_idx(p, node_id)
 
             if priority_idx == 1 || basin_priority_idx == priority_idx
diff --git a/core/src/callback.jl b/core/src/callback.jl
index 51dd3d029..048ba7530 100644
--- a/core/src/callback.jl
+++ b/core/src/callback.jl
@@ -91,7 +91,7 @@ Integrate flows over the last timestep
 """
 function integrate_flows!(u, t, integrator)::Nothing
     (; p, dt) = integrator
-    (; graph, user_demand, basin) = p
+    (; graph, user_demand, basin, allocation) = p
     (; flow, flow_dict, flow_prev, flow_integrated) = graph[]
     (; vertical_flux, vertical_flux_prev, vertical_flux_integrated, vertical_flux_bmi) =
         basin
@@ -106,12 +106,31 @@ function integrate_flows!(u, t, integrator)::Nothing
     @. vertical_flux_integrated += 0.5 * (vertical_flux + vertical_flux_prev) * dt
     @. vertical_flux_bmi += 0.5 * (vertical_flux + vertical_flux_prev) * dt
 
+    # UserDemand realized flows for BMI
     for (i, id) in enumerate(user_demand.node_id)
         src_id = inflow_id(graph, id)
         flow_idx = flow_dict[src_id, id]
         user_demand.realized_bmi[i] += 0.5 * (flow[flow_idx] + flow_prev[flow_idx]) * dt
     end
 
+    # Allocation source flows
+    for (edge, value) in allocation.mean_flows
+        if edge[1] == edge[2]
+            # Vertical fluxes
+            _, basin_idx = id_index(basin.node_id, edge[1])
+            value[] +=
+                0.5 *
+                (get_influx(basin, basin_idx) + get_influx(basin, basin_idx; prev = true)) *
+                dt
+        else
+            # Horizontal flows
+            value[] +=
+                0.5 *
+                (get_flow(graph, edge..., 0) + get_flow(graph, edge..., 0; prev = true)) *
+                dt
+        end
+    end
+
     copyto!(flow_prev, flow)
     copyto!(vertical_flux_prev, vertical_flux)
     return nothing
@@ -446,7 +465,14 @@ end
 function update_allocation!(integrator)::Nothing
     (; p, t, u) = integrator
     (; allocation) = p
-    (; allocation_models) = allocation
+    (; allocation_models, mean_flows) = allocation
+    (; Δt_allocation) = allocation_models[1]
+
+    # Divide by the allocation Δt to obtain the mean flows
+    # from the integrated flows
+    for value in values(mean_flows)
+        value[] /= Δt_allocation
+    end
 
     # If a main network is present, collect demands of subnetworks
     if has_main_network(allocation)
@@ -462,6 +488,11 @@ function update_allocation!(integrator)::Nothing
     for allocation_model in allocation_models
         allocate!(p, allocation_model, t, u, OptimizationType.allocate)
     end
+
+    # Reset the mean source flows
+    for value in values(mean_flows)
+        value[] = 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 19215b250..edcd8f3cd 100644
--- a/core/src/graph.jl
+++ b/core/src/graph.jl
@@ -33,16 +33,15 @@ function create_graph(db::DB, config::Config, chunk_sizes::Vector{Int})::MetaGra
         node_id = NodeID(row.node_type, row.node_id)
         # Process allocation network ID
         if ismissing(row.subnetwork_id)
-            allocation_network_id = 0
+            subnetwork_id = 0
         else
-            allocation_network_id = row.subnetwork_id
-            if !haskey(node_ids, allocation_network_id)
-                node_ids[allocation_network_id] = Set{NodeID}()
+            subnetwork_id = row.subnetwork_id
+            if !haskey(node_ids, subnetwork_id)
+                node_ids[subnetwork_id] = Set{NodeID}()
             end
-            push!(node_ids[allocation_network_id], node_id)
+            push!(node_ids[subnetwork_id], node_id)
         end
-        graph[node_id] =
-            NodeMetadata(Symbol(snake_case(row.node_type)), allocation_network_id)
+        graph[node_id] = NodeMetadata(Symbol(snake_case(row.node_type)), subnetwork_id)
     end
 
     errors = false
@@ -178,9 +177,16 @@ end
 """
 Get the flow over the given edge (val is needed for get_tmp from ForwardDiff.jl).
 """
-function get_flow(graph::MetaGraph, id_src::NodeID, id_dst::NodeID, val)::Number
-    (; flow_dict, flow) = graph[]
-    return get_tmp(flow, val)[flow_dict[id_src, id_dst]]
+function get_flow(
+    graph::MetaGraph,
+    id_src::NodeID,
+    id_dst::NodeID,
+    val;
+    prev::Bool = false,
+)::Number
+    (; flow_dict, flow, flow_prev) = graph[]
+    flow_vector = prev ? flow_prev : flow
+    return get_tmp(flow_vector, val)[flow_dict[id_src, id_dst]]
 end
 
 """
diff --git a/core/src/parameter.jl b/core/src/parameter.jl
index 2566eabb4..d16b793a8 100644
--- a/core/src/parameter.jl
+++ b/core/src/parameter.jl
@@ -43,7 +43,7 @@ const VectorInterpolation =
 """
 Store information for a subnetwork used for allocation.
 
-allocation_network_id: The ID of this allocation network
+subnetwork_id: The ID of this allocation network
 capacity: The capacity per edge of the allocation network, as constrained by nodes that have a max_flow_rate
 problem: The JuMP.jl model for solving the allocation problem
 Δt_allocation: The time interval between consecutive allocation solves
@@ -57,13 +57,15 @@ end
 
 """
 Object for all information about allocation
-allocation_network_ids: The unique sorted allocation network IDs
+subnetwork_ids: The unique sorted allocation network IDs
 allocation models: The allocation models for the main network and subnetworks corresponding to
-    allocation_network_ids
+    subnetwork_ids
 main_network_connections: (from_id, to_id) from the main network to the subnetwork per subnetwork
 priorities: All used priority values.
 subnetwork_demands: The demand of an edge from the main network to a subnetwork
-record_demand: A record of demands and allocated flows for nodes that have these.
+subnetwork_allocateds: The allocated flow of an edge from the main network to a subnetwork
+mean_flows: Flows averaged over Δt_allocation over edges that are allocation sources
+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
 """
@@ -74,6 +76,7 @@ struct Allocation
     priorities::Vector{Int32}
     subnetwork_demands::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
     subnetwork_allocateds::Dict{Tuple{NodeID, NodeID}, Vector{Float64}}
+    mean_flows::Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}
     record_demand::@NamedTuple{
         time::Vector{Float64},
         subnetwork_id::Vector{Int32},
@@ -103,7 +106,7 @@ is_active(allocation::Allocation) = !isempty(allocation.allocation_models)
 """
 Type for storing metadata of nodes in the graph
 type: type of the node
-allocation_network_id: Allocation network ID (0 if not in subnetwork)
+subnetwork_id: Allocation network ID (0 if not in subnetwork)
 """
 struct NodeMetadata
     type::Symbol
diff --git a/core/src/read.jl b/core/src/read.jl
index 7d3ff0a69..9b84fed30 100644
--- a/core/src/read.jl
+++ b/core/src/read.jl
@@ -951,7 +951,7 @@ function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid
     return Subgrid(basin_ids, interpolations, fill(NaN, length(basin_ids)))
 end
 
-function Allocation(db::DB, config::Config)::Allocation
+function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation
     record_demand = (
         time = Float64[],
         subnetwork_id = Int32[],
@@ -976,18 +976,34 @@ function Allocation(db::DB, config::Config)::Allocation
         optimization_type = String[],
     )
 
-    allocation = Allocation(
+    mean_flows = Dict{Tuple{NodeID, NodeID}, Base.RefValue{Float64}}()
+
+    # 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_flows[edge] = Ref(0.0)
+        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_flows[(node_id, node_id)] = Ref(0.0)
+        end
+    end
+
+    return Allocation(
         Int32[],
         AllocationModel[],
         Vector{Tuple{NodeID, NodeID}}[],
         get_all_priorities(db, config),
-        Dict{Tuple{NodeID, NodeID}, Float64}(),
-        Dict{Tuple{NodeID, NodeID}, Float64}(),
+        Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(),
+        Dict{Tuple{NodeID, NodeID}, Vector{Float64}}(),
+        mean_flows,
         record_demand,
         record_flow,
     )
-
-    return allocation
 end
 
 """
@@ -1007,7 +1023,7 @@ 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 = Allocation(db, config)
+    allocation = Allocation(db, config, graph)
 
     if !valid_edges(graph)
         error("Invalid edge(s) found.")
diff --git a/core/src/util.jl b/core/src/util.jl
index a52bdb277..b1f816d9c 100644
--- a/core/src/util.jl
+++ b/core/src/util.jl
@@ -711,10 +711,11 @@ function get_influx(basin::Basin, node_id::NodeID)::Float64
     return get_influx(basin, basin_idx)
 end
 
-function get_influx(basin::Basin, basin_idx::Int)::Float64
-    (; vertical_flux) = basin
+function get_influx(basin::Basin, basin_idx::Int; prev::Bool = false)::Float64
+    (; vertical_flux, vertical_flux_prev) = basin
     vertical_flux = get_tmp(vertical_flux, 0)
-    (; precipitation, evaporation, drainage, infiltration) = vertical_flux
+    flux_vector = prev ? vertical_flux_prev : vertical_flux
+    (; precipitation, evaporation, drainage, infiltration) = flux_vector
     return precipitation[basin_idx] - evaporation[basin_idx] + drainage[basin_idx] -
            infiltration[basin_idx]
 end
diff --git a/core/test/allocation_test.jl b/core/test/allocation_test.jl
index 68c1a262b..d49f0512d 100644
--- a/core/test/allocation_test.jl
+++ b/core/test/allocation_test.jl
@@ -11,10 +11,10 @@
     db = SQLite.DB(db_path)
 
     p = Ribasim.Parameters(db, cfg)
-    graph = p.graph
+    (; graph, allocation) = p
     close(db)
 
-    Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5) # Source flow
+    allocation.mean_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5
     allocation_model = p.allocation.allocation_models[1]
     u = ComponentVector(; storage = zeros(length(p.basin.node_id)))
     Ribasim.allocate!(p, allocation_model, 0.0, u, OptimizationType.allocate)
@@ -210,8 +210,13 @@ end
     close(db)
 
     (; allocation, user_demand, graph, basin) = p
-    (; allocation_models, subnetwork_demands, subnetwork_allocateds, record_flow) =
-        allocation
+    (;
+        allocation_models,
+        subnetwork_demands,
+        subnetwork_allocateds,
+        record_flow,
+        mean_flows,
+    ) = allocation
     t = 0.0
 
     # Collecting demands
@@ -241,7 +246,8 @@ end
     @test F_abs_user_demand[NodeID(:Pump, 38)] ∈ objective_variables
 
     # Running full allocation algorithm
-    Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 1), NodeID(:Basin, 2), 4.5)
+    (; Δt_allocation) = allocation_models[1]
+    mean_flows[(NodeID(:FlowBoundary, 1), NodeID(:Basin, 2))][] = 4.5 * Δt_allocation
     u = ComponentVector(; storage = zeros(length(p.basin.node_id)))
     Ribasim.update_allocation!((; p, t, u))
 
@@ -283,12 +289,13 @@ end
     close(db)
 
     (; allocation, user_demand, graph, basin) = p
-    (; allocation_models, subnetwork_demands, subnetwork_allocateds) = allocation
+    (; allocation_models, subnetwork_demands, subnetwork_allocateds, mean_flows) =
+        allocation
     t = 0.0
 
-    # Set flows of sources in subnetworks
-    Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 58), NodeID(:Basin, 16), 1.0)
-    Ribasim.set_flow!(graph, NodeID(:FlowBoundary, 59), NodeID(:Basin, 44), 1e-3)
+    # Set flows of sources in
+    mean_flows[(NodeID(:FlowBoundary, 58), NodeID(:Basin, 16))][] = 1.0
+    mean_flows[(NodeID(:FlowBoundary, 59), NodeID(:Basin, 44))][] = 1e-3
 
     # Collecting demands
     u = ComponentVector(; storage = zeros(length(basin.node_id)))
@@ -372,7 +379,7 @@ end
     @test storage[stage_6] ≈ u_stage_6.(t[stage_6]) rtol = 1e-4
 end
 
-@testitem "flow_demand" begin
+@testitem "Flow demand" begin
     using JuMP
     using Ribasim: NodeID, OptimizationType
 
@@ -422,6 +429,9 @@ end
     t = 0.0
     (; u) = model.integrator
     optimization_type = OptimizationType.internal_sources
+    for (edge, value) in allocation.mean_flows
+        value[] = Ribasim.get_flow(graph, edge..., 0)
+    end
     Ribasim.set_initial_values!(allocation_model, p, u, t)
 
     # Priority 1
@@ -484,8 +494,8 @@ end
     )
     # Get demand from buffers
     d = user_demand.demand_itp[3][4](t)
-    @assert JuMP.value(F[(NodeID(:UserDemand, 4), NodeID(:Basin, 7))]) +
-            JuMP.value(F[(NodeID(:UserDemand, 6), NodeID(:Basin, 7))]) == d
+    @test JuMP.value(F[(NodeID(:UserDemand, 4), NodeID(:Basin, 7))]) +
+          JuMP.value(F[(NodeID(:UserDemand, 6), NodeID(:Basin, 7))]) == d
 end
 
 @testitem "flow_demand_with_max_flow_rate" begin
diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd
index 9a5947482..1399513b3 100644
--- a/docs/core/allocation.qmd
+++ b/docs/core/allocation.qmd
@@ -42,7 +42,11 @@ Sources are indicated by a set of edges in the subnetwork
 $$
 E_S^\text{source} \subset E.
 $$
-That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork.
+That is, if $(i,j) \in E_S^\text{source}$, then the average over the 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'
+$$
+ is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork. For the definition of $Q_{ij}$ see [the formal model description](equations.qmd#formal-model-description).
 
 ### User demands
 
@@ -71,12 +75,12 @@ for all $i \in FD_S$. Here $d^{p_{\text{df}}}$ is given by the original flow dem
 
 ### Vertical fluxes and local storage
 
-Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the sum of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin
+Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the average over the allocation interval $\Delta t_{\text{alloc}}$ of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin:
 $$
-    \phi_i(t), \quad \forall i \in B_S.
+    \phi_i(t) = \frac{1}{\Delta t_{\text{alloc}}}\int_{t - \Delta t_{\text{alloc}}}^t \left[Q_{P,i}(t') - Q_{E,i}(t') + Q_{\text{drn},i}(t') - Q_{\text{inf},i}(t') \right] dt', \quad \forall i \in B_S.
 $$
 
-We consider fluxes into the basin to be positive and out of the basin to be negative.
+We consider fluxes into the basin to be positive and out of the basin to be negative. For more information see [the natural water balance terms](equations.qmd#natural-water-balance-terms).
 
 Secondly, there is either a supply or demand from the storage in the basin. Given a minimum level $\ell_{\min, i}$ and a maximum level $\ell_{\max, i}$ which correspond to a minimum storage $s_{\min, i}$ and maximum storage $s_{\max, i}$ respectively, we get a flow supply of
 $$
@@ -88,7 +92,7 @@ $$
     d^p_i = \max\left(0.0, \frac{s_{\min,i} - u_i(t)}{\Delta t_{\text{alloc}}} - \phi_i(t)\right),
 $$
 
-for all $i \in B_S$. Here $\Delta t_{\text{alloc}}$ is the simulated time between two consecutive allocation solves. Note that the basin demand has only a single priority, so for other priorities this demand is $0$.
+for all $i \in B_S$. Note that the basin demand has only a single priority, so for other priorities this demand is $0$.
 
 ### Constraining factors
 
diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd
index 1f395cd58..a2719a638 100644
--- a/docs/core/equations.qmd
+++ b/docs/core/equations.qmd
@@ -137,7 +137,7 @@ ax.hlines([0,1],x_min,x_max, color = "k", ls = ":", zorder=-1)
 ax.vlines([0,p], *y_lim, color = "k", ls = ":")
 ax.set_xlim(x_min,x_max)
 ax.set_xlabel("$x$", fontsize=fontsize)
-ax.set_ylabel("$\phi(x;p)$", fontsize=fontsize)
+ax.set_ylabel(r"$\phi(x;p)$", fontsize=fontsize)
 ax.set_ylim(y_lim)
 
 fig.tight_layout()