diff --git a/.github/workflows/core_tests.yml b/.github/workflows/core_tests.yml
index 8b3367920..04822cfae 100644
--- a/.github/workflows/core_tests.yml
+++ b/.github/workflows/core_tests.yml
@@ -28,7 +28,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Setup Micromamba
         uses: mamba-org/setup-micromamba@v1
         with:
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index e4dc493ed..26a58be4e 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -24,7 +24,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - uses: julia-actions/setup-julia@v1
         with:
diff --git a/.github/workflows/pre-commit_auto_update.yml b/.github/workflows/pre-commit_auto_update.yml
index bf4d418ae..bfe8f0621 100644
--- a/.github/workflows/pre-commit_auto_update.yml
+++ b/.github/workflows/pre-commit_auto_update.yml
@@ -11,7 +11,7 @@ jobs:
   auto-update:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
         with:
           ssh-key: ${{ secrets.SSH_PRIVATE_KEY }}
       - uses: actions/setup-python@v4
diff --git a/.github/workflows/pre-commit_check.yml b/.github/workflows/pre-commit_check.yml
index fbad75fd6..2876b910a 100644
--- a/.github/workflows/pre-commit_check.yml
+++ b/.github/workflows/pre-commit_check.yml
@@ -10,7 +10,7 @@ jobs:
   check:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - uses: actions/setup-python@v4
         with:
           python-version: "3.11"
diff --git a/.github/workflows/python_lint.yml b/.github/workflows/python_lint.yml
index de6c69f28..859391309 100644
--- a/.github/workflows/python_lint.yml
+++ b/.github/workflows/python_lint.yml
@@ -18,7 +18,7 @@ jobs:
       run:
         shell: bash -l {0}
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
       - name: Setup Micromamba
         uses: mamba-org/setup-micromamba@v1
         with:
diff --git a/.github/workflows/python_tests.yml b/.github/workflows/python_tests.yml
index fcf173ba8..1400ea31b 100644
--- a/.github/workflows/python_tests.yml
+++ b/.github/workflows/python_tests.yml
@@ -30,7 +30,7 @@ jobs:
         arch:
           - x86
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Setup Micromamba
         uses: mamba-org/setup-micromamba@v1
diff --git a/.github/workflows/qgis.yml b/.github/workflows/qgis.yml
index 6863ecd4e..788a7e059 100644
--- a/.github/workflows/qgis.yml
+++ b/.github/workflows/qgis.yml
@@ -17,7 +17,7 @@ jobs:
         steps:
 
         - name: Check out repository
-          uses: actions/checkout@v3
+          uses: actions/checkout@v4
         - name: Launching docker compose
           run: ./start.sh
         - name: Running tests
diff --git a/.teamcity/Ribasim_Ribasim/buildTypes/BuildRibasimCliWindows.xml b/.teamcity/Ribasim_Ribasim/buildTypes/BuildRibasimCliWindows.xml
index 5c080c704..3e8b67cb4 100644
--- a/.teamcity/Ribasim_Ribasim/buildTypes/BuildRibasimCliWindows.xml
+++ b/.teamcity/Ribasim_Ribasim/buildTypes/BuildRibasimCliWindows.xml
@@ -9,7 +9,7 @@
       <option name="publishArtifactCondition" value="SUCCESSFUL" />
     </options>
     <parameters>
-      <param name="conda_env_path" value="%system.teamcity.build.checkoutDir%/env" />
+      <param name="conda_env_path" value="%system.teamcity.build.checkoutDir%/env_build_ribasim_cli" />
       <param name="julia" value="%env.JULIA_1_9%" />
     </parameters>
     <build-runners>
diff --git a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_RibasimJl_BuildLibribasimWindows.xml b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_RibasimJl_BuildLibribasimWindows.xml
index 23852a530..5bdfb285c 100644
--- a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_RibasimJl_BuildLibribasimWindows.xml
+++ b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_RibasimJl_BuildLibribasimWindows.xml
@@ -12,7 +12,7 @@
       <setting-ref ref="RQ_271" />
     </disabled-settings>
     <parameters>
-      <param name="conda_env_path" value="%system.teamcity.build.checkoutDir%/env" />
+      <param name="conda_env_path" value="%system.teamcity.build.checkoutDir%/env_build_libribasim" />
       <param name="julia" value="%env.JULIA_1_9%" />
     </parameters>
     <build-runners>
diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl
index c1bc06678..8821ef2c8 100644
--- a/core/src/Ribasim.jl
+++ b/core/src/Ribasim.jl
@@ -12,7 +12,7 @@ using ComponentArrays: ComponentVector
 using DataInterpolations: LinearInterpolation, derivative
 using Dates
 using DBInterface: execute, prepare
-using Dictionaries: Indices, Dictionary, gettoken, gettokenvalue, dictionary
+using Dictionaries: Indices, Dictionary, gettoken, dictionary
 using ForwardDiff: pickchunksize
 using DiffEqCallbacks
 using Graphs: DiGraph, add_edge!, adjacency_matrix, inneighbors, outneighbors
@@ -20,7 +20,7 @@ using Legolas: Legolas, @schema, @version, validate, SchemaVersion, declared
 using Logging: current_logger, min_enabled_level, with_logger
 using LoggingExtras: EarlyFilteredLogger, LevelOverrideLogger
 using OrdinaryDiffEq
-using PreallocationTools: DiffCache, get_tmp
+using PreallocationTools: DiffCache, FixedSizeDiffCache, get_tmp
 using SciMLBase
 using SparseArrays
 using SQLite: SQLite, DB, Query, esc_id
diff --git a/core/src/create.jl b/core/src/create.jl
index 8221f4504..c75a6c0e8 100644
--- a/core/src/create.jl
+++ b/core/src/create.jl
@@ -202,10 +202,11 @@ 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)
-    nonzeros(flow) .= 0.0
+    flow .= 0.0
 
     if config.solver.autodiff
-        flow = DiffCache(flow, chunk_size)
+        # FixedSizeDiffCache performs better for sparse matrix
+        flow = FixedSizeDiffCache(flow, chunk_size)
     end
 
     return Connectivity(
diff --git a/core/src/solve.jl b/core/src/solve.jl
index 9c684d35f..227ced273 100644
--- a/core/src/solve.jl
+++ b/core/src/solve.jl
@@ -401,21 +401,21 @@ struct PidControl{T} <: AbstractParameterNode
 end
 
 # TODO Automatically add all nodetypes here
-struct Parameters
+struct Parameters{T, TSparse, C1, C2}
     starttime::DateTime
-    connectivity::Connectivity
-    basin::Basin
+    connectivity::Connectivity{TSparse}
+    basin::Basin{T, C1}
     linear_resistance::LinearResistance
     manning_resistance::ManningResistance
-    tabulated_rating_curve::TabulatedRatingCurve
+    tabulated_rating_curve::TabulatedRatingCurve{C2}
     fractional_flow::FractionalFlow
     level_boundary::LevelBoundary
     flow_boundary::FlowBoundary
-    pump::Pump
-    outlet::Outlet
+    pump::Pump{T}
+    outlet::Outlet{T}
     terminal::Terminal
     discrete_control::DiscreteControl
-    pid_control::PidControl
+    pid_control::PidControl{T}
     lookup::Dict{Int, Symbol}
 end
 
@@ -508,8 +508,8 @@ function formulate!(
     du::AbstractVector,
     basin::Basin,
     storage::AbstractVector,
-    current_area,
-    current_level,
+    current_area::AbstractVector,
+    current_level::AbstractVector,
     t::Float64,
 )::Nothing
     for i in eachindex(storage)
@@ -553,13 +553,13 @@ end
 function continuous_control!(
     u::ComponentVector,
     du::ComponentVector,
-    current_area,
+    current_area::AbstractVector,
     pid_control::PidControl,
     p::Parameters,
     integral_value::SubArray,
     current_level::AbstractVector,
-    flow::SparseMatrixCSC,
-    pid_error,
+    flow::AbstractMatrix,
+    pid_error::AbstractVector,
     t::Float64,
 )::Nothing
     (; connectivity, pump, outlet, basin, fractional_flow) = p
@@ -700,11 +700,11 @@ end
 """
 Directed graph: outflow is positive!
 """
-function formulate!(
+function formulate_flow!(
     linear_resistance::LinearResistance,
     p::Parameters,
-    current_level,
-    flow,
+    current_level::AbstractVector,
+    flow::AbstractMatrix,
     t::Float64,
 )::Nothing
     (; connectivity) = p
@@ -733,11 +733,11 @@ end
 """
 Directed graph: outflow is positive!
 """
-function formulate!(
+function formulate_flow!(
     tabulated_rating_curve::TabulatedRatingCurve,
     p::Parameters,
-    current_level,
-    flow::SparseMatrixCSC,
+    current_level::AbstractVector,
+    flow::AbstractMatrix,
     t::Float64,
 )::Nothing
     (; connectivity) = p
@@ -800,11 +800,11 @@ The average of the upstream and downstream water depth is used to compute cross-
 hydraulic radius. This ensures that a basin can receive water after it has gone
 dry.
 """
-function formulate!(
+function formulate_flow!(
     manning_resistance::ManningResistance,
     p::Parameters,
-    current_level,
-    flow,
+    current_level::AbstractVector,
+    flow::AbstractMatrix,
     t::Float64,
 )::Nothing
     (; basin, connectivity) = p
@@ -859,7 +859,11 @@ function formulate!(
     return nothing
 end
 
-function formulate!(fractional_flow::FractionalFlow, flow, p::Parameters)::Nothing
+function formulate_flow!(
+    fractional_flow::FractionalFlow,
+    flow::AbstractMatrix,
+    p::Parameters,
+)::Nothing
     (; connectivity) = p
     (; graph_flow) = connectivity
     (; node_id, fraction) = fractional_flow
@@ -871,7 +875,12 @@ function formulate!(fractional_flow::FractionalFlow, flow, p::Parameters)::Nothi
     return nothing
 end
 
-function formulate!(flow_boundary::FlowBoundary, p::Parameters, flow, t::Float64)::Nothing
+function formulate_flow!(
+    flow_boundary::FlowBoundary,
+    p::Parameters,
+    flow::AbstractMatrix,
+    t::Float64,
+)::Nothing
     (; connectivity) = p
     (; graph_flow) = connectivity
     (; node_id, active, flow_rate) = flow_boundary
@@ -892,10 +901,10 @@ function formulate!(flow_boundary::FlowBoundary, p::Parameters, flow, t::Float64
     end
 end
 
-function formulate!(
+function formulate_flow!(
     node::Union{Pump, Outlet},
     p::Parameters,
-    flow,
+    flow::AbstractMatrix,
     storage::AbstractVector,
 )::Nothing
     (; connectivity, basin) = p
@@ -938,7 +947,7 @@ end
 function formulate!(
     du::ComponentVector,
     connectivity::Connectivity,
-    flow::SparseMatrixCSC,
+    flow::AbstractMatrix,
     basin::Basin,
 )::Nothing
     # loop over basins
@@ -959,8 +968,8 @@ end
 function formulate_flows!(
     p::Parameters,
     storage::AbstractVector,
-    current_level,
-    flow::SparseMatrixCSC,
+    current_level::AbstractVector,
+    flow::AbstractMatrix,
     t::Float64,
 )::Nothing
     (;
@@ -973,13 +982,13 @@ function formulate_flows!(
         outlet,
     ) = p
 
-    formulate!(linear_resistance, p, current_level, flow, t)
-    formulate!(manning_resistance, p, current_level, flow, t)
-    formulate!(tabulated_rating_curve, p, current_level, flow, t)
-    formulate!(flow_boundary, p, flow, t)
-    formulate!(fractional_flow, flow, p)
-    formulate!(pump, p, flow, storage)
-    formulate!(outlet, p, flow, storage)
+    formulate_flow!(linear_resistance, p, current_level, flow, t)
+    formulate_flow!(manning_resistance, p, current_level, flow, t)
+    formulate_flow!(tabulated_rating_curve, p, current_level, flow, t)
+    formulate_flow!(flow_boundary, p, flow, t)
+    formulate_flow!(fractional_flow, flow, p)
+    formulate_flow!(pump, p, flow, storage)
+    formulate_flow!(outlet, p, flow, storage)
 
     return nothing
 end
@@ -1000,7 +1009,8 @@ function water_balance!(
 
     du .= 0.0
     flow = get_tmp(connectivity.flow, u)
-    nonzeros(flow) .= 0.0
+    # use parent to avoid materializing the ReinterpretArray from FixedSizeDiffCache
+    parent(flow) .= 0.0
 
     current_area = get_tmp(basin.current_area, u)
     current_level = get_tmp(basin.current_level, u)
diff --git a/core/src/utils.jl b/core/src/utils.jl
index 13da0e4b6..d4b80baaa 100644
--- a/core/src/utils.jl
+++ b/core/src/utils.jl
@@ -406,34 +406,36 @@ end
 Get the current water level of a node ID.
 The ID can belong to either a Basin or a LevelBoundary.
 """
-function get_level(p::Parameters, node_id::Int, current_level, t::Float64)::Real
+function get_level(
+    p::Parameters,
+    node_id::Int,
+    current_level::AbstractVector,
+    t::Float64,
+)::Real
     (; basin, level_boundary) = p
-    # since the node_id fields are already Indices, Dictionary creation is instant
-    basin = Dictionary(basin.node_id, current_level)
-    hasindex, token = gettoken(basin, node_id)
+    hasindex, i = id_index(basin.node_id, node_id)
     return if hasindex
-        gettokenvalue(basin, token)
+        current_level[i]
     else
-        boundary = Dictionary(level_boundary.node_id, level_boundary.level)
-        boundary[node_id](t)
+        i = findsorted(level_boundary.node_id, node_id)::Int
+        level_boundary.level[i](t)
     end
 end
 
 "Get the index of an ID in a set of indices."
-function id_index(ids::Indices{Int}, id::Int)
-    # There might be a better approach for this, this feels too internal
-    # the second return is the token, a Tuple{Int, Int}
-    hasindex, (_, idx) = gettoken(ids, id)
-    return hasindex, idx
+function id_index(ids::Indices{Int}, id::Int)::Tuple{Bool, Int}
+    # We avoid creating Dictionary here since it converts the values to a Vector,
+    # leading to allocations when used with PreallocationTools's ReinterpretArrays.
+    hasindex, (_, i) = gettoken(ids, id)
+    return hasindex, i
 end
 
 "Return the bottom elevation of the basin with index i, or nothing if it doesn't exist"
 function basin_bottom(basin::Basin, node_id::Int)::Union{Float64, Nothing}
-    basin = Dictionary(basin.node_id, basin.level)
-    hasindex, token = gettoken(basin, node_id)
+    hasindex, i = id_index(basin.node_id, node_id)
     return if hasindex
         # get level(storage) interpolation function
-        level_discrete = gettokenvalue(basin, token)
+        level_discrete = basin.level[i]
         # and return the first level in this vector, representing the bottom
         first(level_discrete)
     else
diff --git a/core/test/run_models.jl b/core/test/run_models.jl
index e95197086..466bbe2bc 100644
--- a/core/test/run_models.jl
+++ b/core/test/run_models.jl
@@ -32,6 +32,11 @@ end
     end
 
     @test model isa Ribasim.Model
+    p = model.integrator.p
+    @test p isa Ribasim.Parameters
+    @test isconcretetype(typeof(p))
+    @test all(isconcretetype, fieldtypes(typeof(p)))
+
     @test successful_retcode(model)
     @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip =
         Sys.isapple() atol = 1.5
diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd
index b925824ce..3c15e841f 100644
--- a/docs/contribute/addnode.qmd
+++ b/docs/contribute/addnode.qmd
@@ -61,10 +61,10 @@ end
 
 ## Node behavior
 
-In general if the new node type dictates flow, the behaviour of the new node in the Ribasim core is defined in a method of the `formulate!` function, which is called within the `water_balance!` (both in `solve.jl`) function being the right hand side of the system of differential equations solved by Ribasim. Here the details depend highly on the specifics of the node type. An example structure of a `formulate!` method is given below.
+In general if the new node type dictates flow, the behaviour of the new node in the Ribasim core is defined in a method of the `formulate_flow!` function, which is called within the `water_balance!` (both in `solve.jl`) function being the right hand side of the system of differential equations solved by Ribasim. Here the details depend highly on the specifics of the node type. An example structure of a `formulate_flow!` method is given below.
 
 ```julia
-function formulate!(new_node_type::NewNodeType, p::Parameters)::Nothing
+function formulate_flow!(new_node_type::NewNodeType, p::Parameters)::Nothing
     # Retrieve relevant parameters
     (; connectivity) = p
     (; flow) = connectivity