diff --git a/.github/workflows/python_tests.yml b/.github/workflows/python_tests.yml
index b5e850664..e66bfeece 100644
--- a/.github/workflows/python_tests.yml
+++ b/.github/workflows/python_tests.yml
@@ -11,7 +11,7 @@ concurrency:
   cancel-in-progress: true
 jobs:
   test:
-    name: Python ${{ matrix.python_version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    name: Python ${{ matrix.python-version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
@@ -20,8 +20,7 @@ jobs:
           - ubuntu-latest
           - macOS-latest
           - windows-latest
-        python_version:
-          - "3.9"
+        python-version:
           - "3.10"
           - "3.11"
         arch:
@@ -29,15 +28,18 @@ jobs:
     steps:
       - uses: actions/checkout@v4
 
-      - uses: prefix-dev/setup-pixi@v0.3.0
+      - name: Set up Python ${{ matrix.python-version }}
+        uses: actions/setup-python@v4
         with:
-          pixi-version: "latest"
-          cache: true
-      - name: Prepare pixi
-        run: pixi run install-without-pre-commit
+          python-version: "${{ matrix.python-version }}"
+
+      - name: Install test dependencies
+        run: |
+          pip install --editable "python/ribasim_testmodels"
+          pip install --editable "python/ribasim[tests]"
 
       - name: Run tests
-        run: pixi run test-ribasim-python-cov
+        run: pytest --numprocesses=auto --cov=ribasim --cov-report=xml python/ribasim/tests
 
       - name: Upload coverage to Codecov
         uses: codecov/codecov-action@v3
diff --git a/core/Project.toml b/core/Project.toml
index 7650af4b3..07e5e9321 100644
--- a/core/Project.toml
+++ b/core/Project.toml
@@ -44,6 +44,7 @@ CodecZstd = "0.7,0.8"
 ComponentArrays = "0.13.14, 0.14, 0.15"
 Configurations = "0.17"
 DBInterface = "2.4"
+DataFrames = "1.4"
 DataInterpolations = "3.7, 4"
 DataStructures = "0.18"
 Dictionaries = "0.3.25"
@@ -68,6 +69,7 @@ julia = "1.9"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IOCapture = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
@@ -78,4 +80,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestReports = "dcd651b4-b50a-5b6b-8f22-87e9f253a252"
 
 [targets]
-test = ["Aqua", "CSV", "Documenter", "IOCapture", "Logging", "SafeTestsets", "TerminalLoggers", "Test", "TestReports", "TOML"]
+test = ["Aqua", "CSV", "DataFrames", "Documenter", "IOCapture", "Logging", "SafeTestsets", "TerminalLoggers", "Test", "TestReports", "TOML"]
diff --git a/core/src/bmi.jl b/core/src/bmi.jl
index c35e820ed..13ce7a70e 100644
--- a/core/src/bmi.jl
+++ b/core/src/bmi.jl
@@ -159,10 +159,25 @@ end
 Write all output to the configured output files.
 """
 function BMI.finalize(model::Model)::Model
-    compress = get_compressor(model.config.output)
-    write_basin_output(model, compress)
-    write_flow_output(model, compress)
-    write_discrete_control_output(model, compress)
+    (; config) = model
+    (; output) = model.config
+    compress = get_compressor(output)
+
+    # basin
+    table = basin_table(model)
+    path = output_path(config, output.basin)
+    write_arrow(path, table, compress)
+
+    # flow
+    table = flow_table(model)
+    path = output_path(config, output.flow)
+    write_arrow(path, table, compress)
+
+    # discrete control
+    table = discrete_control_table(model)
+    path = output_path(config, output.control)
+    write_arrow(path, table, compress)
+
     @debug "Wrote output."
     return model
 end
@@ -198,20 +213,20 @@ function create_callbacks(
     saveat,
 )::Tuple{CallbackSet, SavedValues{Float64, Vector{Float64}}}
     (; starttime, basin, tabulated_rating_curve, discrete_control) = parameters
+    callbacks = SciMLBase.DECallback[]
 
     tstops = get_tstops(basin.time.time, starttime)
     basin_cb = PresetTimeCallback(tstops, update_basin)
+    push!(callbacks, basin_cb)
 
     tstops = get_tstops(tabulated_rating_curve.time.time, starttime)
     tabulated_rating_curve_cb = PresetTimeCallback(tstops, update_tabulated_rating_curve!)
+    push!(callbacks, tabulated_rating_curve_cb)
 
-    # add a single time step's contribution to the water balance step's totals
-    # trackwb_cb = FunctionCallingCallback(track_waterbalance!)
-    # flows: save the flows over time, as a Vector of the nonzeros(flow)
-
+    # save the flows over time, as a Vector of the nonzeros(flow)
     saved_flow = SavedValues(Float64, Vector{Float64})
-
     save_flow_cb = SavingCallback(save_flow, saved_flow; saveat, save_start = false)
+    push!(callbacks, save_flow_cb)
 
     n_conditions = length(discrete_control.node_id)
     if n_conditions > 0
@@ -221,15 +236,9 @@ function create_callbacks(
             discrete_control_affect_downcrossing!,
             n_conditions,
         )
-        callback = CallbackSet(
-            save_flow_cb,
-            basin_cb,
-            tabulated_rating_curve_cb,
-            discrete_control_cb,
-        )
-    else
-        callback = CallbackSet(save_flow_cb, basin_cb, tabulated_rating_curve_cb)
+        push!(callbacks, discrete_control_cb)
     end
+    callback = CallbackSet(callbacks...)
 
     return callback, saved_flow
 end
diff --git a/core/src/create.jl b/core/src/create.jl
index 640e936bb..8543d238d 100644
--- a/core/src/create.jl
+++ b/core/src/create.jl
@@ -33,6 +33,7 @@ function parse_static_and_time(
     vals_out = []
 
     node_ids = get_ids(db, nodetype)
+    node_names = get_names(db, nodetype)
     n_nodes = length(node_ids)
 
     # Initialize the vectors for the output
@@ -91,7 +92,7 @@ function parse_static_and_time(
     t_end = seconds_since(config.endtime, config.starttime)
     trivial_timespan = [nextfloat(-Inf), prevfloat(Inf)]
 
-    for (node_idx, node_id) in enumerate(node_ids)
+    for (node_idx, (node_id, node_name)) in enumerate(zip(node_ids, node_names))
         if node_id in static_node_ids
             # The interval of rows of the static table that have the current node_id
             rows = searchsorted(static.node_id, node_id)
@@ -153,7 +154,7 @@ function parse_static_and_time(
                     )
                     if !is_valid
                         errors = true
-                        @error "A $parameter_name time series for $nodetype node #$node_id has repeated times, this can not be interpolated."
+                        @error "A $parameter_name time series for $nodetype node $(repr(node_name)) (#$node_id) has repeated times, this can not be interpolated."
                     end
                 else
                     # Activity of transient nodes is assumed to be true
@@ -167,7 +168,7 @@ function parse_static_and_time(
                 getfield(out, parameter_name)[node_idx] = val
             end
         else
-            @error "$nodetype node #$node_id data not in any table."
+            @error "$nodetype node  $(repr(node_name)) (#$node_id) data not in any table."
             errors = true
         end
     end
@@ -179,10 +180,11 @@ function static_and_time_node_ids(
     static::StructVector,
     time::StructVector,
     node_type::String,
-)::Tuple{Set{Int}, Set{Int}, Vector{Int}, Bool}
+)::Tuple{Set{Int}, Set{Int}, Vector{Int}, Vector{String}, Bool}
     static_node_ids = Set(static.node_id)
     time_node_ids = Set(time.node_id)
     node_ids = get_ids(db, node_type)
+    node_names = get_names(db, node_type)
     doubles = intersect(static_node_ids, time_node_ids)
     errors = false
     if !isempty(doubles)
@@ -193,9 +195,12 @@ function static_and_time_node_ids(
         errors = true
         @error "$node_type node IDs don't match."
     end
-    return static_node_ids, time_node_ids, node_ids, !errors
+    return static_node_ids, time_node_ids, node_ids, node_names, !errors
 end
 
+const nonconservative_nodetypes =
+    Set{String}(["Basin", "LevelBoundary", "FlowBoundary", "Terminal", "User"])
+
 function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity
     if !valid_edge_types(db)
         error("Invalid edge types found.")
@@ -208,6 +213,16 @@ function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity
     edge_ids_flow_inv = Dictionary(values(edge_ids_flow), keys(edge_ids_flow))
 
     flow = adjacency_matrix(graph_flow, Float64)
+    # Add a self-loop, i.e. an entry on the diagonal, for all non-conservative node types.
+    # This is used to store the gain (positive) or loss (negative) for the water balance.
+    # Note that this only affects the sparsity structure.
+    # We want to do it here to avoid changing that during the simulation and keeping it predictable,
+    # e.g. if we wouldn't do this, inactive nodes can appear if control turns them on during runtime.
+    for (i, nodetype) in enumerate(get_nodetypes(db))
+        if nodetype in nonconservative_nodetypes
+            flow[i, i] = 1.0
+        end
+    end
     flow .= 0.0
 
     if config.solver.autodiff
@@ -253,7 +268,7 @@ function TabulatedRatingCurve(db::DB, config::Config)::TabulatedRatingCurve
     static = load_structvector(db, config, TabulatedRatingCurveStaticV1)
     time = load_structvector(db, config, TabulatedRatingCurveTimeV1)
 
-    static_node_ids, time_node_ids, node_ids, valid =
+    static_node_ids, time_node_ids, node_ids, node_names, valid =
         static_and_time_node_ids(db, static, time, "TabulatedRatingCurve")
 
     if !valid
@@ -267,7 +282,7 @@ function TabulatedRatingCurve(db::DB, config::Config)::TabulatedRatingCurve
     active = BitVector()
     errors = false
 
-    for node_id in node_ids
+    for (node_id, node_name) in zip(node_ids, node_names)
         if node_id in static_node_ids
             # Loop over all static rating curves (groups) with this node_id.
             # If it has a control_state add it to control_mapping.
@@ -298,11 +313,11 @@ function TabulatedRatingCurve(db::DB, config::Config)::TabulatedRatingCurve
             push!(interpolations, interpolation)
             push!(active, true)
         else
-            @error "TabulatedRatingCurve node #$node_id data not in any table."
+            @error "TabulatedRatingCurve node $(repr(node_name)) (#$node_id) data not in any table."
             errors = true
         end
         if !is_valid
-            @error "A Q(h) relationship for TabulatedRatingCurve #$node_id from the $source table has repeated levels, this can not be interpolated."
+            @error "A Q(h) relationship for TabulatedRatingCurve $(repr(node_name)) (#$node_id) from the $source table has repeated levels, this can not be interpolated."
             errors = true
         end
     end
@@ -353,7 +368,7 @@ function LevelBoundary(db::DB, config::Config)::LevelBoundary
     static = load_structvector(db, config, LevelBoundaryStaticV1)
     time = load_structvector(db, config, LevelBoundaryTimeV1)
 
-    static_node_ids, time_node_ids, node_ids, valid =
+    static_node_ids, time_node_ids, node_ids, node_names, valid =
         static_and_time_node_ids(db, static, time, "LevelBoundary")
 
     if !valid
@@ -381,7 +396,7 @@ function FlowBoundary(db::DB, config::Config)::FlowBoundary
     static = load_structvector(db, config, FlowBoundaryStaticV1)
     time = load_structvector(db, config, FlowBoundaryTimeV1)
 
-    static_node_ids, time_node_ids, node_ids, valid =
+    static_node_ids, time_node_ids, node_ids, node_names, valid =
         static_and_time_node_ids(db, static, time, "FlowBoundary")
 
     if !valid
@@ -401,7 +416,7 @@ function FlowBoundary(db::DB, config::Config)::FlowBoundary
     for itp in parsed_parameters.flow_rate
         if any(itp.u .< 0.0)
             @error(
-                "Currently negative flow rates are not supported, found some for dynamic flow boundary #$node_id."
+                "Currently negative flow rates are not supported, found some in dynamic flow boundary."
             )
             valid = false
         end
@@ -568,7 +583,7 @@ function PidControl(db::DB, config::Config, chunk_size::Int)::PidControl
     static = load_structvector(db, config, PidControlStaticV1)
     time = load_structvector(db, config, PidControlTimeV1)
 
-    static_node_ids, time_node_ids, node_ids, valid =
+    static_node_ids, time_node_ids, node_ids, node_names, valid =
         static_and_time_node_ids(db, static, time, "PidControl")
 
     if !valid
@@ -630,7 +645,7 @@ function User(db::DB, config::Config)::User
     static = load_structvector(db, config, UserStaticV1)
     time = load_structvector(db, config, UserTimeV1)
 
-    static_node_ids, time_node_ids, node_ids, valid =
+    static_node_ids, time_node_ids, node_ids, node_names, valid =
         static_and_time_node_ids(db, static, time, "User")
 
     if !valid
diff --git a/core/src/io.jl b/core/src/io.jl
index afe50db34..ea4be2277 100644
--- a/core/src/io.jl
+++ b/core/src/io.jl
@@ -1,3 +1,7 @@
+function get_nodetypes(db::DB)::Vector{String}
+    return only(execute(columntable, db, "SELECT type FROM Node ORDER BY fid"))
+end
+
 function get_ids(db::DB)::Vector{Int}
     return only(execute(columntable, db, "SELECT fid FROM Node ORDER BY fid"))
 end
@@ -7,6 +11,15 @@ function get_ids(db::DB, nodetype)::Vector{Int}
     return only(execute(columntable, db, sql))
 end
 
+function get_names(db::DB)::Vector{String}
+    return only(execute(columntable, db, "SELECT name FROM Node ORDER BY fid"))
+end
+
+function get_names(db::DB, nodetype)::Vector{String}
+    sql = "SELECT name FROM Node where type = $(esc_id(nodetype)) ORDER BY fid"
+    return only(execute(columntable, db, sql))
+end
+
 function exists(db::DB, tablename::String)
     query = execute(
         db,
@@ -157,59 +170,57 @@ function get_storages_and_levels(
     return (; time = tsteps, node_id, storage, level)
 end
 
-function write_basin_output(model::Model, compress)
-    (; config, integrator) = model
-    (; p) = integrator
-
+"Create the basin result table from the saved data"
+function basin_table(model::Model)::NamedTuple
     data = get_storages_and_levels(model)
     nbasin = length(data.node_id)
     ntsteps = length(data.time)
 
-    time = convert.(Arrow.DATETIME, repeat(data.time; inner = nbasin))
+    time = repeat(data.time; inner = nbasin)
     node_id = repeat(data.node_id; outer = ntsteps)
 
-    basin = (; time, node_id, storage = vec(data.storage), level = vec(data.level))
-    path = output_path(config, config.output.basin)
-    mkpath(dirname(path))
-    Arrow.write(path, basin; compress)
+    return (; time, node_id, storage = vec(data.storage), level = vec(data.level))
 end
 
-function write_flow_output(model::Model, compress)
+"Create a flow result table from the saved data"
+function flow_table(model::Model)::NamedTuple
     (; config, saved_flow, integrator) = model
     (; t, saveval) = saved_flow
     (; connectivity) = integrator.p
 
     I, J, _ = findnz(get_tmp(connectivity.flow, integrator.u))
-    unique_edge_ids = [connectivity.edge_ids_flow[ij] for ij in zip(I, J)]
+    # self-loops have no edge ID
+    unique_edge_ids = [get(connectivity.edge_ids_flow, ij, missing) for ij in zip(I, J)]
     nflow = length(I)
     ntsteps = length(t)
 
-    time =
-        convert.(
-            Arrow.DATETIME,
-            repeat(datetime_since.(t, config.starttime); inner = nflow),
-        )
-
+    time = repeat(datetime_since.(t, config.starttime); inner = nflow)
     edge_id = repeat(unique_edge_ids; outer = ntsteps)
     from_node_id = repeat(I; outer = ntsteps)
     to_node_id = repeat(J; outer = ntsteps)
     flow = FlatVector(saveval)
 
-    table = (; time, edge_id, from_node_id, to_node_id, flow)
-    path = output_path(config, config.output.flow)
-    mkpath(dirname(path))
-    Arrow.write(path, table; compress)
+    return (; time, edge_id, from_node_id, to_node_id, flow)
 end
 
-function write_discrete_control_output(model::Model, compress)
-    config = model.config
-    record = model.integrator.p.discrete_control.record
+"Create a discrete control result table from the saved data"
+function discrete_control_table(model::Model)::NamedTuple
+    (; config) = model
+    (; record) = model.integrator.p.discrete_control
 
-    time = convert.(Arrow.DATETIME, datetime_since.(record.time, config.starttime))
-
-    table = (; time, record.control_node_id, record.truth_state, record.control_state)
+    time = datetime_since.(record.time, config.starttime)
+    return (; time, record.control_node_id, record.truth_state, record.control_state)
+end
 
-    path = output_path(config, config.output.control)
+"Write a result table to disk as an Arrow file"
+function write_arrow(
+    path::AbstractString,
+    table::NamedTuple,
+    compress::TranscodingStreams.Codec,
+)
+    # ensure DateTime is encoded in a compatible manner
+    # https://github.com/apache/arrow-julia/issues/303
+    table = merge(table, (; time = convert.(Arrow.DATETIME, table.time)))
     mkpath(dirname(path))
     Arrow.write(path, table; compress)
 end
diff --git a/core/src/lib.jl b/core/src/lib.jl
index b18804bfb..ab33f6f44 100644
--- a/core/src/lib.jl
+++ b/core/src/lib.jl
@@ -29,7 +29,13 @@ function Model(config::Config)::Model
     return BMI.initialize(Model, config::Config)
 end
 
-timesteps(model::Model) = model.integrator.sol.t
+"Get all saved times in seconds since start"
+timesteps(model::Model)::Vector{Float64} = model.integrator.sol.t
+
+"Get all saved times as a Vector{DateTime}"
+function datetimes(model::Model)::Vector{DateTime}
+    return datetime_since.(timesteps(model), model.config.starttime)
+end
 
 function Base.show(io::IO, model::Model)
     (; config, integrator) = model
diff --git a/core/src/solve.jl b/core/src/solve.jl
index 440472049..640863469 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -551,12 +551,17 @@ end
 Smoothly let the evaporation flux go to 0 when at small water depths
 Currently at less than 0.1 m.
 """
-function formulate!(du::AbstractVector, basin::Basin, storage::AbstractVector)::Nothing
-    (; current_level, current_area) = basin
+function formulate_basins!(
+    du::AbstractVector,
+    basin::Basin,
+    flow::AbstractMatrix,
+    storage::AbstractVector,
+)::Nothing
+    (; node_id, current_level, current_area) = basin
     current_level = get_tmp(current_level, storage)
     current_area = get_tmp(current_area, storage)
 
-    for i in eachindex(storage)
+    for (i, id) in enumerate(node_id)
         # add all precipitation that falls within the profile
         level = current_level[i]
         area = current_area[i]
@@ -571,7 +576,9 @@ function formulate!(du::AbstractVector, basin::Basin, storage::AbstractVector)::
         drainage = basin.drainage[i]
         infiltration = factor * basin.infiltration[i]
 
-        du.storage[i] += precipitation - evaporation + drainage - infiltration
+        influx = precipitation - evaporation + drainage - infiltration
+        du.storage[i] += influx
+        flow[id, id] = influx
     end
     return nothing
 end
@@ -774,8 +781,6 @@ function formulate_flow!(
         dst_id = only(outneighbors(graph_flow, id))
 
         if !active[i]
-            flow[src_id, id] = 0.0
-            flow[id, dst_id] = 0.0
             continue
         end
 
@@ -802,6 +807,7 @@ function formulate_flow!(
 
         # Return flow is immediate
         flow[id, dst_id] = q * return_factor[i]
+        flow[id, id] = -q * (1 - return_factor[i])
     end
     return nothing
 end
@@ -831,9 +837,6 @@ function formulate_flow!(
                 ) / resistance[i]
             flow[basin_a_id, id] = q
             flow[id, basin_b_id] = q
-        else
-            flow[basin_a_id, id] = 0.0
-            flow[id, basin_b_id] = 0.0
         end
     end
     return nothing
@@ -928,8 +931,6 @@ function formulate_flow!(
         basin_b_id = only(outneighbors(graph_flow, id))
 
         if !active[i]
-            flow[basin_a_id, id] = 0.0
-            flow[id, basin_b_id] = 0.0
             continue
         end
 
@@ -990,6 +991,50 @@ function formulate_flow!(
     return nothing
 end
 
+function formulate_flow!(
+    terminal::Terminal,
+    p::Parameters,
+    storage::AbstractVector,
+    t::Float64,
+)::Nothing
+    (; connectivity) = p
+    (; graph_flow, flow) = connectivity
+    (; node_id) = terminal
+    flow = get_tmp(flow, storage)
+
+    for id in node_id
+        for upstream_id in inneighbors(graph_flow, id)
+            q = flow[upstream_id, id]
+            flow[id, id] -= q
+        end
+    end
+    return nothing
+end
+
+function formulate_flow!(
+    level_boundary::LevelBoundary,
+    p::Parameters,
+    storage::AbstractVector,
+    t::Float64,
+)::Nothing
+    (; connectivity) = p
+    (; graph_flow, flow) = connectivity
+    (; node_id) = level_boundary
+    flow = get_tmp(flow, storage)
+
+    for id in node_id
+        for in_id in inneighbors(graph_flow, id)
+            q = flow[in_id, id]
+            flow[id, id] -= q
+        end
+        for out_id in outneighbors(graph_flow, id)
+            q = flow[id, out_id]
+            flow[id, id] += q
+        end
+    end
+    return nothing
+end
+
 function formulate_flow!(
     flow_boundary::FlowBoundary,
     p::Parameters,
@@ -1005,7 +1050,6 @@ function formulate_flow!(
         # Requirement: edge points away from the flow boundary
         for dst_id in outneighbors(graph_flow, id)
             if !active[i]
-                flow[id, dst_id] = 0.0
                 continue
             end
 
@@ -1013,6 +1057,7 @@ function formulate_flow!(
 
             # Adding water is always possible
             flow[id, dst_id] = rate
+            flow[id, id] = rate
         end
     end
 end
@@ -1034,8 +1079,6 @@ function formulate_flow!(
         dst_id = only(outneighbors(graph_flow, id))
 
         if !isactive || pid_controlled
-            flow[src_id, id] = 0.0
-            flow[id, dst_id] = 0.0
             continue
         end
 
@@ -1070,8 +1113,6 @@ function formulate_flow!(
         dst_id = only(outneighbors(graph_flow, id))
 
         if !active[i] || is_pid_controlled[i]
-            flow[src_id, id] = 0.0
-            flow[id, dst_id] = 0.0
             continue
         end
 
@@ -1104,7 +1145,7 @@ function formulate_flow!(
     return nothing
 end
 
-function formulate!(
+function formulate_du!(
     du::ComponentVector,
     connectivity::Connectivity,
     flow::AbstractMatrix,
@@ -1131,10 +1172,12 @@ function formulate_flows!(p::Parameters, storage::AbstractVector, t::Float64)::N
         manning_resistance,
         tabulated_rating_curve,
         flow_boundary,
+        level_boundary,
         pump,
         outlet,
         user,
         fractional_flow,
+        terminal,
     ) = p
 
     formulate_flow!(linear_resistance, p, storage, t)
@@ -1145,8 +1188,10 @@ function formulate_flows!(p::Parameters, storage::AbstractVector, t::Float64)::N
     formulate_flow!(outlet, p, storage, t)
     formulate_flow!(user, p, storage, t)
 
-    # FractionalFlow must be done last as it relies on formulated input flows
+    # do these last since they rely on formulated input flows
     formulate_flow!(fractional_flow, p, storage, t)
+    formulate_flow!(level_boundary, p, storage, t)
+    formulate_flow!(terminal, p, storage, t)
 end
 
 """
@@ -1172,13 +1217,13 @@ function water_balance!(
     set_current_basin_properties!(basin, storage)
 
     # Basin forcings
-    formulate!(du, basin, storage)
+    formulate_basins!(du, basin, flow, storage)
 
     # First formulate intermediate flows
     formulate_flows!(p, storage, t)
 
     # Now formulate du
-    formulate!(du, connectivity, flow, basin)
+    formulate_du!(du, connectivity, flow, basin)
 
     # PID control (changes the du of PID controlled basins)
     continuous_control!(u, du, pid_control, p, integral, t)
diff --git a/core/src/validation.jl b/core/src/validation.jl
index 1dc79af3f..1cc358867 100644
--- a/core/src/validation.jl
+++ b/core/src/validation.jl
@@ -127,11 +127,13 @@ n_neighbor_bounds_control(nodetype) =
 # TODO NodeV1 and EdgeV1 are not yet used
 @version NodeV1 begin
     fid::Int
+    name::String = isnothing(s) ? "" : String(s)
     type::String = in(Symbol(type), nodetypes) ? type : error("Unknown node type $type")
 end
 
 @version EdgeV1 begin
     fid::Int
+    name::String = isnothing(s) ? "" : String(s)
     from_node_id::Int
     to_node_id::Int
     edge_type::String
diff --git a/core/test/run_models.jl b/core/test/run_models.jl
index 6a59c0cc6..116031f4c 100644
--- a/core/test/run_models.jl
+++ b/core/test/run_models.jl
@@ -1,3 +1,4 @@
+using Dates
 using Logging: Debug, with_logger
 using Test
 using Ribasim
@@ -5,6 +6,7 @@ import BasicModelInterface as BMI
 using SciMLBase: successful_retcode
 import Tables
 using PreallocationTools: get_tmp
+using DataFrames: DataFrame
 
 @testset "trivial model" begin
     toml_path = normpath(@__DIR__, "../../generated_testmodels/trivial/trivial.toml")
@@ -38,13 +40,21 @@ end
     @test all(isconcretetype, fieldtypes(typeof(p)))
 
     @test successful_retcode(model)
-    @test allunique(model.integrator.sol.t)
+    @test allunique(Ribasim.timesteps(model))
     @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip =
         Sys.isapple() atol = 1.5
 
     @test length(logger.logs) == 7
     @test logger.logs[1].level == Debug
     @test logger.logs[1].message == "Read database into memory."
+
+    table = Ribasim.flow_table(model)
+    @test Tables.schema(table) == Tables.Schema(
+        (:time, :edge_id, :from_node_id, :to_node_id, :flow),
+        (DateTime, Union{Int, Missing}, Int, Int, Float64),
+    )
+    # flows are recorded at the end of each period, and are undefined at the start
+    @test unique(table.time) == Ribasim.datetimes(model)[2:end]
 end
 
 @testset "basic transient model" begin
@@ -158,13 +168,15 @@ end
     level = level[1]
 
     timesteps = model.saved_flow.t
-    outlet_flow = [saveval[1] for saveval in model.saved_flow.saveval]
+    flow = DataFrame(Ribasim.flow_table(model))
+    outlet_flow =
+        filter([:from_node_id, :to_node_id] => (from, to) -> from === 2 && to === 3, flow)
 
     t_min_crest_level =
         level.t[2] * (outlet.min_crest_level[1] - level.u[1]) / (level.u[2] - level.u[1])
 
     # No outlet flow when upstream level is below minimum crest level
-    @test all(@. outlet_flow[timesteps <= t_min_crest_level] == 0)
+    @test all(@. outlet_flow.flow[timesteps <= t_min_crest_level] == 0)
 
     timesteps = Ribasim.timesteps(model)
     t_maximum_level = level.t[2]
@@ -258,7 +270,7 @@ end
     # slope is. See e.g.:
     # https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/latest/theoretical-basis-for-one-dimensional-and-two-dimensional-hydrodynamic-calculations/1d-steady-flow-water-surface-profiles/friction-loss-evaluation
     @test all(isapprox.(h_expected, h_actual; atol = 0.02))
-    # Test for conservation of mass
-    @test all(isapprox.(model.saved_flow.saveval[end], 5.0, atol = 0.001)) skip =
-        Sys.isapple()
+    # Test for conservation of mass, flow at the beginning == flow at the end
+    @test model.saved_flow.saveval[end][2] ≈ 5.0 atol = 0.001 skip = Sys.isapple()
+    @test model.saved_flow.saveval[end][end - 1] ≈ 5.0 atol = 0.001 skip = Sys.isapple()
 end
diff --git a/core/test/time.jl b/core/test/time.jl
index 6c5f89b3e..c04764c6b 100644
--- a/core/test/time.jl
+++ b/core/test/time.jl
@@ -1,5 +1,7 @@
 using Ribasim
 using SQLite
+using Dates
+using DataFrames: DataFrame
 
 @testset "Time dependent flow boundary" begin
     toml_path = normpath(
@@ -9,14 +11,24 @@ using SQLite
     @test ispath(toml_path)
     model = Ribasim.run(toml_path)
 
-    flow = [flows[1] for flows in model.saved_flow.saveval]
-    i_start = searchsortedlast(flow, 1)
-    i_end = searchsortedfirst(flow, 2)
-
-    t = model.saved_flow.t[i_start:i_end]
+    flow = DataFrame(Ribasim.flow_table(model))
+    # only from March to September the FlowBoundary varies
+    sin_timestamps = 3 .<= month.(unique(flow.time)) .< 10
+
+    flow_added_1 = filter(
+        [:from_node_id, :to_node_id] => (from, to) -> from === 1 && to === 1,
+        flow,
+    ).flow[sin_timestamps]
+    flow_1_to_2 = filter(
+        [:from_node_id, :to_node_id] => (from, to) -> from === 1 && to === 2,
+        flow,
+    ).flow[sin_timestamps]
+    @test flow_added_1 == flow_1_to_2
+
+    t = model.saved_flow.t[sin_timestamps]
     flow_expected = @. 1 + sin(0.5 * π * (t - t[1]) / (t[end] - t[1]))^2
-
-    @test isapprox(flow[i_start:i_end], flow_expected, rtol = 0.001)
+    # some difference is expected since the modeled flow is for the period up to t
+    @test isapprox(flow_added_1, flow_expected, rtol = 0.005)
 end
 
 @testset "User demand interpolation" begin
diff --git a/core/test/validation.jl b/core/test/validation.jl
index 2eacb4963..1f8fcb42e 100644
--- a/core/test/validation.jl
+++ b/core/test/validation.jl
@@ -43,10 +43,10 @@ end
     @test length(logger.logs) == 2
     @test logger.logs[1].level == Error
     @test logger.logs[1].message ==
-          "A Q(h) relationship for TabulatedRatingCurve #1 from the static table has repeated levels, this can not be interpolated."
+          "A Q(h) relationship for TabulatedRatingCurve \"\" (#1) from the static table has repeated levels, this can not be interpolated."
     @test logger.logs[2].level == Error
     @test logger.logs[2].message ==
-          "A Q(h) relationship for TabulatedRatingCurve #2 from the time table has repeated levels, this can not be interpolated."
+          "A Q(h) relationship for TabulatedRatingCurve \"\" (#2) from the time table has repeated levels, this can not be interpolated."
 end
 
 @testset "Neighbor count validation" begin
diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd
index d55e98f52..491743c91 100644
--- a/docs/contribute/addnode.qmd
+++ b/docs/contribute/addnode.qmd
@@ -83,6 +83,10 @@ function formulate_flow!(new_node_type::NewNodeType, p::Parameters)::Nothing
 end
 ```
 
+If the new node type is non-conservative, meaning it either adds or removes water from the model, these boundary flows also need to be recorded.
+This is done by storing it on the diagonal of the `flow[from, to]` matrix, e.g. `flow[id, id] = q`, where `q` is positive for water added to the model.
+Non-conservative node types need to be added to the `nonconservative_nodetypes` set such that this diagonal is set to a nonzero on creating the flow sparse matrix in the `Connectivity` constructor.
+
 ## The Jacobian
 
 See [Equations](../core/equations.qmd#the-jacobian) for a mathematical description of the Jacobian.
diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd
index 94135cd5a..9865146c5 100644
--- a/docs/core/usage.qmd
+++ b/docs/core/usage.qmd
@@ -286,15 +286,17 @@ The table is sorted by time, and per time it is sorted by `node_id`.
 The flow table contains outputs of the flow on every edge in the model, for each solver
 timestep.
 
-column        | type     | unit
-------------- | -------- | ----
-time          | DateTime | -
-edge_id       | Int      | -
-from_node_id  | Int      | -
-to_node_id    | Int      | -
-flow          | Float64  | $m^3 s^{-1}$
-
-The table is sorted by time, and per time the same edge_id order is used, though not sorted.
+column        | type                | unit
+------------- | ------------------- | ----
+time          | DateTime            | -
+edge_id       | Union{Int, Missing} | -
+from_node_id  | Int                 | -
+to_node_id    | Int                 | -
+flow          | Float64             | $m^3 s^{-1}$
+
+The table is sorted by time, and per time the same `edge_id` order is used, though not sorted.
+Flows that are added to the model at a node, have a missing `edge_id`, and identical `from_node_id` and `to_node_id`.
+Flows out of the model always have a negative sign, and additions a positive sign.
 
 # FractionalFlow
 
diff --git a/docs/schema/Edge.schema.json b/docs/schema/Edge.schema.json
index 319a947ba..7cea96291 100644
--- a/docs/schema/Edge.schema.json
+++ b/docs/schema/Edge.schema.json
@@ -9,6 +9,10 @@
             "format": "default",
             "type": "integer"
         },
+        "name": {
+            "format": "default",
+            "type": "string"
+        },
         "from_node_id": {
             "format": "default",
             "type": "integer"
@@ -30,6 +34,7 @@
     },
     "required": [
         "fid",
+        "name",
         "from_node_id",
         "to_node_id",
         "edge_type"
diff --git a/docs/schema/Node.schema.json b/docs/schema/Node.schema.json
index 81a47a013..18a7dc643 100644
--- a/docs/schema/Node.schema.json
+++ b/docs/schema/Node.schema.json
@@ -9,6 +9,10 @@
             "format": "default",
             "type": "integer"
         },
+        "name": {
+            "format": "default",
+            "type": "string"
+        },
         "type": {
             "format": "default",
             "type": "string"
@@ -22,6 +26,7 @@
     },
     "required": [
         "fid",
+        "name",
         "type"
     ]
 }
diff --git a/pixi.toml b/pixi.toml
index 1790f5cc1..69b4c653b 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -60,7 +60,6 @@ build-libribasim = "cd build/create_binaries && julia --project create_lib.jl"
 build = { depends_on = ["build-ribasim-cli", "build-libribasim"] }
 # Test
 test-ribasim-python = "pytest --numprocesses=auto python/ribasim/tests"
-test-ribasim-python-cov = "pytest --numprocesses=auto --cov=ribasim --cov-report=xml python/ribasim/tests"
 test-ribasim-api = "pytest --numprocesses=auto python/ribasim_api/tests"
 test-ribasim-cli = "pytest --numprocesses=auto --basetemp=build/ribasim_cli/tests/temp --junitxml=report.xml build/ribasim_cli/tests"
 test-ribasim-core = { cmd = "julia --project=core --eval 'using Pkg; Pkg.instantiate(); Pkg.test()'", depends_on = [
diff --git a/python/ribasim/pyproject.toml b/python/ribasim/pyproject.toml
index d77d40dc6..aa1e3142b 100644
--- a/python/ribasim/pyproject.toml
+++ b/python/ribasim/pyproject.toml
@@ -28,7 +28,7 @@ dependencies = [
 dynamic = ["version"]
 
 [project.optional-dependencies]
-tests = ["pytest", "ribasim_testmodels"]
+tests = ["pytest", "pytest-xdist", "pytest-cov", "ribasim_testmodels"]
 
 [tool.setuptools]
 zip-safe = true
diff --git a/python/ribasim/ribasim/geometry/edge.py b/python/ribasim/ribasim/geometry/edge.py
index 6a07ac313..d5075d331 100644
--- a/python/ribasim/ribasim/geometry/edge.py
+++ b/python/ribasim/ribasim/geometry/edge.py
@@ -18,10 +18,14 @@
 
 
 class StaticSchema(pa.SchemaModel):
+    name: Series[str] = pa.Field(default="")
     from_node_id: Series[int] = pa.Field(coerce=True)
     to_node_id: Series[int] = pa.Field(coerce=True)
     geometry: GeoSeries[Any]
 
+    class Config:
+        add_missing_columns = True
+
 
 class Edge(TableModel):
     """
diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py
index c604b781d..37145be82 100644
--- a/python/ribasim/ribasim/geometry/node.py
+++ b/python/ribasim/ribasim/geometry/node.py
@@ -16,9 +16,13 @@
 
 
 class StaticSchema(pa.SchemaModel):
+    name: Series[str] = pa.Field(default="")
     type: Series[str]
     geometry: GeoSeries[Any]
 
+    class Config:
+        add_missing_columns = True
+
 
 class Node(TableModel):
     """
diff --git a/python/ribasim/ribasim/models.py b/python/ribasim/ribasim/models.py
index e68def42d..eba64fd86 100644
--- a/python/ribasim/ribasim/models.py
+++ b/python/ribasim/ribasim/models.py
@@ -61,6 +61,7 @@ class DiscreteControlLogic(BaseModel):
 
 class Edge(BaseModel):
     fid: int
+    name: str
     from_node_id: int
     to_node_id: int
     edge_type: str
@@ -123,6 +124,7 @@ class ManningResistanceStatic(BaseModel):
 
 class Node(BaseModel):
     fid: int
+    name: str
     type: str
     remarks: str = Field("", description="a hack for pandera")
 
diff --git a/python/ribasim/ribasim/utils.py b/python/ribasim/ribasim/utils.py
index 6a3db761e..259bfab87 100644
--- a/python/ribasim/ribasim/utils.py
+++ b/python/ribasim/ribasim/utils.py
@@ -1,3 +1,5 @@
+import random
+import string
 from typing import Any, Sequence, Tuple
 
 import numpy as np
@@ -76,3 +78,8 @@ def connectivity_from_geometry(
     from_id = node_index[edge_node_id[:, 0]].to_numpy()
     to_id = node_index[edge_node_id[:, 1]].to_numpy()
     return from_id, to_id
+
+
+def random_string(length=3):
+    letters = string.ascii_lowercase
+    return "".join(random.choice(letters) for i in range(length))
diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py
index fbbc9967f..1ff43127f 100644
--- a/python/ribasim_testmodels/ribasim_testmodels/basic.py
+++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py
@@ -159,7 +159,10 @@ def basic_model() -> ribasim.Model:
     # Make sure the feature id starts at 1: explicitly give an index.
     node = ribasim.Node(
         static=gpd.GeoDataFrame(
-            data={"type": node_type},
+            data={
+                "type": node_type,
+                "name": [ribasim.utils.random_string() for _ in range(len(node_id))],
+            },
             index=pd.Index(node_id, name="fid"),
             geometry=node_xy,
             crs="EPSG:28992",
@@ -177,6 +180,7 @@ def basic_model() -> ribasim.Model:
     edge = ribasim.Edge(
         static=gpd.GeoDataFrame(
             data={
+                "name": [ribasim.utils.random_string() for _ in range(len(from_id))],
                 "from_node_id": from_id,
                 "to_node_id": to_id,
                 "edge_type": len(from_id) * ["flow"],
diff --git a/qgis/core/nodes.py b/qgis/core/nodes.py
index d02115dbc..93f82a0a9 100644
--- a/qgis/core/nodes.py
+++ b/qgis/core/nodes.py
@@ -28,6 +28,7 @@
 from ribasim_qgis.core import geopackage
 
 from qgis.core import (
+    Qgis,
     QgsCategorizedSymbolRenderer,
     QgsEditorWidgetSetup,
     QgsField,
@@ -128,7 +129,10 @@ def set_editor_widget(self) -> None:
 class Node(Input):
     input_type = "Node"
     geometry_type = "Point"
-    attributes = (QgsField("type", QVariant.String),)
+    attributes = (
+        QgsField("name", QVariant.String),
+        QgsField("type", QVariant.String),
+    )
 
     @classmethod
     def is_spatial(cls):
@@ -206,7 +210,8 @@ def renderer(self) -> QgsCategorizedSymbolRenderer:
     @property
     def labels(self) -> Any:
         pal_layer = QgsPalLayerSettings()
-        pal_layer.fieldName = "fid"
+        pal_layer.fieldName = """concat("name", ' (#', "fid", ')')"""
+        pal_layer.isExpression = True
         pal_layer.enabled = True
         pal_layer.dist = 2.0
         labels = QgsVectorLayerSimpleLabeling(pal_layer)
@@ -217,6 +222,7 @@ class Edge(Input):
     input_type = "Edge"
     geometry_type = "Linestring"
     attributes = [
+        QgsField("name", QVariant.String),
         QgsField("from_node_id", QVariant.Int),
         QgsField("to_node_id", QVariant.Int),
         QgsField("edge_type", QVariant.String),
@@ -281,6 +287,15 @@ def renderer(self) -> QgsCategorizedSymbolRenderer:
         )
         return renderer
 
+    @property
+    def labels(self) -> Any:
+        pal_layer = QgsPalLayerSettings()
+        pal_layer.fieldName = "name"
+        pal_layer.enabled = True
+        pal_layer.placement = Qgis.LabelPlacement.Line
+        labels = QgsVectorLayerSimpleLabeling(pal_layer)
+        return labels
+
 
 class BasinProfile(Input):
     input_type = "Basin / profile"