diff --git a/Manifest.toml b/Manifest.toml
index a72a3f92b..7b5174c1c 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.10.5"
 manifest_format = "2.0"
-project_hash = "acf9296ffbaffd3f335bdb6005badef66b3e10fe"
+project_hash = "986b47155507036c06218d23750fe8de28eb001a"
 
 [[deps.AbstractTrees]]
 git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
@@ -11,9 +11,9 @@ version = "0.4.5"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
+git-tree-sha1 = "d80af0733c99ea80575f612813fa6aa71022d33a"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.4"
+version = "4.1.0"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -78,15 +78,15 @@ version = "0.1.6"
 
 [[deps.Blosc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "19b98ee7e3db3b4eff74c5c9c72bf32144e24f10"
+git-tree-sha1 = "ef12cdd1c7fb7e1dfd6fa8fd60d4db6bc61d2f23"
 uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
-version = "1.21.5+0"
+version = "1.21.6+0"
 
 [[deps.Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd"
+git-tree-sha1 = "8873e196c2eb87962a2048b3b8e08946535864a1"
 uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
-version = "1.0.8+1"
+version = "1.0.8+2"
 
 [[deps.CFTime]]
 deps = ["Dates", "Printf"]
@@ -239,9 +239,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "f389674c99bfcde17dc57454011aa44d5a260a40"
+git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.6.0"
+version = "1.6.1"
 
 [[deps.LRUCache]]
 git-tree-sha1 = "b3cc6698599b10e652832c2f23db3cab99d51b59"
@@ -337,9 +337,9 @@ version = "1.10.0+1"
 
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "19d4bd098928a3263693991500d05d74dbdc2004"
+git-tree-sha1 = "7715e65c47ba3941c502bffb7f266a41a7f54423"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.2.2+0"
+version = "4.2.3+0"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
@@ -349,9 +349,9 @@ version = "0.1.11"
 
 [[deps.MPItrampoline_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "fde81c9f9c94fe5fbeaed7b3f1330305cf9a327c"
+git-tree-sha1 = "70e830dab5d0775183c99fc75e4c24c614ed7142"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
-version = "5.5.0+0"
+version = "5.5.1+0"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
@@ -559,9 +559,9 @@ weakdeps = ["OffsetArrays", "StaticArrays"]
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
+git-tree-sha1 = "777657803913ffc7e8cc20f0fd04b634f871af8f"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.7"
+version = "1.9.8"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -668,9 +668,9 @@ version = "5.11.0+0"
 
 [[deps.libzip_jll]]
 deps = ["Artifacts", "Bzip2_jll", "GnuTLS_jll", "JLLWrappers", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "3282b7d16ae7ac3e57ec2f3fa8fafb564d8f9f7f"
+git-tree-sha1 = "668ac0297e6bd8f4d53dfdcd3ace71f2e00f4a35"
 uuid = "337d8026-41b4-5cde-a456-74a10e5b31d1"
-version = "1.10.1+0"
+version = "1.11.1+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
diff --git a/Project.toml b/Project.toml
index bbbf9ec42..4bd017e4e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -18,12 +18,10 @@ NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 Aqua = "0.8"
@@ -50,7 +48,6 @@ Statistics = "1"
 TOML = "1"
 TerminalLoggers = "0.1.5"
 Test = "1"
-UnPack = "1"
 julia = "1.6"
 
 [extras]
@@ -58,7 +55,8 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Downloads", "Polynomials", "QuadGK", "Test"]
+test = ["Aqua", "Downloads", "Polynomials", "QuadGK", "Test", "Random"]
diff --git a/src/Wflow.jl b/src/Wflow.jl
index f4a8bd8f5..f4f898e0d 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -1,26 +1,25 @@
 module Wflow
 
-using Dates
-using TOML
-using Graphs
-using NCDatasets
-using StaticArrays
-using Statistics
-using UnPack
-using Random
-using BasicModelInterface
-using FieldMetadata
-using Parameters
-using DelimitedFiles
-using ProgressLogging
+import BasicModelInterface as BMI
+
+using Base.Threads: nthreads
+using CFTime: CFTime, monthday, dayofyear
+using Dates: Dates, Second, Minute, Hour, Day, Month, year, TimeType, DatePeriod, TimePeriod, Date, DateTime, now, isleapyear, datetime2unix
+using DelimitedFiles: readdlm
+using FieldMetadata: @metadata
+using Glob: glob
+using Graphs: Graphs, Graph, DiGraph, add_edge!, is_cyclic, inneighbors, outneighbors, edges, topological_sort_by_dfs, src, dst, vertices, nv, ne, induced_subgraph, add_vertex!
+using IfElse: IfElse
 using LoggingExtras
+using LoopVectorization: @tturbo
+using NCDatasets: NCDatasets, NCDataset, dimnames, dimsize, nomissing, defDim, defVar
+using Parameters: @with_kw
+using Polyester: @batch
+using ProgressLogging: @progress
+using StaticArrays: SVector, pushfirst, setindex
+using Statistics: mean, median, quantile!
 using TerminalLoggers
-using CFTime
-using Base.Threads
-using Glob
-using Polyester
-using LoopVectorization
-using IfElse
+using TOML: TOML
 
 @metadata get_units "mm dt-1" String
 # metadata for BMI grid
@@ -28,7 +27,6 @@ using IfElse
 @metadata grid_type "unstructured" String
 @metadata grid_location "node" String
 
-const BMI = BasicModelInterface
 const Float = Float64
 const CFDataset = Union{NCDataset,NCDatasets.MFDataset}
 const CFVariable_MF = Union{NCDatasets.CFVariable,NCDatasets.MFCFVariable}
@@ -213,7 +211,7 @@ function run_timestep(model::Model; update_func = update, write_model_output = t
 end
 
 function run(model::Model; close_files = true)
-    @unpack network, config, writer, clock = model
+    (; config, writer, clock) = model
 
     model_type = config.model.type::String
 
diff --git a/src/bmi.jl b/src/bmi.jl
index a2a7b9899..68093f03a 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -48,7 +48,6 @@ Update the model for a single timestep.
 - `run = nothing`: to update a model partially.
 """
 function BMI.update(model::Model; run = nothing)
-    @unpack clock, network, config = model
     if isnothing(run)
         model = run_timestep(model)
     elseif run == "sbm_until_recharge"
@@ -64,7 +63,6 @@ function BMI.update(model::Model; run = nothing)
 end
 
 function BMI.update_until(model::Model, time::Float64)
-    @unpack clock, network, config = model
     t = BMI.get_current_time(model)
     _div, _rem = divrem(time - t, model.clock.dt.value)
     steps = Int(_div)
@@ -85,7 +83,7 @@ end
 
 "Write state output to netCDF and close files."
 function BMI.finalize(model::Model)
-    @unpack config, writer, clock = model
+    (; config, writer) = model
     # it is possible that the state dataset has been closed by `save_state`
     if !isnothing(writer.state_dataset) && isopen(writer.state_dataset)
         write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
@@ -95,8 +93,7 @@ function BMI.finalize(model::Model)
 end
 
 function BMI.get_component_name(model::Model)
-    @unpack config = model
-    return config.model.type
+    return model.config.model.type
 end
 
 function BMI.get_input_item_count(model::Model)
@@ -115,7 +112,7 @@ This `API` sections contains a list of `Model` components for which variables ca
 exchanged.
 """
 function BMI.get_input_var_names(model::Model)
-    @unpack config = model
+    (; config) = model
     if haskey(config, "API")
         var_names = Vector{String}()
         for c in config.API.components
@@ -211,7 +208,7 @@ function BMI.get_var_location(model::Model, name::String)
 end
 
 function BMI.get_current_time(model::Model)
-    @unpack config = model
+    (; config) = model
     calendar = get(config, "calendar", "standard")::String
     starttime = cftime(config.starttime, calendar)
     return 0.001 * Dates.value(model.clock.time - starttime)
@@ -222,7 +219,7 @@ function BMI.get_start_time(model::Model)
 end
 
 function BMI.get_end_time(model::Model)
-    @unpack config = model
+    (; config) = model
     calendar = get(config, "calendar", "standard")::String
     starttime = cftime(config.starttime, calendar)
     endtime = cftime(config.endtime, calendar)
@@ -243,7 +240,7 @@ function BMI.get_value(model::Model, name::String, dest::Vector{T}) where {T<:Ab
 end
 
 function BMI.get_value_ptr(model::Model, name::String)
-    @unpack network = model
+    (; network) = model
     s = split(name, "[")
     key = symbols(first(s))
     if exchange(param(model, key[1:end-1]), key[end]) == 1
@@ -322,9 +319,9 @@ function BMI.get_grid_rank(model::Model, grid::Int)
 end
 
 function BMI.get_grid_x(model::Model, grid::Int, x::Vector{T}) where {T<:AbstractFloat}
-    @unpack reader, config = model
-    @unpack dataset = reader
-    sel = active_indices(model.network, grids[grid])
+    (; reader, network) = model
+    (; dataset) = reader
+    sel = active_indices(network, grids[grid])
     inds = [sel[i][1] for i in eachindex(sel)]
     x_nc = read_x_axis(dataset)
     x .= x_nc[inds]
@@ -332,9 +329,9 @@ function BMI.get_grid_x(model::Model, grid::Int, x::Vector{T}) where {T<:Abstrac
 end
 
 function BMI.get_grid_y(model::Model, grid::Int, y::Vector{T}) where {T<:AbstractFloat}
-    @unpack reader, config = model
-    @unpack dataset = reader
-    sel = active_indices(model.network, grids[grid])
+    (; reader, network) = model
+    (; dataset) = reader
+    sel = active_indices(network, grids[grid])
     inds = [sel[i][2] for i in eachindex(sel)]
     y_nc = read_y_axis(dataset)
     y .= y_nc[inds]
@@ -350,7 +347,7 @@ function BMI.get_grid_size(model::Model, grid::Int)
 end
 
 function BMI.get_grid_edge_count(model::Model, grid::Int)
-    @unpack network = model
+    (; network) = model
     if grid == 3
         return ne(network.river.graph)
     elseif grid == 4
@@ -365,7 +362,7 @@ function BMI.get_grid_edge_count(model::Model, grid::Int)
 end
 
 function BMI.get_grid_edge_nodes(model::Model, grid::Int, edge_nodes::Vector{Int})
-    @unpack network = model
+    (; network) = model
     n = length(edge_nodes)
     m = div(n, 2)
     # inactive nodes (boundary/ghost points) are set at -999
@@ -402,7 +399,7 @@ function load_state(model::Model)
 end
 
 function save_state(model::Model)
-    @unpack config, writer, clock = model
+    (; config, writer) = model
     if haskey(config, "state") && haskey(config.state, "path_output")
         @info "Write output states to netCDF file `$(model.writer.state_nc_path)`."
     end
diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
index 039759940..8ea773d63 100644
--- a/src/flextopo_model.jl
+++ b/src/flextopo_model.jl
@@ -665,7 +665,7 @@ function initialize_flextopo_model(config::Config)
 end
 
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network, config) = model
 
     inds_riv = network.index_river
 
@@ -716,7 +716,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
 end
 
 function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
-    @unpack lateral, config = model
+    (; lateral, config) = model
     reinit = get(config.model, "reinit", true)::Bool
     # read and set states in model object if reinit=true
     if reinit == false
diff --git a/src/flow.jl b/src/flow.jl
index c778b9f6c..f7578cc03 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -189,7 +189,7 @@ end
 
 
 function update(sf::SurfaceFlowLand, network, frac_toriver)
-    @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
+    (; subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes) =
         network
 
     ns = length(subdomain_order)
@@ -256,7 +256,7 @@ end
 
 function update(sf::SurfaceFlowRiver, network, doy)
 
-    @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
+    (; graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes) =
         network
 
     ns = length(subdomain_order)
@@ -429,7 +429,7 @@ end
 end
 
 function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
-    @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes, area =
+    (; subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes, area) =
         network
 
 
@@ -749,7 +749,7 @@ end
 
 function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, update_h)
 
-    @unpack nodes_at_link, links_at_node = network
+    (; nodes_at_link, links_at_node) = network
 
     sw.q0 .= sw.q
     if !isnothing(sw.floodplain)
@@ -952,8 +952,6 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
 end
 
 function update(sw::ShallowWaterRiver{T}, network, doy; update_h = true) where {T}
-    @unpack nodes_at_link, links_at_node = network
-
     if !isnothing(sw.reservoir)
         sw.reservoir.inflow .= 0.0
         sw.reservoir.totaloutflow .= 0.0
@@ -1192,9 +1190,6 @@ function update(
     doy;
     update_h = false,
 ) where {T}
-
-    @unpack nodes_at_link, links_at_node = network.river
-
     if !isnothing(swr.reservoir)
         swr.reservoir.inflow .= 0.0
         swr.reservoir.totaloutflow .= 0.0
@@ -1236,7 +1231,7 @@ function shallowwater_update(
     indices = network.land.staggered_indices
     inds_riv = network.land.index_river
 
-    @unpack nodes_at_link, links_at_node = network.river
+    (; links_at_node) = network.river
 
     sw.qx0 .= sw.qx
     sw.qy0 .= sw.qy
@@ -1637,7 +1632,7 @@ function set_river_inwater(
     model::Model{N,L,V,R,W,T},
     ssf_toriver,
 ) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
-    @unpack lateral, vertical, network, config = model
+    (; lateral, vertical, network, config) = model
     inds = network.index_river
     do_water_demand = haskey(config.model, "water_demand")
     if do_water_demand
@@ -1667,7 +1662,7 @@ end
 Set `inwater` of the lateral river component (based on overland flow).
 """
 function set_river_inwater(model, ssf_toriver)
-    @unpack lateral, network = model
+    (; lateral, network) = model
     inds = network.index_river
     lateral.river.inwater .= lateral.land.to_river[inds]
 end
@@ -1678,7 +1673,7 @@ end
 Set `inwater` of the lateral land component for the `SbmGwfModel` type.
 """
 function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
-    @unpack lateral, vertical, network, config = model
+    (; lateral, vertical, network, config) = model
     do_water_demand = haskey(config.model, "water_demand")
     do_drains = get(config.model, "drains", false)::Bool
     drainflux = zeros(vertical.n)
@@ -1703,7 +1698,7 @@ end
 Set `inwater` of the lateral land component for the `SbmModel` type.
 """
 function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
-    @unpack lateral, vertical, network, config = model
+    (; lateral, vertical, network, config) = model
     do_water_demand = haskey(config.model, "water_demand")
     if do_water_demand
         @. lateral.land.inwater =
@@ -1722,7 +1717,7 @@ end
 Set `inwater` of the lateral land component, based on `runoff` of the `vertical` concept.
 """
 function set_land_inwater(model)
-    @unpack lateral, vertical, network = model
+    (; lateral, vertical, network) = model
     @. lateral.land.inwater =
         (vertical.runoff * network.land.area * 0.001) / lateral.land.dt
 end
@@ -1746,8 +1741,8 @@ Set inflow from the subsurface and land components to a water body (reservoir or
 function set_inflow_waterbody(
     model::Model{N,L,V,R,W,T},
 ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,SurfaceFlow}},V,R,W,T}
-    @unpack lateral, network = model
-    @unpack subsurface, land, river = lateral
+    (; lateral, network) = model
+    (; subsurface, land, river) = lateral
     inds = network.index_river
 
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
@@ -1772,8 +1767,8 @@ Set inflow from the subsurface and land components to a water body (reservoir or
 function set_inflow_waterbody(
     model::Model{N,L,V,R,W,T},
 ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,ShallowWaterRiver}},V,R,W,T}
-    @unpack lateral, network = model
-    @unpack subsurface, land, river = lateral
+    (; lateral, network) = model
+    (; subsurface, land, river) = lateral
     inds = network.index_river
 
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
@@ -1801,8 +1796,8 @@ Set inflow from the subsurface and land components to a water body (reservoir or
 function set_inflow_waterbody(
     model::Model{N,L,V,R,W,T},
 ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T}
-    @unpack lateral, network = model
-    @unpack subsurface, land, river = lateral
+    (; lateral, network) = model
+    (; subsurface, land) = lateral
     inds = network.index_river
 
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
@@ -1820,7 +1815,7 @@ Run surface routing (land and river). Kinematic wave for overland flow and kinem
 local inertial model for river flow.
 """
 function surface_routing(model; ssf_toriver = 0.0)
-    @unpack lateral, network, clock = model
+    (; lateral, network, clock) = model
 
     # run kinematic wave for overland flow
     set_land_inwater(model)
@@ -1846,7 +1841,7 @@ function surface_routing(
     ssf_toriver = 0.0,
 ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T}
 
-    @unpack lateral, vertical, network, clock = model
+    (; lateral, vertical, network, clock) = model
 
     @. lateral.land.runoff = (
         (vertical.net_runoff / 1000.0) * network.land.area / vertical.dt +
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index 988dc4249..7c63bc252 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -385,9 +385,7 @@ function initialize_hbv_model(config::Config)
 end
 
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
-    @unpack lateral, vertical, network, clock, config = model
-
-    inds_riv = network.index_river
+    (; lateral, vertical, network, config) = model
 
     # vertical hbv concept is updated until snow state, after that (optional)
     # snow transport is possible
@@ -412,7 +410,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
 end
 
 function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
-    @unpack lateral, config = model
+    (; lateral, config) = model
 
     reinit = get(config.model, "reinit", true)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -422,7 +420,7 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
         @info "Set initial conditions from state file `$instate_path`."
         set_states(instate_path, model; type = Float)
         # update kinematic wave volume for river and land domain
-        @unpack lateral = model
+        (; lateral) = model
         # makes sure land cells with zero flow width are set to zero q and h
         for i in eachindex(lateral.land.width)
             if lateral.land.width[i] <= 0.0
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index 2d10333ce..5977a667c 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -274,7 +274,7 @@ such as that returned by `Graphs.topological_sort_by_dfs`.
 Returns the material state after transport.
 """
 function accucapacitystate!(material, network, capacity)
-    @unpack graph, order = network
+    (; graph, order) = network
     for v in order
         downstream_nodes = outneighbors(graph, v)
         n = length(downstream_nodes)
@@ -313,7 +313,7 @@ interface, and the order is a valid topological ordering such as that returned b
 Returns the flux (material leaving each cell), and material (left after transport).
 """
 function accucapacityflux!(flux, material, network, capacity)
-    @unpack graph, order = network
+    (; graph, order) = network
     for v in order
         downstream_nodes = outneighbors(graph, v)
         n = length(downstream_nodes)
diff --git a/src/io.jl b/src/io.jl
index 16d89baed..3b7c63499 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -195,8 +195,8 @@ end
 mover_params = (symbols"vertical.precipitation", symbols"vertical.potential_evaporation")
 
 function load_fixed_forcing(model)
-    @unpack reader, network, config = model
-    @unpack forcing_parameters = reader
+    (; reader, network, config) = model
+    (; forcing_parameters) = reader
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -241,8 +241,8 @@ end
 
 "Get dynamic netCDF input for the given time"
 function update_forcing!(model)
-    @unpack vertical, clock, reader, network, config = model
-    @unpack dataset, dataset_times, forcing_parameters = reader
+    (; clock, reader, network, config) = model
+    (; dataset, dataset_times, forcing_parameters) = reader
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -331,8 +331,8 @@ end
 
 "Get cyclic netCDF input for the given time"
 function update_cyclic!(model)
-    @unpack vertical, clock, reader, network, config = model
-    @unpack cyclic_dataset, cyclic_times, cyclic_parameters = reader
+    (; clock, reader, network) = model
+    (; cyclic_dataset, cyclic_times, cyclic_parameters) = reader
 
     # pick up the data that is valid for the past model time step
     month_day = monthday(clock.time - clock.dt)
@@ -1058,7 +1058,7 @@ end
 
 "Write a new timestep with scalar data to a netCDF file"
 function write_netcdf_timestep(model, dataset)
-    @unpack writer, clock, config = model
+    (; writer, clock, config) = model
 
     time_index = add_time(dataset, clock.time)
     for (nt, nc) in zip(writer.nc_scalar, config.netcdf.variable)
@@ -1090,13 +1090,13 @@ end
 
 "Write a new timestep with grid data to a netCDF file"
 function write_netcdf_timestep(model, dataset, parameters)
-    @unpack vertical, clock, reader, network = model
+    (; clock, network) = model
 
     time_index = add_time(dataset, clock.time)
 
     buffer = zeros(Union{Float,Missing}, size(model.network.land.reverse_indices))
     for (key, val) in parameters
-        @unpack par, vector = val
+        (; par, vector) = val
         sel = active_indices(network, par)
         # write the active cells vector to the 2d buffer matrix
         elemtype = eltype(vector)
@@ -1128,8 +1128,8 @@ write_netcdf_timestep(model, dataset::Nothing) = model
 
 "Write model output"
 function write_output(model)
-    @unpack vertical, clock, reader, network, writer = model
-    @unpack dataset, dataset_scalar, parameters = writer
+    (; writer) = model
+    (; dataset, dataset_scalar, parameters) = writer
 
     write_csv_row(model)
     write_netcdf_timestep(model, dataset, parameters)
@@ -1183,7 +1183,7 @@ end
 
 "Close input and output datasets that are opened on model initialization"
 function close_files(model; delete_output::Bool = false)
-    @unpack reader, writer, config = model
+    (; reader, writer, config) = model
 
     close(reader.dataset)
     if haskey(config.input, "cyclic")
@@ -1307,7 +1307,7 @@ function reducer(col, rev_inds, x_nc, y_nc, config, dataset, fileformat)
 end
 
 function write_csv_row(model)
-    @unpack writer, clock, config = model
+    (; writer, clock, config) = model
     isnothing(writer.csv_path) && return nothing
     io = writer.csv_io
     print(io, string(clock.time))
@@ -1630,7 +1630,7 @@ end
 
 "Get `index` for dimension name `layer` or `classes` based on `model`"
 function get_index_dimension(var, model)::Int
-    @unpack vertical = model
+    (; vertical) = model
     if haskey(var, "layer")
         inds = collect(1:vertical.maxlayers)
         index = inds[var["layer"]]
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 2dfe4a2ae..2e4913c83 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -526,7 +526,7 @@ end
 
 "update the sbm_gwf model for a single timestep"
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network, config) = model
 
     do_water_demand = haskey(config.model, "water_demand")
     inds_riv = network.index_river
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index ec1032894..0ad9054f0 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -434,7 +434,7 @@ end
 "update SBM model for a single timestep"
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
 
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network, config) = model
     do_water_demand = haskey(config.model, "water_demand")
     ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
 
@@ -465,7 +465,7 @@ Update SBM model until recharge for a single timestep. This function is also acc
 through BMI, to couple the SBM model to an external groundwater model.
 """
 function update_until_recharge(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network, config) = model
 
     do_water_demand = haskey(config.model, "water_demand")
 
@@ -511,7 +511,7 @@ accessible through BMI, to couple the SBM model to an external groundwater model
 function update_after_subsurfaceflow(
     model::Model{N,L,V,R,W,T},
 ) where {N,L,V,R,W,T<:SbmModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical) = model
 
     # update vertical sbm concept (runoff, ustorelayerdepth and satwaterdepth)
     update_after_subsurfaceflow(
@@ -532,7 +532,7 @@ Update of the total water storage at the end of each timestep per model cell.
 This is done here at model level.
 """
 function update_total_water_storage(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network) = model
 
     # Update the total water storage based on vertical states
     # TODO Maybe look at routing in the near future
@@ -549,7 +549,7 @@ end
 function set_states(
     model::Model{N,L,V,R,W,T},
 ) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
-    @unpack lateral, vertical, network, config = model
+    (; lateral, vertical, network, config) = model
 
     reinit = get(config.model, "reinit", true)::Bool
     routing_options = ("kinematic-wave", "local-inertial")
diff --git a/src/sediment.jl b/src/sediment.jl
index dd708a32a..13a919114 100644
--- a/src/sediment.jl
+++ b/src/sediment.jl
@@ -1133,7 +1133,7 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
 end
 
 function update(rs::RiverSediment, network, config)
-    @unpack graph, order = network
+    (; graph, order) = network
     tcmethod = get(config.model, "rivtransportmethod", "bagnold")::String
 
     # River sediment loads are separated into different particle class.
diff --git a/src/sediment_model.jl b/src/sediment_model.jl
index e2229a1bd..2c5d5ebd6 100644
--- a/src/sediment_model.jl
+++ b/src/sediment_model.jl
@@ -155,7 +155,7 @@ function initialize_sediment_model(config::Config)
 end
 
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SedimentModel}
-    @unpack lateral, vertical, network, clock, config = model
+    (; lateral, vertical, network, config) = model
 
     update_until_ols(vertical, config)
     update_until_oltransport(vertical, config)
@@ -194,7 +194,7 @@ end
 
 function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SedimentModel}
     # read and set states in model object if reinit=false
-    @unpack config = model
+    (; config) = model
     reinit = get(config.model, "reinit", true)::Bool
     if reinit == false
         instate_path = input_path(config, config.state.path_input)
diff --git a/src/utils.jl b/src/utils.jl
index 1f557b1cb..b81248c9e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -137,7 +137,7 @@ and set states in `model` object. Active cells are selected with the correspondi
 - `type = nothing`: type to convert data to after reading. By default no conversion is done.
 """
 function set_states(instate_path, model; type = nothing, dimname = nothing)
-    @unpack network, config = model
+    (; network, config) = model
 
     # Check if required states are covered
     state_ncnames = check_states(config)
diff --git a/src/water_demand.jl b/src/water_demand.jl
index 32e364e0c..432db1b69 100644
--- a/src/water_demand.jl
+++ b/src/water_demand.jl
@@ -536,7 +536,7 @@ function update_water_allocation(
     model::Model{N,L,V,R,W,T},
 ) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
 
-    @unpack network, lateral, vertical = model
+    (; network, lateral, vertical) = model
 
     river = lateral.river
     index_river = network.land.index_river_wb
diff --git a/test/io.jl b/test/io.jl
index f692efed3..25e142de3 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -3,7 +3,6 @@ using Dates
 using TOML
 using CFTime
 using Random
-using UnPack
 using LoggingExtras
 
 tomlpath = joinpath(@__DIR__, "sbm_config.toml")
@@ -196,7 +195,7 @@ model = Wflow.initialize_sbm_model(config)
 Wflow.advance!(model.clock)
 Wflow.load_dynamic_input!(model)
 
-@unpack vertical, clock, reader, writer = model
+(; vertical, clock, reader, writer) = model
 
 @testset "output and state names" begin
     ncdims = ("lon", "lat", "layer", "time")
@@ -242,8 +241,8 @@ end
 end
 
 @testset "network" begin
-    @unpack network = model
-    @unpack indices, reverse_indices = model.network.land
+    (; network) = model
+    (; indices, reverse_indices) = model.network.land
     # test if the reverse index reverses the index
     linear_index = 100
     cartesian_index = indices[linear_index]
@@ -252,7 +251,7 @@ end
 end
 
 @testset "initial parameter values" begin
-    @unpack vertical = model
+    (; vertical) = model
     @test vertical.cfmax[1] ≈ 3.7565300464630127
     @test vertical.soilthickness[1] ≈ 2000.0
     @test vertical.precipitation[49951] ≈ 2.2100000381469727
@@ -280,7 +279,7 @@ Wflow.advance!(model.clock)
 Wflow.load_dynamic_input!(model)
 
 @testset "changed parameter values" begin
-    @unpack vertical = model
+    (; vertical) = model
     @test vertical.cfmax[1] == 2.0
     @test vertical.soilthickness[1] ≈ 2000.0 * 3.0 + 100.0
     @test vertical.precipitation[49951] ≈ 1.5 * 2.2100000381469727
diff --git a/test/run_flextopo.jl b/test/run_flextopo.jl
index cfb4bbffb..974fc085c 100644
--- a/test/run_flextopo.jl
+++ b/test/run_flextopo.jl
@@ -3,7 +3,7 @@ tomlpath = joinpath(@__DIR__, "flextopo_config.toml")
 config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_flextopo_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
@@ -82,7 +82,7 @@ config["model"]["select_fast"] = ["fast_no_storage", "fast_storage", "fast_stora
 config["model"]["select_slow"] = ["common_slow_storage"]
 
 model = Wflow.initialize_flextopo_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
index 12d6eb169..b0706f316 100644
--- a/test/run_hbv.jl
+++ b/test/run_hbv.jl
@@ -3,7 +3,7 @@ tomlpath = joinpath(@__DIR__, "hbv_config.toml")
 config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_hbv_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index d50ccb752..fe81c53ec 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -4,7 +4,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_sbm_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
@@ -446,7 +446,7 @@ Wflow.close_files(model, delete_output = false)
 
     @testset "exponential profile" begin
         model = Wflow.initialize_sbm_model(config)
-        @unpack vertical = model
+        (; vertical) = model
         z = vertical.zi[i]
         kv_z = Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential")
         @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv_0[i] * exp(-vertical.f[i] * z)
@@ -465,7 +465,7 @@ Wflow.close_files(model, delete_output = false)
     @testset "exponential constant profile" begin
         config.input.vertical.ksat_profile = "exponential_constant"
         model = Wflow.initialize_sbm_model(config)
-        @unpack vertical = model
+        (; vertical) = model
         z = vertical.zi[i]
         kv_z =
             Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential_constant")
@@ -500,7 +500,7 @@ Wflow.close_files(model, delete_output = false)
     @testset "layered profile" begin
         config.input.vertical.ksat_profile = "layered"
         model = Wflow.initialize_sbm_model(config)
-        @unpack vertical = model
+        (; vertical) = model
         z = vertical.zi[i]
         @test Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "layered") ≈
               vertical.kv[100][2]
@@ -513,7 +513,7 @@ Wflow.close_files(model, delete_output = false)
     @testset "layered exponential profile" begin
         config.input.vertical.ksat_profile = "layered_exponential"
         model = Wflow.initialize_sbm_model(config)
-        @unpack vertical = model
+        (; vertical) = model
         z = vertical.zi[i]
         @test Wflow.hydraulic_conductivity_at_depth(
             vertical,
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index b984fed2b..3af207643 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -3,7 +3,7 @@ tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
 config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_sbm_gwf_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
@@ -144,7 +144,7 @@ config = Wflow.Config(tomlpath)
 config.model.reinit = false
 
 model = Wflow.initialize_sbm_gwf_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 model = Wflow.run_timestep(model)
diff --git a/test/run_sediment.jl b/test/run_sediment.jl
index 7d009bddc..f8b323996 100644
--- a/test/run_sediment.jl
+++ b/test/run_sediment.jl
@@ -4,7 +4,7 @@ tomlpath = joinpath(@__DIR__, "sediment_config.toml")
 config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_sediment_model(config)
-@unpack network = model
+(; network) = model
 
 model = Wflow.run_timestep(model)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 357608a23..67947e377 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,6 @@ using NCDatasets
 using StaticArrays
 using Statistics
 using Test
-using UnPack
 using Wflow
 using Base.MathConstants: eulergamma
 using Base.Threads
diff --git a/test/subdomains.jl b/test/subdomains.jl
index bb5517640..12b5170b3 100644
--- a/test/subdomains.jl
+++ b/test/subdomains.jl
@@ -3,7 +3,7 @@ config = Wflow.Config(tomlpath)
 
 model = Wflow.initialize_sbm_model(config)
 
-@unpack network = model
+(; network) = model
 
 min_sto_river = get(config.model, "min_streamorder_river", 6)
 min_sto_land = get(config.model, "min_streamorder_land", 5)