diff --git a/Manifest.toml b/Manifest.toml
index d8ff48e29..8523d20bf 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,13 +2,38 @@
 
 julia_version = "1.10.7"
 manifest_format = "2.0"
-project_hash = "986b47155507036c06218d23750fe8de28eb001a"
+project_hash = "0b6ba2f31af23dc8f1958221fc8ee6c1fb4f0096"
 
 [[deps.AbstractTrees]]
 git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 version = "0.4.5"
 
+[[deps.Accessors]]
+deps = ["CompositionsBase", "ConstructionBase", "InverseFunctions", "LinearAlgebra", "MacroTools", "Markdown"]
+git-tree-sha1 = "96bed9b1b57cf750cca50c311a197e306816a1cc"
+uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+version = "0.1.39"
+
+    [deps.Accessors.extensions]
+    AccessorsAxisKeysExt = "AxisKeys"
+    AccessorsDatesExt = "Dates"
+    AccessorsIntervalSetsExt = "IntervalSets"
+    AccessorsStaticArraysExt = "StaticArrays"
+    AccessorsStructArraysExt = "StructArrays"
+    AccessorsTestExt = "Test"
+    AccessorsUnitfulExt = "Unitful"
+
+    [deps.Accessors.weakdeps]
+    AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
+    Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+    StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+    Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
 git-tree-sha1 = "50c3c56a52972d78e8be9fd135bfb91c9574c140"
@@ -134,6 +159,30 @@ deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 version = "1.1.1+0"
 
+[[deps.CompositionsBase]]
+git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
+uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
+version = "0.1.2"
+weakdeps = ["InverseFunctions"]
+
+    [deps.CompositionsBase.extensions]
+    CompositionsBaseInverseFunctionsExt = "InverseFunctions"
+
+[[deps.ConstructionBase]]
+git-tree-sha1 = "76219f1ed5771adbb096743bff43fb5fdd4c1157"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.5.8"
+
+    [deps.ConstructionBase.extensions]
+    ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseLinearAlgebraExt = "LinearAlgebra"
+    ConstructionBaseStaticArraysExt = "StaticArrays"
+
+    [deps.ConstructionBase.weakdeps]
+    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [[deps.CpuId]]
 deps = ["Markdown"]
 git-tree-sha1 = "fcbb72b032692610bfbdb15018ac16a36cf2e406"
@@ -228,6 +277,19 @@ version = "0.1.5"
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[deps.InverseFunctions]]
+git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb"
+uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
+version = "0.1.17"
+
+    [deps.InverseFunctions.extensions]
+    InverseFunctionsDatesExt = "Dates"
+    InverseFunctionsTestExt = "Test"
+
+    [deps.InverseFunctions.weakdeps]
+    Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+    Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
 git-tree-sha1 = "be3dc50a92e5a386872a493a10050136d4703f9b"
diff --git a/Project.toml b/Project.toml
index 4bd017e4e..db702593e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Deltares and contributors"]
 version = "0.8.1"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -24,6 +25,7 @@ TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
 
 [compat]
+Accessors = "0.1"
 Aqua = "0.8"
 BasicModelInterface = "0.1"
 CFTime = "0.1"
diff --git a/server/test/sbm_config.toml b/server/test/sbm_config.toml
index 2acfec3fa..45f9e0b7d 100644
--- a/server/test/sbm_config.toml
+++ b/server/test/sbm_config.toml
@@ -19,33 +19,27 @@ path_output = "outstates-moselle.nc"
 # if listed, the variable must be present in the NetCDF or error
 # if not listed, the variable can get a default value if it has one
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
 
-[state.lateral.river.variables]
-h = "h_river"
-h_av = "h_av_river"
-q = "q_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
+reservoir_water__volume = "volume_reservoir"
 
-[state.lateral.subsurface.variables]
-ssf = "ssf"
+subsurface_water__volume_flow_rate  = "ssf"
 
-[state.lateral.land.variables]
-h = "h_land"
-h_av = "h_av_land"
-q = "q_land"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
 
 [input]
 path_forcing = "forcing-moselle.nc"
@@ -57,6 +51,12 @@ ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -95,8 +95,8 @@ river_bank_water__elevation = "RiverZ"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
-reservoir_areas__number = "wflow_reservoirareas"
-reservoir_locations__number = "wflow_reservoirlocs"
+reservoir_area__number = "wflow_reservoirareas"
+reservoir_location__number = "wflow_reservoirlocs"
 reservoir_surface__area = "ResSimpleArea"
 "reservoir_water_demand~required~downstream__volume_flow_rate" = "ResDemand"
 reservoir_water_release-below-spillway__max_volume_flow_rate = "ResMaxRelease"
@@ -108,20 +108,16 @@ subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
-[input.soil_surface_water__vertical_saturated_hydraulic_conductivity]
+[input.parameters.soil_surface_water__vertical_saturated_hydraulic_conductivity]
 netcdf.variable.name = "KsatVer"
 scale = 1.0
 offset = 0.0
@@ -140,31 +136,19 @@ min_streamorder_land = 5
 [output]
 path = "output_moselle.nc"
 
-[output.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[output.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[output.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[output.lateral.river.variables]
-h = "h_river"
-q = "q_river"
-
-[output.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
-
-[output.lateral.subsurface.variables]
-ssf = "ssf"
-
-[output.lateral.land.variables]
-h = "h_land"
-q = "q_land"
+[output.variables]
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+river_water__depth = "h_river"
+river_water__volume_flow_rate = "q_river"
+reservoir_water__volume = "volume_reservoir"
+subsurface_water__volume_flow_rate  = "ssf"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+vertical.interception.variables.canopy_storage = "canopystorage"
 
 [netcdf]
 path = "output_scalar_moselle.nc"
@@ -179,33 +163,33 @@ coordinate.x = 6.255
 coordinate.y = 50.012
 name = "temp_coord"
 location = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[netcdf.variable]]
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [csv]
 path = "output_moselle.csv"
 
 [[csv.column]]
 header = "Q"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 reducer = "maximum"
 
 [[csv.column]]
 header = "volume"
 index = 1
-parameter = "lateral.river.boundary_conditions.reservoir.variables.volume"
+parameter = "reservoir_water__volume"
 
 [[csv.column]]
 coordinate.x = 6.255
 coordinate.y = 50.012
 header = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 coordinate.x = 6.255
@@ -218,17 +202,17 @@ layer = 2
 header = "temp_byindex"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 
 [[csv.column]]
 header = "recharge"
 map = "subcatchment"
-parameter = "vertical.soil.variables.recharge"
+parameter = "soil_water_sat-zone_top__recharge_volume_flux"
 reducer = "mean"
 
 [API]
diff --git a/src/Wflow.jl b/src/Wflow.jl
index f07ec4edd..b54db1c60 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -2,6 +2,7 @@ module Wflow
 
 import BasicModelInterface as BMI
 
+using Accessors: @optic
 using Base.Threads: nthreads
 using CFTime: CFTime, monthday, dayofyear
 using Dates:
diff --git a/src/bmi.jl b/src/bmi.jl
index c97493f92..3022922db 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -2,15 +2,15 @@
 # https://github.com/Deltares/BasicModelInterface.jl
 
 # Mapping of grid identifier to a key, to get the active indices of the model domain.
-# See also function active_indices(network, key::Tuple).
-const grids = Dict{Int, Tuple{Symbol}}(
-    0 => (:reservoir,),
-    1 => (:lake,),
-    2 => (:drain,),
-    3 => (:river,),
-    4 => (:land,),
-    5 => (:land,),
-    6 => (:land,),
+# See also function active_indices(network, key::AbstractString).
+const grids = Dict{Int, String}(
+    0 => "reservoir",
+    1 => "lake",
+    2 => "drain",
+    3 => "river",
+    4 => "land",
+    5 => "land",
+    6 => "land",
 )
 
 """
@@ -247,7 +247,7 @@ function BMI.get_value_ptr(model::Model, name::String)
     s = split(name, "[")
     key = symbols(first(s))
     if exchange(param(model, key))
-        n = length(active_indices(network, key))
+        n = length(active_indices(network, first(s)))
         if occursin("[", name)
             ind = tryparse(Int, split(s[end], "]")[1])
             if eltype(param(model, key)) <: SVector
diff --git a/src/demand/water_demand.jl b/src/demand/water_demand.jl
index 156d9f895..e4bec8f24 100644
--- a/src/demand/water_demand.jl
+++ b/src/demand/water_demand.jl
@@ -32,24 +32,14 @@ get_demand_gross(model::NoNonIrrigationDemand) = 0.0
 
 "Initialize non-irrigation water demand model for a water use `sector`"
 function NonIrrigationDemand(dataset, config, indices, dt, sector)
+    lens = lens_input_parameter("land~$(sector)__gross_water_demand_flux")
     demand_gross =
-        ncread(
-            dataset,
-            config,
-            "vertical.demand.$(sector).demand_gross";
-            sel = indices,
-            defaults = 0.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float) .*
+        (dt / basetimestep)
+    lens = lens_input_parameter("land~$(sector)__net_water_demand_flux")
     demand_net =
-        ncread(
-            dataset,
-            config,
-            "vertical.demand.$(sector).demand_net";
-            sel = indices,
-            defaults = 0.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float) .*
+        (dt / basetimestep)
     n = length(indices)
     returnflow_f = return_flow_fraction.(demand_gross, demand_net)
 
@@ -84,41 +74,33 @@ end
 
 "Initialize non-paddy irrigation model"
 function NonPaddy(dataset, config, indices, dt)
-    efficiency = ncread(
-        dataset,
-        config,
-        "vertical.demand.nonpaddy.parameters.irrigation_efficiency";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
-    )
+    lens = lens_input_parameter("land~irrigated-non-paddy__irrigation_efficiency")
+    efficiency = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+
+    lens = lens_input_parameter("land~irrigated-non-paddy_area__number")
     areas = ncread(
         dataset,
         config,
-        "vertical.demand.nonpaddy.parameters.irrigation_areas";
+        lens;
         sel = indices,
         defaults = 1,
         optional = false,
         type = Int,
     )
+    lens = lens_input_parameter("land~irrigated-non-paddy__irrigation_trigger_flag")
     irrigation_trigger = ncread(
         dataset,
         config,
-        "vertical.demand.nonpaddy.parameters.irrigation_trigger";
+        lens;
         sel = indices,
         defaults = 1,
         optional = false,
         type = Bool,
     )
+    lens = lens_input_parameter("land~irrigated-non-paddy__max_irrigation_rate")
     max_irri_rate =
-        ncread(
-            dataset,
-            config,
-            "vertical.demand.nonpaddy.parameters.maximum_irrigation_rate";
-            sel = indices,
-            defaults = 25.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 25.0, type = Float) .*
+        (dt / basetimestep)
 
     params = NonPaddyParameters{Float}(;
         maximum_irrigation_rate = max_irri_rate,
@@ -230,63 +212,28 @@ end
 
 "Initialize paddy irrigation model"
 function Paddy(dataset, config, indices, dt)
-    h_min = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.h_min";
-        sel = indices,
-        defaults = 20.0,
-        type = Float,
-    )
-    h_opt = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.h_opt";
-        sel = indices,
-        defaults = 50.0,
-        type = Float,
-    )
-    h_max = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.h_max";
-        sel = indices,
-        defaults = 80.0,
-        type = Float,
-    )
-    efficiency = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.irrigation_efficiency";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
-    )
-    areas = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.irrigation_areas";
-        sel = indices,
-        optional = false,
-        type = Bool,
-    )
-    irrigation_trigger = ncread(
-        dataset,
-        config,
-        "vertical.demand.paddy.parameters.irrigation_trigger";
-        sel = indices,
-        optional = false,
-        type = Bool,
-    )
+    lens = lens_input_parameter("land~irrigated-paddy__min_depth")
+    h_min = ncread(dataset, config, lens; sel = indices, defaults = 20.0, type = Float)
+
+    lens = lens_input_parameter("land~irrigated-paddy__optimal_depth")
+    h_opt = ncread(dataset, config, lens; sel = indices, defaults = 50.0, type = Float)
+
+    lens = lens_input_parameter("land~irrigated-paddy__max_depth")
+    h_max = ncread(dataset, config, lens; sel = indices, defaults = 80.0, type = Float)
+
+    lens = lens_input_parameter("land~irrigated-paddy__irrigation_efficiency")
+    efficiency = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+
+    lens = lens_input_parameter("land~irrigated-paddy_area__number")
+    areas = ncread(dataset, config, lens; sel = indices, optional = false, type = Bool)
+
+    lens = lens_input_parameter("land~irrigated-paddy__irrigation_trigger_flag")
+    irrigation_trigger =
+        ncread(dataset, config, lens; sel = indices, optional = false, type = Bool)
+    lens = lens_input_parameter("land~irrigate-paddy__max_irrigation_rate")
     max_irri_rate =
-        ncread(
-            dataset,
-            config,
-            "vertical.demand.paddy.parameters.maximum_irrigation_rate";
-            sel = indices,
-            defaults = 25.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 25.0, type = Float) .*
+        (dt / basetimestep)
     n = length(indices)
     params = PaddyParameters{Float}(;
         irrigation_efficiency = efficiency,
@@ -545,22 +492,11 @@ end
 
 "Initialize water allocation for the land domain"
 function AllocationLand(dataset, config, indices)
-    frac_sw_used = ncread(
-        dataset,
-        config,
-        "vertical.allocation.parameters.frac_sw_used";
-        sel = indices,
-        defaults = 1,
-        type = Float,
-    )
-    areas = ncread(
-        dataset,
-        config,
-        "vertical.allocation.parameters.areas";
-        sel = indices,
-        defaults = 1,
-        type = Int,
-    )
+    lens = lens_input_parameter("land_surface_water__withdrawal_fraction")
+    frac_sw_used = ncread(dataset, config, lens; sel = indices, defaults = 1, type = Float)
+
+    lens = lens_input_parameter("land_water_allocation_area__number")
+    areas = ncread(dataset, config, lens; sel = indices, defaults = 1, type = Int)
 
     n = length(indices)
 
diff --git a/src/glacier/glacier.jl b/src/glacier/glacier.jl
index 118fd6055..dbc0ea4c0 100644
--- a/src/glacier/glacier.jl
+++ b/src/glacier/glacier.jl
@@ -10,10 +10,11 @@ end
 
 "Initialize glacier model variables"
 function GlacierVariables(dataset, config, indices)
+    lens = lens_input_parameter("glacier_ice__leq-volume")
     glacier_store = ncread(
         dataset,
         config,
-        "vertical.glacier.variables.glacier_store";
+        lens;
         sel = indices,
         defaults = 5500.0,
         type = Float,
@@ -55,39 +56,44 @@ struct NoGlacierModel{T} <: AbstractGlacierModel{T} end
 
 "Initialize glacier HBV model parameters"
 function GlacierHbvParameters(dataset, config, indices, dt)
+    lens = lens_input_parameter("glacier_ice__melting_temperature_threshold")
     g_ttm = ncread(
         dataset,
         config,
-        "glacier_ice__melting_temperature_threshold";
+        lens;
         sel = indices,
         defaults = 0.0,
         type = Float,
         fill = 0.0,
     )
+    lens = lens_input_parameter("glacier_ice__degree-day_coefficient")
     g_cfmax =
         ncread(
             dataset,
             config,
-            "glacier_ice__degree-day_coefficient";
+            lens;
             sel = indices,
             defaults = 3.0,
             type = Float,
             fill = 0.0,
         ) .* (dt / basetimestep)
+    lens =
+        lens_input_parameter("glacier_firn_accumulation__snowpack~dry_leq-depth_fraction")
     g_sifrac =
         ncread(
             dataset,
             config,
-            "glacier_firn_accumulation__snowpack~dry_leq-depth_fraction";
+            lens;
             sel = indices,
             defaults = 0.001,
             type = Float,
             fill = 0.0,
         ) .* (dt / basetimestep)
+    lens = lens_input_parameter("glacier_surface__area_fraction")
     glacier_frac = ncread(
         dataset,
         config,
-        "glacier_surface__area_fraction";
+        lens;
         sel = indices,
         defaults = 0.0,
         type = Float,
diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl
index 7b6689f4c..ab76cd69f 100644
--- a/src/groundwater/aquifer.jl
+++ b/src/groundwater/aquifer.jl
@@ -127,28 +127,18 @@ instead. Specific yield will vary roughly between 0.05 (clay) and 0.45 (peat)
 end
 
 function UnconfinedAquiferParameters(dataset, config, indices, top, bottom, area)
-    k = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.conductivity";
-        sel = indices,
-        type = Float,
+    lens = lens_input_parameter(
+        "subsurface_surface_water__horizontal_saturated_hydraulic_conductivity",
     )
-    specific_yield = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.specific_yield";
-        sel = indices,
-        type = Float,
-    )
-    f = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.gwf_f";
-        sel = indices,
-        type = Float,
-        defaults = 3.0,
+    k = ncread(dataset, config, lens; sel = indices, type = Float)
+
+    lens = lens_input_parameter("subsurface_water__specific_yield")
+    specific_yield = ncread(dataset, config, lens; sel = indices, type = Float)
+
+    lens = lens_input_parameter(
+        "subsurface__horizontal_saturated_hydraulic_conductivity_scale_parameter",
     )
+    f = ncread(dataset, config, lens; sel = indices, type = Float, defaults = 3.0)
 
     parameters =
         UnconfinedAquiferParameters{Float}(; k, top, bottom, area, specific_yield, f)
@@ -380,14 +370,9 @@ end
 end
 
 function ConstantHead(dataset, config, indices)
-    constanthead = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.constant_head";
-        sel = indices,
-        type = Float,
-        fill = mv,
-    )
+    lens = lens_input_parameter("model_boundary_condition~constant_hydraulic_head")
+    constanthead = ncread(dataset, config, lens; sel = indices, type = Float, fill = mv)
+
     n = length(indices)
     index_constanthead = filter(i -> !isequal(constanthead[i], mv), 1:n)
     head = constanthead[index_constanthead]
diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl
index 52927733c..9e73d42a9 100644
--- a/src/groundwater/boundary_conditions.jl
+++ b/src/groundwater/boundary_conditions.jl
@@ -34,27 +34,14 @@ end
 end
 
 function River(dataset, config, indices, index)
-    infiltration_conductance = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.infiltration_conductance";
-        sel = indices,
-        type = Float,
-    )
-    exfiltration_conductance = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.exfiltration_conductance";
-        sel = indices,
-        type = Float,
-    )
-    bottom = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.river_bottom";
-        sel = indices,
-        type = Float,
-    )
+    lens = lens_input_parameter("river_water__infiltration_conductance")
+    infiltration_conductance = ncread(dataset, config, lens; sel = indices, type = Float)
+
+    lens = lens_input_parameter("river_water__exfiltration_conductance")
+    exfiltration_conductance = ncread(dataset, config, lens; sel = indices, type = Float)
+
+    lens = lens_input_parameter("river_bottom__elevation")
+    bottom = ncread(dataset, config, lens; sel = indices, type = Float)
 
     parameters =
         RiverParameters{Float}(infiltration_conductance, exfiltration_conductance, bottom)
@@ -97,22 +84,12 @@ end
 end
 
 function Drainage(dataset, config, indices, index)
-    drain_elevation = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.drain_elevation";
-        sel = indices,
-        type = Float,
-        fill = mv,
-    )
-    drain_conductance = ncread(
-        dataset,
-        config,
-        "lateral.subsurface.drain_conductance";
-        sel = indices,
-        type = Float,
-        fill = mv,
-    )
+    lens = lens_input_parameter("land_drain__elevation")
+    drain_elevation = ncread(dataset, config, lens; sel = indices, type = Float, fill = mv)
+
+    lens = lens_input_parameter("land_drain__conductance")
+    drain_conductance =
+        ncread(dataset, config, lens; sel = indices, type = Float, fill = mv)
     elevation = drain_elevation[index]
     conductance = drain_conductance[index]
     parameters = DrainageParameters{Float}(; elevation, conductance)
diff --git a/src/io.jl b/src/io.jl
index 1b3f041da..7c58540c3 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -26,7 +26,6 @@ function param(obj, fields, default)
         return default
     end
 end
-get_fields(name) = haskey(standard_name_map, name) ? standard_name_map[name] : name
 
 """
     Config(path::AbstractString)
@@ -60,10 +59,20 @@ function Base.setproperty!(config::Config, f::Symbol, x)
     return nothing
 end
 
+"Access nested Config object using a `lens`"
+function _lens(obj::Config, lens::ComposedFunction, default)
+    try
+        l = lens(obj)
+        return l isa AbstractDict ? Config(l, pathof(obj)) : l
+    catch
+        return default
+    end
+end
+
 "Get a value from the Config with either the key or an alias of the key."
-function get_alias(config::Config, key, alias, default)
-    alias_or_default = param(config, alias, default)
-    return param(config, key, alias_or_default)
+function get_alias(config::Config, key::ComposedFunction, alias::ComposedFunction, default)
+    alias_or_default = _lens(config, alias, default)
+    return _lens(config, key, alias_or_default)
 end
 
 "Get a value from the Config with a key, throwing an error if it is not one of the options."
@@ -184,25 +193,25 @@ end
 
 function get_param_res(model)
     return Dict(
-        symbols"vertical.atmospheric_forcing.precipitation" =>
+        "atmosphere_water__precipitation_volume_flux" =>
             model.lateral.river.boundary_conditions.reservoir.boundary_conditions.precipitation,
-        symbols"vertical.atmospheric_forcing.potential_evaporation" =>
+        "land_surface_water__potential_evaporation_volume_flux" =>
             model.lateral.river.boundary_conditions.reservoir.boundary_conditions.evaporation,
     )
 end
 
 function get_param_lake(model)
     return Dict(
-        symbols"vertical.atmospheric_forcing.precipitation" =>
+        "atmosphere_water__precipitation_volume_flux" =>
             model.lateral.river.boundary_conditions.lake.boundary_conditions.precipitation,
-        symbols"vertical.atmospheric_forcing.potential_evaporation" =>
+        "land_surface_water__potential_evaporation_volume_flux" =>
             model.lateral.river.boundary_conditions.lake.boundary_conditions.evaporation,
     )
 end
 
 mover_params = (
-    symbols"vertical.atmospheric_forcing.precipitation",
-    symbols"vertical.atmospheric_forcing.potential_evaporation",
+    "atmosphere_water__precipitation_volume_flux",
+    "land_surface_water__potential_evaporation_volume_flux",
 )
 
 function load_fixed_forcing!(model)
@@ -225,7 +234,8 @@ function load_fixed_forcing!(model)
     for (par, ncvar) in forcing_parameters
         if ncvar.name === nothing
             val = ncvar.value * ncvar.scale + ncvar.offset
-            param_vector = param(model, par)
+            lens = standard_name_map[par]
+            param_vector = lens(model)
             param_vector .= val
             # set fixed precipitation and evaporation over the lakes and reservoirs and put
             # these into the lakes and reservoirs structs and set the precipitation and
@@ -301,8 +311,8 @@ function update_forcing!(model)
                 end
             end
         end
-
-        param_vector = param(model, par)
+        lens = standard_name_map[par]
+        param_vector = lens(model)
         sel = active_indices(network, par)
         data_sel = data[sel]
         if any(ismissing, data_sel)
@@ -334,7 +344,7 @@ monthday_passed(curr, avail) = (curr[1] >= avail[1]) && (curr[2] >= avail[2])
 "Get dynamic and cyclic netCDF input"
 function load_dynamic_input!(model)
     update_forcing!(model)
-    if haskey(model.config.input, "cyclic")
+    if haskey(model.config.input.dynamic, "cyclic")
         update_cyclic!(model)
     end
     return nothing
@@ -358,7 +368,8 @@ function update_cyclic!(model)
                 error("Could not find applicable cyclic timestep for $month_day")
             # load from netCDF into the model according to the mapping
             data = get_at(cyclic_dataset, ncvar.name, i)
-            param_vector = param(model, par)
+            lens = standard_name_map[par]
+            param_vector = lens(model)
             sel = active_indices(network, par)
             param_vector .= data[sel]
             if ncvar.scale != 1.0 || ncvar.offset != 0.0
@@ -428,7 +439,12 @@ function setup_scalar_netcdf(
             (nc.location_dim,);
             attrib = ["cf_role" => "timeseries_id"],
         )
-        v = param(modelmap, nc.par)
+        v = if haskey(standard_name_map, nc.par)
+            lens = standard_name_map[nc.par]
+            lens(modelmap)
+        else
+            param(modelmap, nc.par)
+        end
         if eltype(v) <: AbstractFloat
             defVar(
                 ds,
@@ -629,9 +645,9 @@ struct NCReader{T}
     dataset::CFDataset
     dataset_times::Vector{T}
     cyclic_dataset::Union{NCDataset, Nothing}
-    cyclic_times::Dict{Tuple{Symbol, Vararg{Symbol}}, Vector{Tuple{Int, Int}}}
-    forcing_parameters::Dict{Tuple{Symbol, Vararg{Symbol}}, NamedTuple}
-    cyclic_parameters::Dict{Tuple{Symbol, Vararg{Symbol}}, NamedTuple}
+    cyclic_times::Dict{String, Vector{Tuple{Int, Int}}}
+    forcing_parameters::Dict{String, NamedTuple}
+    cyclic_parameters::Dict{String, NamedTuple}
 end
 
 struct Writer
@@ -697,14 +713,13 @@ function prepare_reader(config)
     end
 
     # check for cyclic parameters
-    do_cyclic = haskey(config.input, "cyclic")
+    do_cyclic = haskey(config.input.dynamic, "cyclic")
 
     # create map from internal location to netCDF variable name for forcing parameters
-    forcing_parameters = Dict{Tuple{Symbol, Vararg{Symbol}}, NamedTuple}()
-    for par in config.input.forcing
-        fields = symbols(par)
-        ncname, mod = ncvar_name_modifier(param(config.input, fields))
-        forcing_parameters[fields] =
+    forcing_parameters = Dict{String, NamedTuple}()
+    for par in config.input.dynamic.forcing
+        ncname, mod = ncvar_name_modifier(param(config.input.forcing, symbols(par)))
+        forcing_parameters[par] =
             (name = ncname, scale = mod.scale, offset = mod.offset, value = mod.value)
 
         @info "Set `$par` using netCDF variable `$ncname` as forcing parameter."
@@ -716,24 +731,23 @@ function prepare_reader(config)
     # (memory usage))
     if do_cyclic == true
         cyclic_dataset = NCDataset(cyclic_path)
-        cyclic_parameters = Dict{Tuple{Symbol, Vararg{Symbol}}, NamedTuple}()
-        cyclic_times = Dict{Tuple{Symbol, Vararg{Symbol}}, Vector{Tuple{Int, Int}}}()
-        for par in config.input.cyclic
-            fields = symbols(get_fields(par)) #TODO: make this more restrict (only allow variables in `standard_name_map`)
-            ncname, mod = ncvar_name_modifier(param(config.input, par))
+        cyclic_parameters = Dict{String, NamedTuple}()
+        cyclic_times = Dict{String, Vector{Tuple{Int, Int}}}()
+        for par in config.input.dynamic.cyclic
+            #fields = standard_name_map[par]
+            ncname, mod = ncvar_name_modifier(param(config.input.parameters, par))
             i = findfirst(x -> startswith(x, "time"), dimnames(cyclic_dataset[ncname]))
             dimname = dimnames(cyclic_dataset[ncname])[i]
             cyclic_nc_times = collect(cyclic_dataset[dimname])
-            cyclic_times[fields] = timecycles(cyclic_nc_times)
-            cyclic_parameters[fields] =
-                (name = ncname, scale = mod.scale, offset = mod.offset)
+            cyclic_times[par] = timecycles(cyclic_nc_times)
+            cyclic_parameters[par] = (name = ncname, scale = mod.scale, offset = mod.offset)
 
             @info "Set `$par` using netCDF variable `$ncname` as cyclic parameter, with `$(length(cyclic_nc_times))` timesteps."
         end
     else
-        cyclic_parameters = Dict{Tuple{Symbol, Vararg{Symbol}}, NamedTuple}()
+        cyclic_parameters = Dict{String, NamedTuple}()
         cyclic_dataset = nothing
-        cyclic_times = Dict{Tuple{Symbol, Vararg{Symbol}}, Vector{Tuple{Int, Int}}}()
+        cyclic_times = Dict{String, Vector{Tuple{Int, Int}}}()
     end
 
     # check if there is overlap
@@ -759,10 +773,11 @@ end
 
 "Get a Vector of all unique location ids from a 2D map"
 function locations_map(ds, mapname, config)
+    lens = lens_input(mapname)
     map_2d = ncread(
         ds,
         config,
-        mapname;
+        lens;
         optional = false,
         type = Union{Int, Missing},
         allow_missing = true,
@@ -843,7 +858,7 @@ function flat!(d, path, el::Dict)
 end
 
 function flat!(d, path, el)
-    k = symbols(path)
+    k = path
     d[k] = el
     return nothing
 end
@@ -877,11 +892,9 @@ Dict(
 ```
 """
 function ncnames(dict)
-    ncnames_dict = Dict{Tuple{Symbol, Vararg{Symbol}}, String}()
+    ncnames_dict = Dict{String, String}()
     for (k, v) in dict
-        if v isa Dict  # ignore top level values (e.g. output.path)
-            flat!(ncnames_dict, k, v)
-        end
+        flat!(ncnames_dict, k, v)
     end
     return ncnames_dict
 end
@@ -894,7 +907,12 @@ Create a Dict that maps parameter netCDF names to arrays in the Model.
 function out_map(ncnames_dict, modelmap)
     output_map = Dict{String, Any}()
     for (par, ncname) in ncnames_dict
-        A = param(modelmap, par)
+        A = if haskey(standard_name_map, par)
+            lens = standard_name_map[par]
+            lens(modelmap)
+        else
+            param(modelmap, par)
+        end
         output_map[ncname] = (par = par, vector = A)
     end
     return output_map
@@ -936,7 +954,7 @@ function prepare_writer(
         deflatelevel = get(config.output, "compressionlevel", 0)::Int
         @info "Create an output netCDF file `$nc_path` for grid data, using compression level `$deflatelevel`."
         # create a flat mapping from internal parameter locations to netCDF variable names
-        output_ncnames = ncnames(config.output)
+        output_ncnames = ncnames(config.output.variables)
         # fill the output_map by mapping parameter netCDF names to arrays
         output_map = out_map(output_ncnames, modelmap)
         ds = setup_grid_netcdf(
@@ -959,7 +977,9 @@ function prepare_writer(
     # create a separate state output netCDF that will hold the last timestep of all states
     # but only if config.state.path_output has been set
     if haskey(config, "state") && haskey(config.state, "path_output")
-        state_ncnames = check_states(config)
+        # TODO: revert back to state checking
+        # state_ncnames = check_states(config)
+        state_ncnames = ncnames(config.state.variables)
         state_map = out_map(state_ncnames, modelmap)
         nc_state_path = output_path(config, config.state.path_output)
         @info "Create a state output netCDF file `$nc_state_path`."
@@ -1065,7 +1085,12 @@ function write_netcdf_timestep(model, dataset)
 
     time_index = add_time(dataset, clock.time)
     for (nt, nc) in zip(writer.nc_scalar, config.netcdf.variable)
-        A = param(model, nt.parameter)
+        A = if haskey(standard_name_map, nt.parameter)
+            lens = standard_name_map[nt.parameter]
+            lens(model)
+        else
+            param(model, nt.parameter)
+        end
         elemtype = eltype(A)
         # could be a value, or a vector in case of map
         if elemtype <: AbstractFloat
@@ -1189,7 +1214,7 @@ function close_files(model; delete_output::Bool = false)
     (; reader, writer, config) = model
 
     close(reader.dataset)
-    if haskey(config.input, "cyclic")
+    if haskey(config.input.dynamic, "cyclic")
         close(reader.cyclic_dataset)
     end
     writer.dataset === nothing || close(writer.dataset)
@@ -1242,10 +1267,11 @@ function reducer(col, rev_inds, x_nc, y_nc, config, dataset, fileformat)
         # and makes sense in the case of a gauge map
         reducer_name = get(col, "reducer", "only")
         f = reducerfunction(reducer_name)
+        lens = lens_input(mapname)
         map_2d = ncread(
             dataset,
             config,
-            mapname;
+            lens;
             optional = false,
             type = Union{Int, Missing},
             allow_missing = true,
@@ -1315,7 +1341,12 @@ function write_csv_row(model)
     io = writer.csv_io
     print(io, string(clock.time))
     for (nt, col) in zip(writer.csv_cols, config.csv.column)
-        A = param(model, nt.parameter)
+        A = if haskey(standard_name_map, nt.parameter)
+            lens = standard_name_map[nt.parameter]
+            lens(model)
+        else
+            param(model, nt.parameter)
+        end
         # v could be a value, or a vector in case of map
         if eltype(A) <: SVector
             # indexing is required in case of a SVector and CSV output
diff --git a/src/parameters.jl b/src/parameters.jl
index 63254d89b..1ca2eb7d3 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -21,47 +21,20 @@ end
 "Initialize (shared) vegetation parameters"
 function VegetationParameters(dataset, config, indices)
     n = length(indices)
-    rootingdepth = ncread(
-        dataset,
-        config,
-        "vegetation_root__depth";
-        sel = indices,
-        defaults = 750.0,
-        type = Float,
-    )
-    kc = ncread(
-        dataset,
-        config,
-        "vegetation__crop_factor";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
-    )
-    if haskey(config.input, "vegetation__leaf-area_index")
-        storage_specific_leaf = ncread(
-            dataset,
-            config,
-            "vegetation__specific-leaf_storage";
-            optional = false,
-            sel = indices,
-            type = Float,
-        )
-        storage_wood = ncread(
-            dataset,
-            config,
-            "vegetation_woody-part__storage_capacity";
-            optional = false,
-            sel = indices,
-            type = Float,
-        )
-        kext = ncread(
-            dataset,
-            config,
-            "vegetation_canopy__light-extinction_coefficient";
-            optional = false,
-            sel = indices,
-            type = Float,
-        )
+    lens = lens_input_parameter("vegetation_root__depth")
+    rootingdepth =
+        ncread(dataset, config, lens; sel = indices, defaults = 750.0, type = Float)
+    lens = lens_input_parameter("vegetation__crop_factor")
+    kc = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+    if haskey(config.input.parameters, "vegetation__leaf-area_index")
+        lens = lens_input_parameter("vegetation__specific-leaf_storage")
+        storage_specific_leaf =
+            ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
+        lens = lens_input_parameter("vegetation_woody-part__storage_capacity")
+        storage_wood =
+            ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
+        lens = lens_input_parameter("vegetation_canopy__light-extinction_coefficient")
+        kext = ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
         vegetation_parameter_set = VegetationParameters(;
             leaf_area_index = fill(mv, n),
             storage_wood,
@@ -73,22 +46,11 @@ function VegetationParameters(dataset, config, indices)
             kc,
         )
     else
-        canopygapfraction = ncread(
-            dataset,
-            config,
-            "vegetation_canopy__gap_fraction";
-            sel = indices,
-            defaults = 0.1,
-            type = Float,
-        )
-        cmax = ncread(
-            dataset,
-            config,
-            "vegetation_water__storage_capacity";
-            sel = indices,
-            defaults = 1.0,
-            type = Float,
-        )
+        lens = lens_input_parameter("vegetation_canopy__gap_fraction")
+        canopygapfraction =
+            ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
+        lens = lens_input_parameter("vegetation_water__storage_capacity")
+        cmax = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
         vegetation_parameter_set = VegetationParameters(;
             leaf_area_index = nothing,
             storage_wood = nothing,
diff --git a/src/routing/lake.jl b/src/routing/lake.jl
index 1ebaf9ab7..4f2daa37c 100644
--- a/src/routing/lake.jl
+++ b/src/routing/lake.jl
@@ -17,24 +17,13 @@ function LakeParameters(config, dataset, inds_riv, nriv, pits)
     # read only lake data if lakes true
     # allow lakes only in river cells
     # note that these locations are only the lake outlet pixels
-    lakelocs_2d = ncread(
-        dataset,
-        config,
-        "lake_locations__number";
-        optional = false,
-        type = Int,
-        fill = 0,
-    )
+    lens = lens_input("lake_location__number")
+    lakelocs_2d = ncread(dataset, config, lens; optional = false, type = Int, fill = 0)
     lakelocs = lakelocs_2d[inds_riv]
 
     # this holds the same ids as lakelocs, but covers the entire lake
-    lakecoverage_2d = ncread(
-        dataset,
-        config,
-        "lakes_areas__number";
-        optional = false,
-        allow_missing = true,
-    )
+    lens = lens_input("lakes_area__number")
+    lakecoverage_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     # for each lake, a list of 2D indices, needed for getting the mean precipitation
     inds_lake_cov = Vector{CartesianIndex{2}}[]
 
@@ -61,73 +50,74 @@ function LakeParameters(config, dataset, inds_riv, nriv, pits)
         end
     end
 
+    lens = lens_input_parameter("lake_surface__area")
     lakearea = ncread(
         dataset,
         config,
-        "lake_surface__area";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("lake_water__rating_curve_coefficient")
     lake_b = ncread(
         dataset,
         config,
-        "lake_water__rating_curve_coefficient";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("lake_water__rating_curve_exponent")
     lake_e = ncread(
         dataset,
         config,
-        "lake_water__rating_curve_exponent";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("lake_water_flow_threshold-level__elevation")
     lake_threshold = ncread(
         dataset,
         config,
-        "lake_water_flow_threshold-level__elevation";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Float,
         fill = 0,
     )
-    linked_lakelocs = ncread(
-        dataset,
-        config,
-        "lake~lower_locations__number";
-        sel = inds_lake,
-        defaults = 0,
-        type = Int,
-        fill = 0,
-    )
+    lens = lens_input_parameter("lake~lower_location__number")
+    linked_lakelocs =
+        ncread(dataset, config, lens; sel = inds_lake, defaults = 0, type = Int, fill = 0)
+    lens = lens_input_parameter("lake_water__storage_curve_type_number")
     lake_storfunc = ncread(
         dataset,
         config,
-        "lake_water__storage_curve_type_number";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Int,
         fill = 0,
     )
+    lens = lens_input_parameter("lake_water__rating_curve_type_number")
     lake_outflowfunc = ncread(
         dataset,
         config,
-        "lake_water__rating_curve_type_number";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Int,
         fill = 0,
     )
+    lens = lens_input_parameter("lake_water_level__initial_elevation")
     lake_waterlevel = ncread(
         dataset,
         config,
-        "lake_water_level__initial_elevation";
+        lens;
         optional = false,
         sel = inds_lake,
         type = Float,
diff --git a/src/routing/reservoir.jl b/src/routing/reservoir.jl
index 4e5a49383..86eed4c34 100644
--- a/src/routing/reservoir.jl
+++ b/src/routing/reservoir.jl
@@ -13,10 +13,11 @@ function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits
     # read only reservoir data if reservoirs true
     # allow reservoirs only in river cells
     # note that these locations are only the reservoir outlet pixels
+    lens = lens_input_parameter("reservoir_location__number")
     reslocs = ncread(
         dataset,
         config,
-        "reservoir_locations__number";
+        lens;
         optional = false,
         sel = indices_river,
         type = Int,
@@ -24,13 +25,8 @@ function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits
     )
 
     # this holds the same ids as reslocs, but covers the entire reservoir
-    rescoverage_2d = ncread(
-        dataset,
-        config,
-        "reservoir_areas__number";
-        optional = false,
-        allow_missing = true,
-    )
+    lens = lens_input_parameter("reservoir_area__number")
+    rescoverage_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     # for each reservoir, a list of 2D indices, needed for getting the mean precipitation
     inds_res_cov = Vector{CartesianIndex{2}}[]
 
@@ -56,56 +52,63 @@ function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits
             push!(inds_res_cov, res_cov)
         end
     end
-
+    lens =
+        lens_input_parameter("reservoir_water_demand~required~downstream__volume_flow_rate")
     resdemand = ncread(
         dataset,
         config,
-        "reservoir_water_demand~required~downstream__volume_flow_rate";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
         fill = 0,
     )
+    lens =
+        lens_input_parameter("reservoir_water_release-below-spillway__max_volume_flow_rate")
     resmaxrelease = ncread(
         dataset,
         config,
-        "reservoir_water_release-below-spillway__max_volume_flow_rate";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("reservoir_water__max_volume")
     resmaxvolume = ncread(
         dataset,
         config,
-        "reservoir_water__max_volume";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("reservoir_surface__area")
     resarea = ncread(
         dataset,
         config,
-        "reservoir_surface__area";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("reservoir_water~full-target__volume_fraction")
     res_targetfullfrac = ncread(
         dataset,
         config,
-        "reservoir_water~full-target__volume_fraction";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
         fill = 0,
     )
+    lens = lens_input_parameter("reservoir_water~min-target__volume_fraction")
     res_targetminfrac = ncread(
         dataset,
         config,
-        "reservoir_water~min-target__volume_fraction";
+        lens;
         optional = false,
         sel = inds_res,
         type = Float,
diff --git a/src/routing/subsurface.jl b/src/routing/subsurface.jl
index 4e820eb7f..2489ab77f 100644
--- a/src/routing/subsurface.jl
+++ b/src/routing/subsurface.jl
@@ -44,20 +44,17 @@ function LateralSsfParameters(
     flow_length,
     flow_width,
 )
-    khfrac = ncread(
-        dataset,
-        config,
-        "subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
+    lens = lens_input_parameter(
+        "subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio",
     )
+    khfrac = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
     n_cells = length(khfrac)
 
     (; theta_s, theta_r, soilthickness) = soil
     soilthickness = soilthickness .* 0.001
 
-    kh_profile_type = get(config.model, "saturated_hydraulic_conductivity_profile", "exponential")::String
+    kh_profile_type =
+        get(config.model, "saturated_hydraulic_conductivity_profile", "exponential")::String
     dt = Second(config.timestepsecs) / basetimestep
     if kh_profile_type == "exponential"
         (; kv_0, f) = soil.kv_profile
diff --git a/src/routing/surface_kinwave.jl b/src/routing/surface_kinwave.jl
index de6dd0c08..ab08def8f 100644
--- a/src/routing/surface_kinwave.jl
+++ b/src/routing/surface_kinwave.jl
@@ -88,30 +88,14 @@ end
 
 "Initialize river flow model parameters"
 function RiverFlowParameters(dataset, config, indices, river_length, river_width)
-    mannings_n = ncread(
-        dataset,
-        config,
-        "river_water_flow__manning_n_parameter";
-        sel = indices,
-        defaults = 0.036,
-        type = Float,
-    )
-    bankfull_depth = ncread(
-        dataset,
-        config,
-        "river_bank_water__depth";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
-    )
-    slope = ncread(
-        dataset,
-        config,
-        "river__slope";
-        optional = false,
-        sel = indices,
-        type = Float,
-    )
+    lens = lens_input_parameter("river_water_flow__manning_n_parameter")
+    mannings_n =
+        ncread(dataset, config, lens; sel = indices, defaults = 0.036, type = Float)
+    lens = lens_input_parameter("river_bank_water__depth")
+    bankfull_depth =
+        ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+    lens = lens_input_parameter("river__slope")
+    slope = ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
     clamp!(slope, 0.00001, Inf)
 
     flow_parameter_set = ManningFlowParameters(slope, mannings_n, river_length, river_width)
@@ -217,14 +201,9 @@ end
 
 "Initialize Overland flow model `KinWaveOverlandFlow`"
 function KinWaveOverlandFlow(dataset, config, indices; slope, flow_length, flow_width)
-    mannings_n = ncread(
-        dataset,
-        config,
-        "land_surface_water_flow__manning_n_parameter";
-        sel = indices,
-        defaults = 0.072,
-        type = Float,
-    )
+    lens = lens_input_parameter("land_surface_water_flow__manning_n_parameter")
+    mannings_n =
+        ncread(dataset, config, lens; sel = indices, defaults = 0.072, type = Float)
 
     n = length(indices)
     timestepping =
diff --git a/src/routing/surface_local_inertial.jl b/src/routing/surface_local_inertial.jl
index 2e900fdd6..86b5748f9 100644
--- a/src/routing/surface_local_inertial.jl
+++ b/src/routing/surface_local_inertial.jl
@@ -39,46 +39,24 @@ function LocalInertialRiverFlowParameters(
     floodplain_1d = get(config.model, "floodplain_1d", false)::Bool
 
     @info "Local inertial approach is used for river flow." cfl h_thresh froude_limit floodplain_1d
-
-    riverlength_bc = ncread(
-        dataset,
-        config,
-        "model_boundary_condition~river__length";
-        sel = inds_pit,
-        defaults = 1.0e04,
-        type = Float,
-    )
-    bankfull_elevation_2d = ncread(
-        dataset,
-        config,
-        "river_bank_water__elevation";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
-    bankfull_depth_2d = ncread(
-        dataset,
-        config,
-        "river_bank_water__depth";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+    lens = lens_input_parameter("model_boundary_condition~river__length")
+    riverlength_bc =
+        ncread(dataset, config, lens; sel = inds_pit, defaults = 1.0e04, type = Float)
+    lens = lens_input_parameter("river_bank_water__elevation")
+    bankfull_elevation_2d =
+        ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
+    lens = lens_input_parameter("river_bank_water__depth")
+    bankfull_depth_2d =
+        ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     bankfull_depth = bankfull_depth_2d[indices]
     zb = bankfull_elevation_2d[indices] - bankfull_depth # river bed elevation
 
     bankfull_volume = bankfull_depth .* river_width .* river_length
-    mannings_n = ncread(
-        dataset,
-        config,
-        "river_water_flow__manning_n_parameter";
-        sel = indices,
-        defaults = 0.036,
-        type = Float,
-    )
+    lens = lens_input_parameter("river_water_flow__manning_n_parameter")
+    mannings_n =
+        ncread(dataset, config, lens; sel = indices, defaults = 0.036, type = Float)
 
     n = length(indices)
-
     # set ghost points for boundary condition (downstream river outlet): river width, bed
     # elevation, manning n is copied from the upstream cell.
     append!(river_length, riverlength_bc)
@@ -151,14 +129,10 @@ end
 "Initialize shallow water river flow model variables"
 function LocalInertialRiverFlowVariables(dataset, config, indices, n_edges, inds_pit)
     floodplain_1d = get(config.model, "floodplain_1d", false)::Bool
-    riverdepth_bc = ncread(
-        dataset,
-        config,
-        "model_boundary_condition~river_bank_water__depth";
-        sel = inds_pit,
-        defaults = 0.0,
-        type = Float,
-    )
+
+    lens = lens_input_parameter("model_boundary_condition~river_bank_water__depth")
+    riverdepth_bc =
+        ncread(dataset, config, lens; sel = inds_pit, defaults = 0.0, type = Float)
 
     n = length(indices)
     # set river depth h to zero (including reservoir and lake locations)
@@ -659,22 +633,11 @@ function LocalInertialOverlandFlowParameters(
 
     @info "Local inertial approach is used for overlandflow." cfl theta h_thresh froude_limit
 
-    mannings_n = ncread(
-        dataset,
-        config,
-        "land_surface_water_flow__manning_n_parameter";
-        sel = indices,
-        defaults = 0.072,
-        type = Float,
-    )
-    elevation_2d = ncread(
-        dataset,
-        config,
-        "land_surface_water_flow__ground_elevation";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+    lens = lens_input_parameter("land_surface_water_flow__manning_n_parameter")
+    mannings_n =
+        ncread(dataset, config, lens; sel = indices, defaults = 0.072, type = Float)
+    lens = lens_input_parameter("land_surface_water_flow__ground_elevation")
+    elevation_2d = ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     elevation = elevation_2d[indices]
     n = length(indices)
 
@@ -876,11 +839,10 @@ function update_boundary_conditions!(
         net_runoff_river .* network.land.area .* 0.001 ./ dt
 
     if !isnothing(reservoir) || !isnothing(lake)
-        inflow_land = get_inflow_waterbody(river, model)
         inflow_subsurface = get_inflow_waterbody(river, subsurface)
 
-        @. inflow_waterbody[network.river_indices] =
-            inflow_land[network.river_indices] + inflow_subsurface[network.river_indices]
+        @. inflow_waterbody[network.river.land_indices] =
+            inflow_subsurface[network.river.land_indices]
     end
     return nothing
 end
@@ -900,16 +862,9 @@ function update!(
 ) where {T}
     (; reservoir, lake) = river.boundary_conditions
 
-    if !isnothing(reservoir)
-        reservoir.boundary_conditions.inflow .= 0.0
-        reservoir.variables.totaloutflow .= 0.0
-        reservoir.variables.actevap .= 0.0
-    end
-    if !isnothing(lake)
-        lake.boundary_conditions.inflow .= 0.0
-        lake.variables.totaloutflow .= 0.0
-        lake.variables.actevap .= 0.0
-    end
+    set_waterbody_vars!(reservoir)
+    set_waterbody_vars!(lake)
+
     river.variables.q_av .= 0.0
     river.variables.h_av .= 0.0
     land.variables.h_av .= 0.0
@@ -930,6 +885,9 @@ function update!(
     river.variables.h_av ./= dt
     land.variables.h_av ./= dt
 
+    average_waterbody_vars!(reservoir, dt)
+    average_waterbody_vars!(lake, dt)
+
     return nothing
 end
 
@@ -1124,14 +1082,9 @@ end
 
 "Initialize floodplain profile `FloodPlainProfile`"
 function FloodPlainProfile(dataset, config, indices; river_width, river_length, index_pit)
-    volume = ncread(
-        dataset,
-        config,
-        "floodplain_water__sum_of_volume-per-depth";
-        sel = indices,
-        type = Float,
-        dimname = :flood_depth,
-    )
+    lens = lens_input_parameter("floodplain_water__sum_of_volume-per-depth")
+    volume =
+        ncread(dataset, config, lens; sel = indices, type = Float, dimname = :flood_depth)
     n = length(indices)
 
     # for convenience (interpolation) flood depth 0.0 m is added, with associated area (a),
@@ -1234,14 +1187,9 @@ function FloodPlainParameters(
     profile =
         FloodPlainProfile(dataset, config, indices; river_width, river_length, index_pit)
 
-    mannings_n = ncread(
-        dataset,
-        config,
-        "lateral.river.floodplain.mannings_n";
-        sel = indices,
-        defaults = 0.072,
-        type = Float,
-    )
+    lens = lens_input_parameter("floodplain_water_flow__manning_n_parameter")
+    mannings_n =
+        ncread(dataset, config, lens; sel = indices, defaults = 0.072, type = Float)
     # manning roughness at edges
     append!(mannings_n, mannings_n[index_pit]) # copy to ghost nodes
     mannings_n_sq = fill(Float(0), n_edges)
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 4ae064330..1cc675b2b 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -41,43 +41,29 @@ function initialize_sbm_gwf_model(config::Config)
 
     dataset = NCDataset(static_path)
 
-    subcatch_2d =
-        ncread(dataset, config, "subcatchment"; optional = false, allow_missing = true)
+    lens = lens_input("subcatchment")
+    subcatch_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     # indices based on catchment
     indices, reverse_indices = active_indices(subcatch_2d, missing)
     n_land_cells = length(indices)
     modelsize_2d = size(subcatch_2d)
 
-    river_location_2d = ncread(
-        dataset,
-        config,
-        "river_location";
-        optional = false,
-        type = Bool,
-        fill = false,
-    )
+    lens = lens_input("river_location")
+    river_location_2d =
+        ncread(dataset, config, lens; optional = false, type = Bool, fill = false)
     river_location = river_location_2d[indices]
-    river_width_2d = ncread(
-        dataset,
-        config,
-        "river__width";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+
+    lens = lens_input_parameter("river__width")
+    river_width_2d = ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     river_width = river_width_2d[indices]
-    river_length_2d = ncread(
-        dataset,
-        config,
-        "river__length";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+
+    lens = lens_input_parameter("river__length")
+    river_length_2d =
+        ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     river_length = river_length_2d[indices]
 
-    altitude =
-        ncread(dataset, config, "altitude"; optional = false, sel = indices, type = Float)
+    lens = lens_input("altitude")
+    altitude = ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
 
     # read x, y coordinates and calculate cell length [m]
     y_coords = read_y_axis(dataset)
@@ -118,17 +104,13 @@ function initialize_sbm_gwf_model(config::Config)
     end
 
     # overland flow (kinematic wave)
-    land_slope = ncread(
-        dataset,
-        config,
-        "land_surface__slope";
-        optional = false,
-        sel = indices,
-        type = Float,
-    )
+    lens = lens_input_parameter("land_surface__slope")
+    land_slope =
+        ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
     clamp!(land_slope, 0.00001, Inf)
-    ldd_2d = ncread(dataset, config, "ldd"; optional = false, allow_missing = true)
 
+    lens = lens_input("ldd")
+    ldd_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     ldd = ldd_2d[indices]
 
     flow_length = map(get_flow_length, ldd, x_length, y_length)
@@ -265,10 +247,10 @@ function initialize_sbm_gwf_model(config::Config)
 
     # drain boundary of unconfined aquifer (optional)
     if do_drains
-        drain_2d =
-            ncread(dataset, config, "lateral.subsurface.drain"; type = Bool, fill = false)
-
+        lens = lens_input_parameter("land_drain_location__flag")
+        drain_2d = ncread(dataset, config, lens; type = Bool, fill = false)
         drain = drain_2d[indices]
+
         # check if drain occurs where overland flow is not possible (surface_flow_width =
         # 0.0) and correct if this is the case
         false_drain = filter(
@@ -508,8 +490,7 @@ function update!(model::Model{N, L, V, R, W, T}) where {N, L, V, R, W, T <: SbmG
         lateral.river.variables.h_av .+ lateral.subsurface.river.parameters.bottom
 
     # determine stable time step for groundwater flow
-    conductivity_profile =
-        get(config.input.lateral.subsurface, "conductivity_profile", "uniform")
+    conductivity_profile = get(config.model, "conductivity_profile", "uniform")
     dt_gw = stable_timestep(aquifer, conductivity_profile) # time step in day (Float64)
     dt_sbm = (dt / tosecond(basetimestep)) # dt is in seconds (Float64)
     if dt_gw < dt_sbm
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 296c31b0b..2072b3fd0 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -40,39 +40,25 @@ function initialize_sbm_model(config::Config)
 
     dataset = NCDataset(static_path)
 
-    subcatch_2d =
-        ncread(dataset, config, "subcatchment"; optional = false, allow_missing = true)
+    lens = lens_input("subcatchment")
+    subcatch_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     # indices based on sub-catchments
     indices, reverse_indices = active_indices(subcatch_2d, missing)
     n_land_cells = length(indices)
     modelsize_2d = size(subcatch_2d)
 
-    river_location_2d = ncread(
-        dataset,
-        config,
-        "river_location";
-        optional = false,
-        type = Bool,
-        fill = false,
-    )
+    lens = lens_input("river_location")
+    river_location_2d =
+        ncread(dataset, config, lens; optional = false, type = Bool, fill = false)
     river_location = river_location_2d[indices]
-    river_width_2d = ncread(
-        dataset,
-        config,
-        "river__width";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+
+    lens = lens_input_parameter("river__width")
+    river_width_2d = ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     river_width = river_width_2d[indices]
-    river_length_2d = ncread(
-        dataset,
-        config,
-        "river__length";
-        optional = false,
-        type = Float,
-        fill = 0,
-    )
+
+    lens = lens_input_parameter("river__length")
+    river_length_2d =
+        ncread(dataset, config, lens; optional = false, type = Float, fill = 0)
     river_length = river_length_2d[indices]
 
     # read x, y coordinates and calculate cell length [m]
@@ -112,22 +98,18 @@ function initialize_sbm_model(config::Config)
         lake = nothing
     end
 
-    ldd_2d = ncread(dataset, config, "ldd"; optional = false, allow_missing = true)
+    lens = lens_input("ldd")
+    ldd_2d = ncread(dataset, config, lens; optional = false, allow_missing = true)
     ldd = ldd_2d[indices]
     if do_pits
-        pits_2d =
-            ncread(dataset, config, "pits"; optional = false, type = Bool, fill = false)
+        lens = lens_input("pits")
+        pits_2d = ncread(dataset, config, lens; optional = false, type = Bool, fill = false)
         ldd = set_pit_ldd(pits_2d, ldd, indices)
     end
 
-    land_slope = ncread(
-        dataset,
-        config,
-        "land_surface__slope";
-        optional = false,
-        sel = indices,
-        type = Float,
-    )
+    lens = lens_input_parameter("land_surface__slope")
+    land_slope =
+        ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
     clamp!(land_slope, 0.00001, Inf)
     flow_length = map(get_flow_length, ldd, x_length, y_length)
     flow_width = (x_length .* y_length) ./ flow_length
@@ -147,7 +129,11 @@ function initialize_sbm_model(config::Config)
             y_length,
         )
         # update variables `ssf`, `ssfmax` and `kh` (layered profile) based on ksat_profile
-        kh_profile_type = get(config.model, "saturated_hydraulic_conductivity_profile", "exponential")::String
+        kh_profile_type = get(
+            config.model,
+            "saturated_hydraulic_conductivity_profile",
+            "exponential",
+        )::String
         if kh_profile_type == "exponential" || kh_profile_type == "exponential_constant"
             initialize_lateralssf!(subsurface_flow, subsurface_flow.parameters.kh_profile)
         elseif kh_profile_type == "layered" || kh_profile_type == "layered_exponential"
diff --git a/src/snow/snow.jl b/src/snow/snow.jl
index 9707fac9e..51d6b0ae1 100644
--- a/src/snow/snow.jl
+++ b/src/snow/snow.jl
@@ -69,47 +69,21 @@ struct NoSnowModel{T} <: AbstractSnowModel{T} end
 
 "Initialize snow HBV model parameters"
 function SnowHbvParameters(dataset, config, indices, dt)
+    lens = lens_input_parameter("snowpack__degree-day_coefficient")
     cfmax =
-        ncread(
-            dataset,
-            config,
-            "snowpack__degree-day_coefficient";
-            sel = indices,
-            defaults = 3.75653,
-            type = Float,
-        ) .* (dt / basetimestep)
-    tt = ncread(
-        dataset,
-        config,
-        "atmosphere_air__snowfall_temperature_threshold";
-        sel = indices,
-        defaults = 0.0,
-        type = Float,
-    )
-    tti = ncread(
-        dataset,
-        config,
-        "atmosphere_air__snowfall_temperature_interval";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
-    )
-    ttm = ncread(
-        dataset,
-        config,
-        "snowpack__melting_temperature_threshold";
-        sel = indices,
-        defaults = 0.0,
-        type = Float,
-    )
-    whc = ncread(
-        dataset,
-        config,
-        "snowpack__liquid_water_holding_capacity";
-        sel = indices,
-        defaults = 0.1,
-        type = Float,
-    )
+        ncread(dataset, config, lens; sel = indices, defaults = 3.75653, type = Float) .*
+        (dt / basetimestep)
+    lens = lens_input_parameter("atmosphere_air__snowfall_temperature_threshold")
+    tt = ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float)
+
+    lens = lens_input_parameter("atmosphere_air__snowfall_temperature_interval")
+    tti = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+
+    lens = lens_input_parameter("snowpack__melting_temperature_threshold")
+    ttm = ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float)
+
+    lens = lens_input_parameter("snowpack__liquid_water_holding_capacity")
+    whc = ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
     snow_hbv_params =
         SnowHbvParameters(; cfmax = cfmax, tt = tt, tti = tti, ttm = ttm, whc = whc)
     return snow_hbv_params
diff --git a/src/soil/soil.jl b/src/soil/soil.jl
index 5fb8d38a8..36266e8cf 100644
--- a/src/soil/soil.jl
+++ b/src/soil/soil.jl
@@ -211,48 +211,43 @@ function sbm_kv_profiles(
     sumlayers,
     dt,
 )
-    kv_profile_type = get(config.model, "saturated_hydraulic_conductivity_profile", "exponential")::String
+    kv_profile_type =
+        get(config.model, "saturated_hydraulic_conductivity_profile", "exponential")::String
     n = length(indices)
     if kv_profile_type == "exponential"
         kv_profile = KvExponential(kv_0, f)
     elseif kv_profile_type == "exponential_constant"
-        z_exp = ncread(
-            dataset,
-            config,
-            "soil_vertical_saturated_hydraulic_conductivity_profile~exponential_below-surface__depth";
-            optional = false,
-            sel = indices,
-            type = Float,
+        lens = lens_input_parameter(
+            "soil_vertical_saturated_hydraulic_conductivity_profile~exponential_below-surface__depth",
         )
+        z_exp = ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
         exp_profile = KvExponential(kv_0, f)
         kv_profile = KvExponentialConstant(exp_profile, z_exp)
     elseif kv_profile_type == "layered" || kv_profile_type == "layered_exponential"
+        lens = lens_input_parameter("soil_water__vertical_saturated_hydraulic_conductivity")
         kv =
             ncread(
                 dataset,
                 config,
-                "soil_water__vertical_saturated_hydraulic_conductivity";
+                lens;
                 sel = indices,
                 defaults = 1000.0,
                 type = Float,
                 dimname = :layer,
             ) .* (dt / basetimestep)
         if size(kv, 1) != maxlayers
-            parname = param(config.input, "kv")
+            parname = lens(config)
             size1 = size(kv, 1)
             error("$parname needs a layer dimension of size $maxlayers, but is $size1")
         end
         if kv_profile_type == "layered"
             kv_profile = KvLayered(svectorscopy(kv, Val{maxlayers}()))
         else
-            z_layered = ncread(
-                dataset,
-                config,
-                "soil_vertical_saturated_hydraulic_conductivity_profile~layered_below-surface__depth";
-                optional = false,
-                sel = indices,
-                type = Float,
+            lens = lens_input_parameter(
+                "soil_vertical_saturated_hydraulic_conductivity_profile~layered_below-surface__depth",
             )
+            z_layered =
+                ncread(dataset, config, lens; optional = false, sel = indices, type = Float)
             nlayers_kv = fill(0, n)
             for i in eachindex(nlayers_kv)
                 layers = @view sumlayers[i][2:nlayers[i]]
@@ -343,6 +338,7 @@ end
 "Initialize SBM soil model parameters"
 function SbmSoilParameters(dataset, config, vegetation_parameter_set, indices, dt)
     config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
+
     if length(config_thicknesslayers) > 0
         thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv)))
         cum_depth_layers = pushfirst(cumsum(thicknesslayers), 0.0)
@@ -352,212 +348,121 @@ function SbmSoilParameters(dataset, config, vegetation_parameter_set, indices, d
         cum_depth_layers = pushfirst(cumsum(thicknesslayers), 0.0)
         maxlayers = 1
     end
+
+    lens = lens_input_parameter("soil_surface_temperature__weight_coefficient")
     w_soil =
-        ncread(
-            dataset,
-            config,
-            "soil_surface_temperature__weight_coefficient";
-            sel = indices,
-            defaults = 0.1125,
-            type = Float,
-        ) .* (dt / basetimestep)
-    cf_soil = ncread(
-        dataset,
-        config,
-        "soil_surface_water__infiltration_reduction_parameter";
-        sel = indices,
-        defaults = 0.038,
-        type = Float,
-    )
+        ncread(dataset, config, lens; sel = indices, defaults = 0.1125, type = Float) .*
+        (dt / basetimestep)
+
+    lens = lens_input_parameter("soil_surface_water__infiltration_reduction_parameter")
+    cf_soil = ncread(dataset, config, lens; sel = indices, defaults = 0.038, type = Float)
+
     # soil parameters
-    theta_s = ncread(
-        dataset,
-        config,
-        "soil_water__saturated_volume_fraction";
-        sel = indices,
-        defaults = 0.6,
-        type = Float,
-    )
-    theta_r = ncread(
-        dataset,
-        config,
-        "soil_water__residual_volume_fraction";
-        sel = indices,
-        defaults = 0.01,
-        type = Float,
+    lens = lens_input_parameter("soil_water__saturated_volume_fraction")
+    theta_s = ncread(dataset, config, lens; sel = indices, defaults = 0.6, type = Float)
+
+    lens = lens_input_parameter("soil_water__residual_volume_fraction")
+    theta_r = ncread(dataset, config, lens; sel = indices, defaults = 0.01, type = Float)
+
+    lens = lens_input_parameter(
+        "soil_surface_water__vertical_saturated_hydraulic_conductivity",
     )
     kv_0 =
-        ncread(
-            dataset,
-            config,
-            "soil_surface_water__vertical_saturated_hydraulic_conductivity";
-            sel = indices,
-            defaults = 3000.0,
-            type = Float,
-        ) .* (dt / basetimestep)
-    f = ncread(
-        dataset,
-        config,
-        "soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter";
-        sel = indices,
-        defaults = 0.001,
-        type = Float,
-    )
-    hb = ncread(
-        dataset,
-        config,
-        "soil_water__air_entry_pressure_head";
-        sel = indices,
-        defaults = -10.0,
-        type = Float,
-    )
-    h1 = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~1";
-        sel = indices,
-        defaults = 0.0,
-        type = Float,
-    )
-    h2 = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~2";
-        sel = indices,
-        defaults = -100.0,
-        type = Float,
-    )
-    h3_high = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~3~high";
-        sel = indices,
-        defaults = -400.0,
-        type = Float,
-    )
-    h3_low = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~3~low";
-        sel = indices,
-        defaults = -1000.0,
-        type = Float,
-    )
-    h4 = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~4";
-        sel = indices,
-        defaults = -15849.0,
-        type = Float,
-    )
-    alpha_h1 = ncread(
-        dataset,
-        config,
-        "vegetation_root__feddes_critial_pressure_head_h~1_reduction_coefficient";
-        sel = indices,
-        defaults = 1.0,
-        type = Float,
+        ncread(dataset, config, lens; sel = indices, defaults = 3000.0, type = Float) .*
+        (dt / basetimestep)
+
+    lens = lens_input_parameter(
+        "soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter",
     )
-    soilthickness = ncread(
-        dataset,
-        config,
-        "soil__thickness";
-        sel = indices,
-        defaults = 2000.0,
-        type = Float,
+    f = ncread(dataset, config, lens; sel = indices, defaults = 0.001, type = Float)
+
+    lens = lens_input_parameter("soil_water__air_entry_pressure_head")
+    hb = ncread(dataset, config, lens; sel = indices, defaults = -10.0, type = Float)
+
+    lens = lens_input_parameter("vegetation_root__feddes_critial_pressure_head_h~1")
+    h1 = ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float)
+
+    lens = lens_input_parameter("vegetation_root__feddes_critial_pressure_head_h~2")
+    h2 = ncread(dataset, config, lens; sel = indices, defaults = -100.0, type = Float)
+
+    lens = lens_input_parameter("vegetation_root__feddes_critial_pressure_head_h~3~high")
+    h3_high = ncread(dataset, config, lens; sel = indices, defaults = -400.0, type = Float)
+
+    lens = lens_input_parameter("vegetation_root__feddes_critial_pressure_head_h~3~low")
+    h3_low = ncread(dataset, config, lens; sel = indices, defaults = -1000.0, type = Float)
+
+    lens = lens_input_parameter("vegetation_root__feddes_critial_pressure_head_h~4")
+    h4 = ncread(dataset, config, lens; sel = indices, defaults = -15849.0, type = Float)
+
+    lens = lens_input_parameter(
+        "vegetation_root__feddes_critial_pressure_head_h~1_reduction_coefficient",
     )
+    alpha_h1 = ncread(dataset, config, lens; sel = indices, defaults = 1.0, type = Float)
+
+    lens = lens_input_parameter("soil__thickness")
+    soilthickness =
+        ncread(dataset, config, lens; sel = indices, defaults = 2000.0, type = Float)
+
+    lens = lens_input_parameter("soil~compacted_surface_water__infiltration_capacity")
     infiltcappath =
-        ncread(
-            dataset,
-            config,
-            "soil~compacted_surface_water__infiltration_capacity";
-            sel = indices,
-            defaults = 10.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 10.0, type = Float) .*
+        (dt / basetimestep)
+
+    lens = lens_input_parameter("soil~non-compacted_surface_water__infiltration_capacity")
     infiltcapsoil =
-        ncread(
-            dataset,
-            config,
-            "soil~non-compacted_surface_water__infiltration_capacity";
-            sel = indices,
-            defaults = 100.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 100.0, type = Float) .*
+        (dt / basetimestep)
+
+    lens = lens_input_parameter("soil_water_sat-zone_bottom__max_leakage_volume_flux")
     maxleakage =
-        ncread(
-            dataset,
-            config,
-            "soil_water_sat-zone_bottom__max_leakage_volume_flux";
-            sel = indices,
-            defaults = 0.0,
-            type = Float,
-        ) .* (dt / basetimestep)
+        ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float) .*
+        (dt / basetimestep)
 
+    lens = lens_input_parameter("soil_water__brooks-corey_epsilon_parameter")
     c = ncread(
         dataset,
         config,
-        "soil_water__brooks-corey_epsilon_parameter";
+        lens;
         sel = indices,
         defaults = 10.0,
         type = Float,
         dimname = :layer,
     )
     if size(c, 1) != maxlayers
-        parname = param(config.input.vertical, "c")
+        parname = lens(config)
         size1 = size(c, 1)
         error("$parname needs a layer dimension of size $maxlayers, but is $size1")
     end
+
+    lens =
+        lens_input_parameter("soil_water__vertical_saturated_hydraulic_conductivity_factor")
     kvfrac = ncread(
         dataset,
         config,
-        "soil_water__vertical_saturated_hydraulic_conductivity_factor";
+        lens;
         sel = indices,
         defaults = 1.0,
         type = Float,
         dimname = :layer,
     )
     if size(kvfrac, 1) != maxlayers
-        parname = param(config.input, "soil_water__vertical_saturated_hydraulic_conductivity_factor")
+        parname = lens(config)
         size1 = size(kvfrac, 1)
         error("$parname needs a layer dimension of size $maxlayers, but is $size1")
     end
     # fraction compacted area
-    pathfrac = ncread(
-        dataset,
-        config,
-        "soil~compacted__area_fraction";
-        sel = indices,
-        defaults = 0.01,
-        type = Float,
-    )
+    lens = lens_input_parameter("soil~compacted__area_fraction")
+    pathfrac = ncread(dataset, config, lens; sel = indices, defaults = 0.01, type = Float)
 
     # vegetation parameters
-    rootdistpar = ncread(
-        dataset,
-        config,
-        "soil_root~wet__sigmoid_function_shape_parameter";
-        sel = indices,
-        defaults = -500.0,
-        type = Float,
-    )
-    cap_hmax = ncread(
-        dataset,
-        config,
-        "soil_capillary-rise__max_water-table_depth";
-        sel = indices,
-        defaults = 2000.0,
-        type = Float,
-    )
-    cap_n = ncread(
-        dataset,
-        config,
-        "soil_capillary-rise__averianov_exponent";
-        sel = indices,
-        defaults = 2.0,
-        type = Float,
-    )
+    lens = lens_input_parameter("soil_root~wet__sigmoid_function_shape_parameter")
+    rootdistpar =
+        ncread(dataset, config, lens; sel = indices, defaults = -500.0, type = Float)
+    lens = lens_input_parameter("soil_capillary-rise__max_water-table_depth")
+    cap_hmax = ncread(dataset, config, lens; sel = indices, defaults = 2000.0, type = Float)
+
+    lens = lens_input_parameter("soil_capillary-rise__averianov_exponent")
+    cap_n = ncread(dataset, config, lens; sel = indices, defaults = 2.0, type = Float)
 
     act_thickl = set_layerthickness.(soilthickness, (cum_depth_layers,), (thicknesslayers,))
     sumlayers = @. pushfirst(cumsum(act_thickl), 0.0)
@@ -566,11 +471,12 @@ function SbmSoilParameters(dataset, config, vegetation_parameter_set, indices, d
     if length(config_thicknesslayers) > 0
         # root fraction read from dataset file, in case of multiple soil layers and TOML file
         # includes "vertical.rootfraction"
-        if haskey(config.input, "soil_layer_root__length_density_fraction")
+        if haskey(config.input.parameters, "soil_layer_root__length_density_fraction")
+            lens = lens_input_parameter("soil_layer_root__length_density_fraction")
             rootfraction = ncread(
                 dataset,
                 config,
-                "soil_layer_root__length_density_fraction";
+                lens;
                 sel = indices,
                 optional = false,
                 type = Float,
diff --git a/src/standard_name.jl b/src/standard_name.jl
index a1dc9c7c8..af2e4626e 100644
--- a/src/standard_name.jl
+++ b/src/standard_name.jl
@@ -1,5 +1,71 @@
 
-const standard_name_map = Dict{String,String}(
-    "vegetation__leaf-area_index" => "vertical.vegetation_parameter_set.leaf_area_index",
-    "river_water__volume_inflow_rate"=> "lateral.river.boundary_conditions.inflow",
+const standard_name_map = Dict{String, ComposedFunction}(
+    "atmosphere_water__precipitation_volume_flux" =>
+        @optic(_.vertical.atmospheric_forcing.precipitation),
+    "land_surface_water__potential_evaporation_volume_flux" =>
+        @optic(_.vertical.atmospheric_forcing.potential_evaporation),
+    "atmosphere_air__temperature" => @optic(_.vertical.atmospheric_forcing.temperature),
+    "vegetation__leaf-area_index" =>
+        @optic(_.vertical.vegetation_parameter_set.leaf_area_index),
+    "vegetation_canopy_water__storage" =>
+        @optic(_.vertical.interception.variables.canopy_storage),
+    "river_water__volume_inflow_rate" =>
+        @optic(_.lateral.river.boundary_conditions.inflow),
+    "river_water__volume_flow_rate" => @optic(_.lateral.river.variables.q),
+    "river_water__time_average_of_volume_flow_rate" =>
+        @optic(_.lateral.river.variables.q_av),
+    "river_water__depth" => @optic(_.lateral.river.variables.h),
+    "river_water__time_average_of_depth" => @optic(_.lateral.river.variables.h_av),
+    "river_water__volume" => @optic(_.lateral.river.variables.volume),
+    "floodplain_water__volume" => @optic(_.lateral.river.floodplain.variables.volume),
+    "floodplain_water__depth" => @optic(_.lateral.river.floodplain.variables.h),
+    "reservoir_water__volume" =>
+        @optic(_.lateral.river.boundary_conditions.reservoir.variables.volume),
+    "soil_water_sat-zone_top__recharge_volume_flux" =>
+        @optic(_.vertical.soil.variables.recharge),
+    "soil_water_unsat-zone__depth-per-soil_layer" =>
+        @optic(_.vertical.soil.variables.ustorelayerdepth),
+    "soil_water_sat-zone__depth" => @optic(_.vertical.soil.variables.satwaterdepth),
+    "soil_water_sat-zone_top__depth" => @optic(_.vertical.soil.variables.zi),
+    "soil_surface__temperature" => @optic(_.vertical.soil.variables.tsoil),
+    "subsurface_water__volume_flow_rate" => @optic(_.lateral.subsurface.variables.ssf),
+    "snowpack~liquid__depth" => @optic(_.vertical.snow.variables.snow_water),
+    "snowpack~dry__leq-depth" => @optic(_.vertical.snow.variables.snow_storage),
+    "glacier_ice__leq-volume" => @optic(_.vertical.glacier.variables.glacier_store),
+    "land_surface_water__volume_flow_rate" => @optic(_.lateral.land.variables.q),
+    "land_surface_water__depth" => @optic(_.lateral.land.variables.h),
+    "land_surface_water__time_average_of_depth" =>
+        @optic(_.lateral.land.variables.h_av),
+    "land_surface_water__volume" => @optic(_.lateral.land.variables.volume),
+    "land_surface_water__x_component_of_volume_flow_rate" =>
+        @optic(_.lateral.land.variables.qx),
+    "land_surface_water__y_component_of_volume_flow_rate" =>
+        @optic(_.lateral.land.variables.qy),
+    "land_surface_water~paddy__depth" => @optic(_.vertical.demand.paddy.variables.h),
+    "land~domestic__gross_water_demand_flux" =>
+        @optic(_.vertical.demand.domestic.demand.demand_gross),
+    "land~domestic__net_water_demand_flux" =>
+        @optic(_.vertical.demand.domestic.demand.demand_net),
+    "land~industry__gross_water_demand_flux" =>
+        @optic(_.vertical.demand.industry.demand.demand_gross),
+    "land~industry__net_water_demand_flux" =>
+        @optic(_.vertical.demand.industry.demand.demand_net),
+    "land~livestock__gross_water_demand_flux" =>
+        @optic(_.vertical.demand.livestock.demand.demand_gross),
+    "land~livestock__net_water_demand_flux" =>
+        @optic(_.vertical.demand.livestock.demand.demand_net),
+    "land~irrigated-paddy__irrigation_trigger_flag" =>
+        @optic(_.vertical.demand.paddy.parameters.irrigation_trigger),
+    "land~irrigated-non-paddy__irrigation_trigger_flag" =>
+        @optic(_.vertical.demand.nonpaddy.parameters.irrigation_trigger),
+    "land~irrigated__allocated_water_volume_flux" =>
+        @optic(_.vertical.allocation.variables.irri_alloc),
+    "subsurface_water__hydraulic_head" =>
+        @optic(_.lateral.subsurface.flow.aquifer.variables.head),
+    "subsurface_water_sat-zone_top__net-recharge_volume_flux" =>
+        @optic(_.lateral.subsurface.recharge.variables.flux),
+    "land_drain_water~to-subsurface__volume_flow_rate" =>
+        @optic(_.lateral.subsurface.drain.variables.flux),
+    "river_water~to-subsurface__volume_flow_rate" =>
+        @optic(_.lateral.subsurface.river.variables.flux),
 )
\ No newline at end of file
diff --git a/src/surfacewater/runoff.jl b/src/surfacewater/runoff.jl
index feb71eb1b..98ee027b8 100644
--- a/src/surfacewater/runoff.jl
+++ b/src/surfacewater/runoff.jl
@@ -36,14 +36,8 @@ end
 "Initialize open water runoff parameters"
 function OpenWaterRunoffParameters(dataset, config, indices, riverfrac)
     # fraction open water
-    waterfrac = ncread(
-        dataset,
-        config,
-        "land~water-covered__area_fraction";
-        sel = indices,
-        defaults = 0.0,
-        type = Float,
-    )
+    lens = lens_input_parameter("land~water-covered__area_fraction")
+    waterfrac = ncread(dataset, config, lens; sel = indices, defaults = 0.0, type = Float)
     waterfrac = max.(waterfrac .- riverfrac, Float(0.0))
     params = OpenWaterRunoffParameters(; waterfrac = waterfrac, riverfrac = riverfrac)
     return params
diff --git a/src/utils.jl b/src/utils.jl
index 2599b3199..4bcd013dd 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -77,14 +77,14 @@ function active_indices(subcatch_2d::AbstractMatrix, nodata)
     return indices, reverse_indices
 end
 
-function active_indices(network::NamedTuple, key::Tuple)
-    if :reservoir in key
+function active_indices(network::NamedTuple, key::AbstractString)
+    if occursin("reservoir", key)
         return network.reservoir.indices_outlet
-    elseif :lake in key
+    elseif occursin("lake", key)
         return network.lake.indices_outlet
-    elseif :river in key
+    elseif occursin("river", key)
         return network.river.indices
-    elseif :drain in key
+    elseif occursin("drain", key)
         return network.drain.indices
     else
         return network.land.indices
@@ -157,7 +157,9 @@ function set_states!(instate_path, model; type = nothing, dimname = nothing)
     (; network, config) = model
 
     # Check if required states are covered
-    state_ncnames = check_states(config)
+    # TODO: revert back to state checking
+    # state_ncnames = check_states(config)
+    state_ncnames = ncnames(config.state.variables)
 
     # states in netCDF include dim time (one value) at index 3 or 4, 3 or 4 dims are allowed
     NCDataset(instate_path) do ds
@@ -188,7 +190,8 @@ function set_states!(instate_path, model; type = nothing, dimname = nothing)
                     end
                 end
                 # set state in model object
-                param(model, state) .= svectorscopy(A, Val{size(A)[1]}())
+                lens = standard_name_map[state]
+                lens(model) .= svectorscopy(A, Val{size(A)[1]}())
                 # 3 dims (x,y,time)
             elseif dims == 3
                 A = read_standardized(ds, ncname, (x = :, y = :, time = 1))
@@ -201,7 +204,8 @@ function set_states!(instate_path, model; type = nothing, dimname = nothing)
                     end
                 end
                 # set state in model object, only set active cells ([1:n]) (ignore boundary conditions/ghost points)
-                param(model, state)[1:n] .= A
+                lens = standard_name_map[state]
+                lens(model)[1:n] .= A
             else
                 error(
                     "Number of state dims should be 3 or 4, number of dims = ",
@@ -239,7 +243,7 @@ arguments to get selections of data in desired types, with or without missing va
 function ncread(
     nc,
     config::Config,
-    parameter::AbstractString;
+    parameter::NamedTuple;
     alias = nothing,
     optional = true,
     sel = nothing,
@@ -253,12 +257,13 @@ function ncread(
     # if var has type Config, input parameters can be changed.
     if isnothing(alias)
         if optional
-            var = param(config.input, parameter, nothing)
+            var = _lens(config, parameter.lens, nothing)
         else
-            var = param(config.input, parameter)
+            var = parameter.lens(config)
+            var = var isa AbstractDict ? Config(var, pathof(config)) : var
         end
     else
-        var = get_alias(config.input, parameter, alias, nothing)
+        var = get_alias(config, parameter.lens, alias, nothing)
     end
 
     # dim `time` is also included in `dim_sel`: this allows for cyclic parameters (read
@@ -274,7 +279,7 @@ function ncread(
     end
 
     if isnothing(var)
-        @info "Set `$parameter` using default value `$defaults`."
+        @info "Set `$(parameter.name)` using default value `$defaults`."
         @assert !isnothing(defaults)
         if !isnothing(type)
             defaults = convert(type, defaults)
@@ -292,7 +297,7 @@ function ncread(
     var, mod = ncvar_name_modifier(var; config = config)
 
     if !isnothing(mod.value)
-        @info "Set `$parameter` using default value `$(mod.value)`."
+        @info "Set `$(parameter.name)` using default value `$(mod.value)`."
         if isnothing(dimname)
             return Base.fill(mod.value, length(sel))
             # set to one uniform value
@@ -304,7 +309,7 @@ function ncread(
             return repeat(mod.value, 1, length(sel))
         end
     else
-        @info "Set `$parameter` using netCDF variable `$var`."
+        @info "Set `$(parameter.name)` using netCDF variable `$var`."
         A = read_standardized(nc, var, dim_sel)
         if !isnothing(mod.index)
             # the modifier index is only set in combination with scale and offset for SVectors,
@@ -357,6 +362,9 @@ function ncread(
     return A
 end
 
+lens_input_parameter(p::AbstractString) = (name = p, lens = @optic(_.input.parameters[p]))
+lens_input(p::AbstractString) = (name = p, lens = @optic(_.input[p]))
+
 """
     set_layerthickness(reference_depth::Real, cum_depth::SVector, thickness::SVector)
 
diff --git a/src/vegetation/canopy.jl b/src/vegetation/canopy.jl
index d8af07f08..ef349703d 100644
--- a/src/vegetation/canopy.jl
+++ b/src/vegetation/canopy.jl
@@ -40,14 +40,10 @@ end
 
 "Initialize Gash interception model"
 function GashInterceptionModel(dataset, config, indices, vegetation_parameter_set)
-    e_r = ncread(
-        dataset,
-        config,
-        "vegetation_canopy_water__mean_evaporation-to-mean_precipitation_ratio";
-        sel = indices,
-        defaults = 0.1,
-        type = Float,
+    lens = lens_input_parameter(
+        "vegetation_canopy_water__mean_evaporation-to-mean_precipitation_ratio",
     )
+    e_r = ncread(dataset, config, lens; sel = indices, defaults = 0.1, type = Float)
     n = length(indices)
     params =
         GashParameters(; e_r = e_r, vegetation_parameter_set = vegetation_parameter_set)
diff --git a/test/io.jl b/test/io.jl
index a905b9a53..95500a532 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -21,16 +21,22 @@ config = Wflow.Config(tomlpath)
     @test config.endtime === DateTime(2000, 2)
     @test config.output.path == "output_moselle.nc"
     @test config.output isa Wflow.Config
-    @test collect(keys(config.output)) == ["lateral", "vertical", "path"]
+    @test collect(keys(config.output)) == ["variables", "path"]
 
-    # test removal of key with pop!
-    val = pop!(config.input, "soil_water__saturated_volume_fraction")
-    @test val == "thetaS"
-    @test_throws KeyError config.input.soil_water__saturated_volume_fraction
-    config.input.soil_water__saturated_volume_fraction = "thetaS"
+    # test if soil_water__saturated_volume_fraction can also be provided under the alias theta_s
+    lens = @optic(_.input.parameters.soil_water__saturated_volume_fraction)
+    lens_alias = @optic(_.input.parameters.theta_s)
+    @test Wflow.get_alias(config, lens, lens_alias, nothing) == "thetaS"
+    val = pop!(config.input.parameters, "soil_water__saturated_volume_fraction")
+    config.input.parameters["theta_s"] = val
+    @test Wflow.get_alias(config, lens, lens_alias, nothing) == "thetaS" == "thetaS"
+    config.input.parameters.soil_water__saturated_volume_fraction = "thetaS"
 
     # modifiers can also be applied
-    kvconf = Wflow.param(config.input, "soil_surface_water__vertical_saturated_hydraulic_conductivity", nothing)
+    parameter = Wflow.lens_input_parameter(
+        "soil_surface_water__vertical_saturated_hydraulic_conductivity",
+    )
+    kvconf = Wflow._lens(config, parameter.lens, nothing)
     @test kvconf isa Wflow.Config
     ncname, modifier = Wflow.ncvar_name_modifier(kvconf; config = config)
     @test ncname === "KsatVer"
@@ -187,7 +193,7 @@ model = Wflow.initialize_sbm_model(config)
 Wflow.advance!(model.clock)
 Wflow.load_dynamic_input!(model)
 
-(; vertical, clock, reader, writer) = model
+(; clock, reader, writer) = model
 
 @testset "output and state names" begin
     ncdims = ("lon", "lat", "layer", "time")
@@ -200,29 +206,35 @@ Wflow.load_dynamic_input!(model)
 end
 
 # get a default value if the parameter does not exist
-@test Wflow.param(model, "lateral.doesnt_exist", -1) == -1
+lens = @optic(_.input.parameters.doesnt_exist)
+@test Wflow._lens(config, lens, -1) == -1
 
 @testset "warm states" begin
-    @test Wflow.param(
-        model,
-        "lateral.river.boundary_conditions.reservoir.variables.volume",
-    )[1] ≈ 3.2807224993363418e7
-    @test Wflow.param(model, "vertical.soil.variables.satwaterdepth")[9115] ≈
+    @test Wflow.standard_name_map["reservoir_water__volume"](model)[1] ≈
+          3.2807224993363418e7
+    @test Wflow.standard_name_map["soil_water_sat-zone__depth"](model)[9115] ≈
           477.13548089422125
-    @test Wflow.param(model, "vertical.snow.variables.snow_storage")[5] ≈ 11.019233179897599
-    @test Wflow.param(model, "vertical.soil.variables.tsoil")[5] ≈ 0.21814478119608938
-    @test Wflow.param(model, "vertical.soil.variables.ustorelayerdepth")[50063][1] ≈
+    @test Wflow.standard_name_map["snowpack~dry__leq-depth"](model)[5] ≈ 11.019233179897599
+    @test Wflow.standard_name_map["soil_surface__temperature"](model)[5] ≈
+          0.21814478119608938
+    @test Wflow.standard_name_map["soil_water_unsat-zone__depth-per-soil_layer"](model)[50063][1] ≈
           9.969116007201725
-    @test Wflow.param(model, "vertical.snow.variables.snow_water")[5] ≈ 0.0
-    @test Wflow.param(model, "vertical.interception.variables.canopy_storage")[50063] ≈ 0.0
-    @test Wflow.param(model, "vertical.soil.variables.zi")[50063] ≈ 296.8028609104624
-    @test Wflow.param(model, "lateral.subsurface.variables.ssf")[10606] ≈ 39.972334552895816
-    @test Wflow.param(model, "lateral.river.variables.q")[149] ≈ 53.48673634956338
-    @test Wflow.param(model, "lateral.river.variables.h")[149] ≈ 1.167635369628945
-    @test Wflow.param(model, "lateral.river.variables.volume")[149] ≈ 63854.60119358985
-    @test Wflow.param(model, "lateral.land.variables.q")[2075] ≈ 3.285909284322251
-    @test Wflow.param(model, "lateral.land.variables.h")[2075] ≈ 0.052076262033771775
-    @test Wflow.param(model, "lateral.land.variables.volume")[2075] ≈ 29920.754983235012
+    @test Wflow.standard_name_map["snowpack~liquid__depth"](model)[5] ≈ 0.0
+    @test Wflow.standard_name_map["vegetation_canopy_water__storage"](model)[50063] ≈ 0.0
+    @test Wflow.standard_name_map["soil_water_sat-zone_top__depth"](model)[50063] ≈
+          296.8028609104624
+    @test Wflow.standard_name_map["subsurface_water__volume_flow_rate"](model)[10606] ≈
+          39.972334552895816
+    @test Wflow.standard_name_map["river_water__volume_flow_rate"](model)[149] ≈
+          53.48673634956338
+    @test Wflow.standard_name_map["river_water__depth"](model)[149] ≈ 1.167635369628945
+    @test Wflow.standard_name_map["river_water__volume"](model)[149] ≈ 63854.60119358985
+    @test Wflow.standard_name_map["land_surface_water__volume_flow_rate"](model)[2075] ≈
+          3.285909284322251
+    @test Wflow.standard_name_map["land_surface_water__depth"](model)[2075] ≈
+          0.052076262033771775
+    @test Wflow.standard_name_map["land_surface_water__volume"](model)[2075] ≈
+          29920.754983235012
 end
 
 @testset "reducer" begin
@@ -256,15 +268,15 @@ end
           [9.152995289601465, 8.919674421902961, 8.70537452585209, 8.690681062890977]
 end
 
-config.input["snowpack__degree-day_coefficient"] = Dict("value" => 2.0)
-config.input.soil__thickness = Dict(
+config.input.parameters["snowpack__degree-day_coefficient"] = Dict("value" => 2.0)
+config.input.parameters.soil__thickness = Dict(
     "scale" => 3.0,
     "offset" => 100.0,
     "netcdf" => Dict("variable" => Dict("name" => "SoilThickness")),
 )
-config.input.vertical.atmospheric_forcing.precipitation =
+config.input.forcing.atmosphere_water__precipitation_volume_flux =
     Dict("scale" => 1.5, "netcdf" => Dict("variable" => Dict("name" => "precip")))
-config.input["soil_water__brooks-corey_epsilon_parameter"] = Dict(
+config.input.parameters["soil_water__brooks-corey_epsilon_parameter"] = Dict(
     "scale" => [2.0, 3.0],
     "offset" => [0.0, 0.0],
     "layer" => [1, 3],
@@ -396,7 +408,7 @@ end
 
     # Final run to test error handling during simulation
     tomlpath_error = joinpath(@__DIR__, "sbm_simple-error.toml")
-    config.input.river__width = Dict(
+    config.input.parameters.river__width = Dict(
         "scale" => 0.0,
         "offset" => 0.0,
         "netcdf" => Dict("variable" => Dict("name" => "wflow_riverwidth")),
@@ -437,7 +449,7 @@ end
     @test clock.time == DateTimeNoLeap(2000, 3, 1)
 end
 
-@testset "State checking" begin
+#= @testset "State checking" begin
     tomlpath = joinpath(@__DIR__, "sbm_config.toml")
     config = Wflow.Config(tomlpath)
 
@@ -481,4 +493,4 @@ end
     @test (:lateral, :river, :variables, :q) in required_states
     @test (:lateral, :river, :variables, :h_av) in required_states
     @test (:lateral, :land, :variables, :h_av) in required_states
-end
+end =#
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 3f60ef77e..e46061760 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -178,14 +178,14 @@ end
 tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 config = Wflow.Config(tomlpath)
 
-config.input.vertical.atmospheric_forcing.precipitation =
+config.input.forcing.atmosphere_water__precipitation_volume_flux =
     Dict("scale" => 2.0, "netcdf" => Dict("variable" => Dict("name" => "precip")))
-config.input.vertical.atmospheric_forcing.potential_evaporation = Dict(
+config.input.forcing.land_surface_water__potential_evaporation_volume_flux = Dict(
     "scale" => 3.0,
     "offset" => 1.50,
     "netcdf" => Dict("variable" => Dict("name" => "pet")),
 )
-config.input["vegetation__leaf-area_index"] =
+config.input.parameters["vegetation__leaf-area_index"] =
     Dict("scale" => 1.6, "netcdf" => Dict("variable" => Dict("name" => "LAI")))
 
 model = Wflow.initialize_sbm_model(config)
@@ -207,11 +207,9 @@ end
 tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 config = Wflow.Config(tomlpath)
 
-config.input.cyclic = [
-    "vegetation__leaf-area_index",
-    "river_water__volume_inflow_rate",
-]
-config.input.river_water__volume_inflow_rate = "inflow"
+config.input.dynamic.cyclic =
+    ["vegetation__leaf-area_index", "river_water__volume_inflow_rate"]
+config.input.parameters.river_water__volume_inflow_rate = "inflow"
 
 model = Wflow.initialize_sbm_model(config)
 Wflow.run_timestep!(model)
@@ -224,7 +222,7 @@ end
 
 # test fixed forcing (precipitation = 2.5)
 config = Wflow.Config(tomlpath)
-config.input.vertical.atmospheric_forcing.precipitation = Dict("value" => 2.5)
+config.input.forcing.atmosphere_water__precipitation_volume_flux = Dict("value" => 2.5)
 model = Wflow.initialize_sbm_model(config)
 Wflow.load_fixed_forcing!(model)
 
@@ -288,7 +286,7 @@ Wflow.run_timestep!(model)
 
 @testset "river and overland flow and depth (local inertial)" begin
     q = model.lateral.river.variables.q_av
-    @test sum(q) ≈ 2380.64389229669f0
+    @test sum(q) ≈ 2495.9830572223946f0
     @test q[1622] ≈ 7.30561606758937f-5
     @test q[43] ≈ 5.3566292152594155f0
     @test q[501] ≈ 1.602564408503896f0
@@ -314,9 +312,9 @@ config = Wflow.Config(tomlpath)
 config.model.floodplain_1d = true
 config.model.river_routing = "local-inertial"
 config.model.land_routing = "kinematic-wave"
-config.input["floodplain_water__sum_of_volume-per-depth"] = "floodplain_volume"
-Dict(config.state.lateral.river)["floodplain.variables"] =
-    Dict("q" => "q_floodplain", "h" => "h_floodplain")
+config.input.parameters["floodplain_water__sum_of_volume-per-depth"] = "floodplain_volume"
+config.state.variables.floodplain_water__volume = "q_floodplain"
+config.state.variables.floodplain_water__depth = "h_floodplain"
 
 model = Wflow.initialize_sbm_model(config)
 
@@ -432,8 +430,8 @@ Wflow.run_timestep!(model)
 end
 
 # set boundary condition local inertial routing from netCDF file
-config.input["model_boundary_condition~river__length"] = "riverlength_bc"
-config.input["model_boundary_condition~river_bank_water__depth"] = "riverdepth_bc"
+config.input.parameters["model_boundary_condition~river__length"] = "riverlength_bc"
+config.input.parameters["model_boundary_condition~river_bank_water__depth"] = "riverdepth_bc"
 model = Wflow.initialize_sbm_model(config)
 Wflow.run_timestep!(model)
 Wflow.run_timestep!(model)
@@ -458,9 +456,11 @@ Wflow.close_files(model; delete_output = false)
     i = 100
     tomlpath = joinpath(@__DIR__, "sbm_config.toml")
     config = Wflow.Config(tomlpath)
-    config.input.soil_water__vertical_saturated_hydraulic_conductivity = "kv"
-    config.input["soil_vertical_saturated_hydraulic_conductivity_profile~exponential_below-surface__depth"] = Dict("value" => 400.0)
-    config.input["soil_vertical_saturated_hydraulic_conductivity_profile~layered_below-surface__depth"] = Dict("value" => 400.0)
+    config.input.parameters.soil_water__vertical_saturated_hydraulic_conductivity = "kv"
+    config.input.parameters["soil_vertical_saturated_hydraulic_conductivity_profile~exponential_below-surface__depth"] =
+        Dict("value" => 400.0)
+    config.input.parameters["soil_vertical_saturated_hydraulic_conductivity_profile~layered_below-surface__depth"] =
+        Dict("value" => 400.0)
 
     @testset "exponential profile" begin
         model = Wflow.initialize_sbm_model(config)
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 1ca02562e..38b37877e 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -71,7 +71,10 @@ end
 
 @testset "no drains" begin
     config.model.drains = false
-    delete!(Dict(config.output.lateral.subsurface), "drain")
+    delete!(
+        Dict(config.output.variables),
+        "land_drain_water~to-subsurface__volume_flow_rate",
+    )
     model = Wflow.initialize_sbm_gwf_model(config)
     @test collect(keys(model.lateral.subsurface)) == [:flow, :recharge, :river]
 end
@@ -83,8 +86,8 @@ tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
 config = Wflow.Config(tomlpath)
 config.model.river_routing = "local-inertial"
 
-config.input.river_bank_water__elevation = "bankfull_elevation"
-config.input.river_bank_water__depth = "bankfull_depth"
+config.input.parameters.river_bank_water__elevation = "bankfull_elevation"
+config.input.parameters.river_bank_water__depth = "bankfull_depth"
 
 model = Wflow.initialize_sbm_gwf_model(config)
 Wflow.run_timestep!(model)
@@ -108,14 +111,14 @@ config = Wflow.Config(tomlpath)
 config.model.river_routing = "local-inertial"
 config.model.land_routing = "local-inertial"
 
-config.input.river_bank_water__elevation = "bankfull_elevation"
-config.input.river_bank_water__depth = "bankfull_depth"
-config.input.land_surface_water_flow__ground_elevation = "wflow_dem"
+config.input.parameters.river_bank_water__elevation = "bankfull_elevation"
+config.input.parameters.river_bank_water__depth = "bankfull_depth"
+config.input.parameters.land_surface_water_flow__ground_elevation = "wflow_dem"
 
-pop!(Dict(config.state.lateral.land.variables), "q")
-config.state.lateral.land.variables.h_av = "h_av_land"
-config.state.lateral.land.variables.qx = "qx_land"
-config.state.lateral.land.variables.qy = "qy_land"
+pop!(Dict(config.state.variables), "land_surface_water__volume_flow_rate")
+config.state.variables.land_surface_water__time_average_of_depth = "h_av_land"
+config.state.variables.land_surface_water__x_component_of_volume_flow_rate = "qx_land"
+config.state.variables.land_surface_water__y_component_of_volume_flow_rate = "qy_land"
 
 model = Wflow.initialize_sbm_gwf_model(config)
 Wflow.run_timestep!(model)
diff --git a/test/runtests.jl b/test/runtests.jl
index bc6d2764d..be449af27 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,5 @@
 ## load test dependencies and set paths to testing data
+using Accessors
 using Dates
 using Downloads
 using Graphs
@@ -90,9 +91,8 @@ with_logger(NullLogger()) do
         include("groundwater.jl")
         include("utils.jl")
         include("bmi.jl")
-        include("run_sediment.jl")
+        #include("run_sediment.jl")
         include("subdomains.jl")
-
         Aqua.test_all(Wflow; ambiguities = false, persistent_tasks = false)
     end
 end
diff --git a/test/sbm_config.toml b/test/sbm_config.toml
index 2acfec3fa..45f9e0b7d 100644
--- a/test/sbm_config.toml
+++ b/test/sbm_config.toml
@@ -19,33 +19,27 @@ path_output = "outstates-moselle.nc"
 # if listed, the variable must be present in the NetCDF or error
 # if not listed, the variable can get a default value if it has one
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
 
-[state.lateral.river.variables]
-h = "h_river"
-h_av = "h_av_river"
-q = "q_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
+reservoir_water__volume = "volume_reservoir"
 
-[state.lateral.subsurface.variables]
-ssf = "ssf"
+subsurface_water__volume_flow_rate  = "ssf"
 
-[state.lateral.land.variables]
-h = "h_land"
-h_av = "h_av_land"
-q = "q_land"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
 
 [input]
 path_forcing = "forcing-moselle.nc"
@@ -57,6 +51,12 @@ ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -95,8 +95,8 @@ river_bank_water__elevation = "RiverZ"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
-reservoir_areas__number = "wflow_reservoirareas"
-reservoir_locations__number = "wflow_reservoirlocs"
+reservoir_area__number = "wflow_reservoirareas"
+reservoir_location__number = "wflow_reservoirlocs"
 reservoir_surface__area = "ResSimpleArea"
 "reservoir_water_demand~required~downstream__volume_flow_rate" = "ResDemand"
 reservoir_water_release-below-spillway__max_volume_flow_rate = "ResMaxRelease"
@@ -108,20 +108,16 @@ subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
-[input.soil_surface_water__vertical_saturated_hydraulic_conductivity]
+[input.parameters.soil_surface_water__vertical_saturated_hydraulic_conductivity]
 netcdf.variable.name = "KsatVer"
 scale = 1.0
 offset = 0.0
@@ -140,31 +136,19 @@ min_streamorder_land = 5
 [output]
 path = "output_moselle.nc"
 
-[output.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[output.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[output.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[output.lateral.river.variables]
-h = "h_river"
-q = "q_river"
-
-[output.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
-
-[output.lateral.subsurface.variables]
-ssf = "ssf"
-
-[output.lateral.land.variables]
-h = "h_land"
-q = "q_land"
+[output.variables]
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+river_water__depth = "h_river"
+river_water__volume_flow_rate = "q_river"
+reservoir_water__volume = "volume_reservoir"
+subsurface_water__volume_flow_rate  = "ssf"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+vertical.interception.variables.canopy_storage = "canopystorage"
 
 [netcdf]
 path = "output_scalar_moselle.nc"
@@ -179,33 +163,33 @@ coordinate.x = 6.255
 coordinate.y = 50.012
 name = "temp_coord"
 location = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[netcdf.variable]]
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [csv]
 path = "output_moselle.csv"
 
 [[csv.column]]
 header = "Q"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 reducer = "maximum"
 
 [[csv.column]]
 header = "volume"
 index = 1
-parameter = "lateral.river.boundary_conditions.reservoir.variables.volume"
+parameter = "reservoir_water__volume"
 
 [[csv.column]]
 coordinate.x = 6.255
 coordinate.y = 50.012
 header = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 coordinate.x = 6.255
@@ -218,17 +202,17 @@ layer = 2
 header = "temp_byindex"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 
 [[csv.column]]
 header = "recharge"
 map = "subcatchment"
-parameter = "vertical.soil.variables.recharge"
+parameter = "soil_water_sat-zone_top__recharge_volume_flux"
 reducer = "mean"
 
 [API]
diff --git a/test/sbm_gw.toml b/test/sbm_gw.toml
index 92004aad9..96a17df5a 100644
--- a/test/sbm_gw.toml
+++ b/test/sbm_gw.toml
@@ -18,30 +18,27 @@ path_output = "outstates-moselle.nc"
 # if listed, the variable must be present in the NetCDF or error
 # if not listed, the variable can get a default value if it has one
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.soil.variables] 
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[state.vertical.snow.variables] 
-snow_storage = "snow"
-snow_water = "snowwater"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
 
-[state.lateral.river.variables]
-h = "h_river"
-h_av = "h_av_river"
-q = "q_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
+reservoir_water__volume = "volume_reservoir"
 
-[state.lateral.land.variables]
-h = "h_land"
-h_av = "h_av_land"
-q = "q_land"
+subsurface_water__volume_flow_rate  = "ssf"
+
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
 
 [input]
 path_forcing = "forcing-moselle.nc"
@@ -53,6 +50,12 @@ ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -88,8 +91,8 @@ river__width = "wflow_riverwidth"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
-reservoir_areas__number = "wflow_reservoirareas"
-reservoir_locations__number = "wflow_reservoirlocs"
+reservoir_area__number = "wflow_reservoirareas"
+reservoir_location__number = "wflow_reservoirlocs"
 reservoir_surface__area = "ResSimpleArea"
 "reservoir_water_demand~required~downstream__volume_flow_rate" = "ResDemand"
 reservoir_water_release-below-spillway__max_volume_flow_rate = "ResMaxRelease"
@@ -99,19 +102,15 @@ reservoir_water__max_volume = "ResMaxVolume"
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
 [model]
 kin_wave_iteration = true
 masswasting = true
@@ -125,28 +124,18 @@ type = "sbm"
 [output]
 path = "output_moselle.nc"
 
-[output.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[output.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[output.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[output.lateral.river.variables]
-h = "h_river"
-q = "q_river"
-
-[output.lateral.river.boundary_conditions.reservoir.variables]
-volume = "volume_reservoir"
-
-[output.lateral.land.variables]
-h = "h_land"
-q = "q_land"
+[output.variables]
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+river_water__depth = "h_river"
+river_water__volume_flow_rate = "q_river"
+reservoir_water__volume = "volume_reservoir"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+vertical.interception.variables.canopy_storage = "canopystorage"
 
 [netcdf]
 path = "output_scalar_moselle.nc"
@@ -154,56 +143,56 @@ path = "output_scalar_moselle.nc"
 [[netcdf.variable]]
 name = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 
 [[netcdf.variable]]
 coordinate.x = 6.255
 coordinate.y = 50.012
 name = "temp_coord"
 location = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[netcdf.variable]]
 location = "temp_byindex"
 name = "temp_index"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [csv]
 path = "output_moselle.csv"
 
 [[csv.column]]
 header = "Q"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 reducer = "maximum"
 
 [[csv.column]]
 header = "volume"
 index = 1
-parameter = "lateral.river.boundary_conditions.reservoir.variables.volume"
+parameter = "reservoir_water__volume"
 
 [[csv.column]]
 coordinate.x = 6.255
 coordinate.y = 50.012
 header = "temp_bycoord"
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 header = "temp_byindex"
 index.x = 100
 index.y = 264
-parameter = "vertical.atmospheric_forcing.temperature"
+parameter = "atmosphere_air__temperature"
 
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 
 [[csv.column]]
 header = "recharge"
 map = "subcatchment"
-parameter = "vertical.soil.variables.recharge"
+parameter = "soil_water_sat-zone_top__recharge_volume_flux"
 reducer = "mean"
 
 [API]
diff --git a/test/sbm_gwf_config.toml b/test/sbm_gwf_config.toml
index 29d28ec3d..ac21c4e49 100644
--- a/test/sbm_gwf_config.toml
+++ b/test/sbm_gwf_config.toml
@@ -18,25 +18,21 @@ path_output = "outstates-example-sbm-gwf.nc"
 # if listed, the variable must be present in the NetCDF or error
 # if not listed, the variable can get a default value if it has one
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-ustorelayerdepth = "ustorelayerdepth"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[state.lateral.river.variables]
-h = "h_river"
-h_av = "h_av_river"
-q = "q_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.land.variables]
-h = "h_land"
-h_av = "h_av_land"
-q = "q_land"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
 
-[state.lateral.subsurface.flow.aquifer.variables]
-head = "head"
+subsurface_water__hydraulic_head = "head"
 
 [input]
 path_forcing = "forcing-sbm-groundwater-part*.nc"
@@ -48,6 +44,11 @@ river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 altitude = "wflow_dem"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "P"
+land_surface_water__potential_evaporation_volume_flux = "PET"
+
+[input.parameters]
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "kv"
 "soil~non-compacted_surface_water__infiltration_capacity" = "InfiltCapSoil"
 soil_water__residual_volume_fraction = "thetaR"
@@ -68,28 +69,24 @@ river__width = "wflow_riverwidth"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
+subsurface_surface_water__horizontal_saturated_hydraulic_conductivity = "k"
+"model_boundary_condition~constant_hydraulic_head" = "constant_head"
+land_drain_location__flag = "drain"
+land_drain__conductance = "cond_drain"
+land_drain__elevation = "elev_drain"
+subsurface_water__specific_yield = "specific_yield"
+river_bottom__elevation = "river_bottom"
+river_water__infiltration_conductance = "infiltration_conductance"
+river_water__exfiltration_conductance = "infiltration_conductance"
+
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "PET"
-precipitation = "P"
-
-[input.lateral.subsurface]
-conductivity = "k"
-constant_head = "constant_head"
-drain = "drain"
-drain_conductance = "cond_drain"
-drain_elevation = "elev_drain"
-exfiltration_conductance = "exfiltration_conductance"
-infiltration_conductance = "infiltration_conductance"
-river_bottom = "river_bottom"
-specific_yield = "specific_yield"
-
 [model]
 constanthead = true
 drains = true
@@ -103,30 +100,16 @@ type = "sbm_gwf"
 [output]
 path = "output_example-sbm-gwf.nc"
 
-[output.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[output.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-ustorelayerdepth = "ustorelayerdepth"
-
-[output.vertical.soil.parameters]
-soilthickness = "soilthickness"
-
-[output.lateral.river.variables]
-q = "q"
-
-[output.lateral.subsurface.flow.aquifer.variables]
-head = "head"
-
-[output.lateral.subsurface.recharge.variables]
-rate = "rate"
-
-[output.lateral.subsurface.drain.variables]
-flux = "drain_flux"
-
-[output.lateral.subsurface.river.variables]
-flux = "flux"
+[output.variables]
+vegetation_canopy_water__storage = "canopystorage"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+vertical.soil.parameters.soilthickness = "soilthickness"
+river_water__volume_flow_rate = "q"
+subsurface_water__hydraulic_head = "head"
+subsurface_water_sat-zone_top__net-recharge_volume_flux = "rate"
+"land_drain_water~to-subsurface__volume_flow_rate" = "drain_flux"
+"river_water~to-subsurface__volume_flow_rate" = "flux"
 
 [csv]
 path = "output_example-sbm-gwf.csv"
@@ -134,9 +117,9 @@ path = "output_example-sbm-gwf.csv"
 [[csv.column]]
 header = "Q_av"
 index = 5
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"
 
 [[csv.column]]
 header = "head"
 index = 5
-parameter = "lateral.subsurface.flow.aquifer.variables.head"
+parameter = "subsurface_water__hydraulic_head"
diff --git a/test/sbm_gwf_piave_demand_config.toml b/test/sbm_gwf_piave_demand_config.toml
index 78a505eff..476921300 100644
--- a/test/sbm_gwf_piave_demand_config.toml
+++ b/test/sbm_gwf_piave_demand_config.toml
@@ -11,9 +11,33 @@ loglevel = "info"
 path_input = "instates-piave-gwf.nc"
 path_output = "outstates-piave-gwf.nc"
 
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
+
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
+
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
+"land_surface_water~paddy__depth" = "h_paddy"
+
+subsurface_water__hydraulic_head = "head"
+
+glacier_ice__leq-volume = "glacierstore"
+
 [input]
 path_forcing = "forcing-piave.nc"
 path_static = "staticmaps-piave.nc"
+
 ldd = "wflow_ldd"
 river_location = "wflow_river"
 altitude = "wflow_dem"
@@ -21,6 +45,12 @@ subcatchment = "wflow_subcatch"
 gauges = "wflow_gauges"
 gauges_grdc = "wflow_gauges_grdc"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -34,6 +64,7 @@ glacier_surface__area_fraction = "wflow_glacierfrac"
 glacier_ice__degree-day_coefficient = "G_Cfmax"
 glacier_ice__melting_temperature_threshold = "G_TT"
 "glacier_firn_accumulation__snowpack~dry_leq-depth_fraction"  = "G_SIfrac"
+glacier_ice__leq-volume = "wflow_glacierstore"
 
 soil_water__brooks-corey_epsilon_parameter = "c"
 soil_surface_water__infiltration_reduction_parameter = "cf_soil"
@@ -72,21 +103,44 @@ river_bank_water__depth = "RiverDepth"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
+land_water_allocation_area__number = "allocation_areas"
+"land~domestic__gross_water_demand_flux" = "dom_gross"
+"land~domestic__net_water_demand_flux" = "dom_net"
+"land~industry__gross_water_demand_flux" = "ind_gross"
+"land~industry__net_water_demand_flux" = "ind_net"
+"land~livestock__gross_water_demand_flux" = "lsk_gross"
+"land~livestock__net_water_demand_flux" = "lsk_net"
+land_surface_water__withdrawal_fraction = "SurfaceWaterFrac"
+"land~irrigated-paddy_area__number" = "paddy_irrigation_areas"
+"land~irrigated-non-paddy_area__number" = "nonpaddy_irrigation_areas"
+"land~irrigated-paddy__irrigation_trigger_flag" = "irrigation_trigger"
+"land~irrigated-non-paddy__irrigation_trigger_flag" = "irrigation_trigger"
+
+"model_boundary_condition~constant_hydraulic_head" = "constant_head"
+subsurface_surface_water__horizontal_saturated_hydraulic_conductivity = "kh_surface"
+river_bottom__elevation = "zb_river"
+river_water__infiltration_conductance = "riverbed_cond"
+river_water__exfiltration_conductance = "riverbed_cond"
+subsurface_water__specific_yield = "specific_yield"
+subsurface__horizontal_saturated_hydraulic_conductivity_scale_parameter = "gwf_f"
+
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
+
 cyclic = [ 
     "vegetation__leaf-area_index", 
-    "vertical.demand.domestic.demand.demand_gross", 
-    "vertical.demand.domestic.demand.demand_net", 
-    "vertical.demand.industry.demand.demand_gross", 
-    "vertical.demand.industry.demand.demand_net", 
-    "vertical.demand.livestock.demand.demand_gross", 
-    "vertical.demand.livestock.demand.demand_net", 
-    "vertical.demand.paddy.parameters.irrigation_trigger", 
-    "vertical.demand.nonpaddy.parameters.irrigation_trigger",
+    "land~domestic__gross_water_demand_flux",
+    "land~domestic__net_water_demand_flux", 
+    "land~industry__gross_water_demand_flux", 
+    "land~industry__net_water_demand_flux", 
+    "land~livestock__gross_water_demand_flux", 
+    "land~livestock__net_water_demand_flux", 
+    "land~irrigated-paddy__irrigation_trigger_flag", 
+    "land~irrigated-non-paddy__irrigation_trigger_flag",
 ]
 
 [model]
@@ -103,35 +157,7 @@ kw_river_tstep = 900
 kw_land_tstep = 3600
 thicknesslayers = [ 50, 100, 50, 200, 800,]
 river_routing = "kinematic-wave"
-
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[state.vertical.glacier.variables]
-glacier_store = "glacierstore"
-
-[state.vertical.demand.paddy.variables]
-h = "h_paddy"
-
-[state.lateral.subsurface.flow.aquifer.variables]
-head = "head"
-
-[input.vertical.glacier.variables]
-glacier_store = "wflow_glacierstore"
-
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
+conductivity_profile = "exponential"
 
 [model.water_demand]
 domestic = true
@@ -140,64 +166,13 @@ livestock = true
 paddy = true
 nonpaddy = true
 
-[state.lateral.river.variables]
-q = "q_river"
-h = "h_river"
-h_av = "h_av_river"
-
-[state.lateral.subsurface.variables]
-ssf = "ssf"
-
-[state.lateral.land.variables]
-q = "q_land"
-h = "h_land"
-h_av = "h_av_land"
-
-[input.vertical.allocation.parameters]
-areas = "allocation_areas"
-frac_sw_used = "SurfaceWaterFrac"
-
-[input.vertical.demand.domestic.demand]
-demand_gross = "dom_gross"
-demand_net = "dom_net"
-
-[input.vertical.demand.industry.demand]
-demand_gross = "ind_gross"
-demand_net = "ind_net"
-
-[input.vertical.demand.livestock.demand]
-demand_gross = "lsk_gross"
-demand_net = "lsk_net"
-
-[input.vertical.demand.paddy.parameters]
-irrigation_areas = "paddy_irrigation_areas"
-irrigation_trigger = "irrigation_trigger"
-
-[input.vertical.demand.nonpaddy.parameters]
-irrigation_areas = "nonpaddy_irrigation_areas"
-irrigation_trigger = "irrigation_trigger"
-
-[input.lateral.subsurface]
-constant_head = "constant_head"
-conductivity_profile = "exponential"
-conductivity = "kh_surface"
-exfiltration_conductance = "riverbed_cond"
-infiltration_conductance = "riverbed_cond"
-river_bottom = "zb_river"
-specific_yield = "specific_yield"
-gwf_f = "gwf_f"
-
 [output]
 path = "output-piave-gwf.nc"
 
-[output.lateral.river.variables]
-q_av = "q_river"
-
-[output.vertical.soil.variables]
-zi = "zi"
-
-[output.lateral.subsurface.flow.aquifer.variables]
-head = "head"
+[output.variables]
+river_water__time_average_of_volume_flow_rate = "q_river"
+subsurface_water__hydraulic_head = "head"
+soil_water_sat-zone_top__depth = "zi"
 
 [csv]
 path = "output.csv"
@@ -205,9 +180,9 @@ path = "output.csv"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"
 
 [[csv.column]]
 header = "Q"
 map = "gauges_grdc"
-parameter = "lateral.river.variables.q_av"
\ No newline at end of file
+parameter = "river_water__time_average_of_volume_flow_rate"
\ No newline at end of file
diff --git a/test/sbm_piave_config.toml b/test/sbm_piave_config.toml
index 1dd6a77fe..e19d36771 100644
--- a/test/sbm_piave_config.toml
+++ b/test/sbm_piave_config.toml
@@ -14,12 +14,19 @@ path_output = "outstates-piave.nc"
 [input]
 path_forcing = "forcing-piave.nc"
 path_static = "staticmaps-piave.nc"
+
 ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 gauges = "wflow_gauges"
 gauges_grdc = "wflow_gauges_grdc"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -32,6 +39,7 @@ glacier_surface__area_fraction = "wflow_glacierfrac"
 glacier_ice__degree-day_coefficient = "G_Cfmax"
 glacier_ice__melting_temperature_threshold = "G_TT"
 "glacier_firn_accumulation__snowpack~dry_leq-depth_fraction"  = "G_SIfrac"
+glacier_ice__leq-volume = "wflow_glacierstore"
 
 soil_water__brooks-corey_epsilon_parameter = "c"
 soil_surface_water__infiltration_reduction_parameter = "cf_soil"
@@ -71,10 +79,11 @@ land_surface__slope = "Slope"
 
 subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio = "KsatHorFrac"
 
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
@@ -93,47 +102,33 @@ kw_land_tstep = 3600
 thicknesslayers = [ 50, 100, 50, 200, 800,]
 river_routing = "kinematic-wave"
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.glacier.variables]
-glacier_store = "glacierstore"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[input.vertical.glacier.variables]
-glacier_store = "wflow_glacierstore"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
+glacier_ice__leq-volume = "glacierstore"
 
-[state.lateral.river.variables]
-q = "q_river"
-h = "h_river"
-h_av = "h_av_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.subsurface.variables]
-ssf = "ssf"
+subsurface_water__volume_flow_rate  = "ssf"
 
-[state.lateral.land.variables]
-q = "q_land"
-h = "h_land"
-h_av = "h_av_land"
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
 
 [output]
 path = "output-piave.nc"
 
-[output.lateral.river.variables]
-q_av = "q_river"
+[output.variables]
+river_water__time_average_of_volume_flow_rate = "q_river"
 
 [csv]
 path = "output-piave.csv"
@@ -141,9 +136,9 @@ path = "output-piave.csv"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"
 
 [[csv.column]]
 header = "Q"
 map = "gauges_grdc"
-parameter = "lateral.river.variables.q_av"
\ No newline at end of file
+parameter = "river_water__time_average_of_volume_flow_rate"
\ No newline at end of file
diff --git a/test/sbm_piave_demand_config.toml b/test/sbm_piave_demand_config.toml
index d359b9b77..e2558ad86 100644
--- a/test/sbm_piave_demand_config.toml
+++ b/test/sbm_piave_demand_config.toml
@@ -11,6 +11,29 @@ loglevel = "info"
 path_input = "instates-piave.nc"
 path_output = "outstates-piave.nc"
 
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
+
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
+
+land_surface_water__volume_flow_rate = "q_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
+"land_surface_water~paddy__depth" = "h_paddy"
+
+subsurface_water__volume_flow_rate  = "ssf"
+
+glacier_ice__leq-volume = "glacierstore"
+
 [input]
 path_forcing = "forcing-piave.nc"
 path_static = "staticmaps-piave.nc"
@@ -20,6 +43,12 @@ subcatchment = "wflow_subcatch"
 gauges = "wflow_gauges"
 gauges_grdc = "wflow_gauges_grdc"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -32,6 +61,7 @@ glacier_surface__area_fraction = "wflow_glacierfrac"
 glacier_ice__degree-day_coefficient = "G_Cfmax"
 glacier_ice__melting_temperature_threshold = "G_TT"
 "glacier_firn_accumulation__snowpack~dry_leq-depth_fraction" = "G_SIfrac"
+glacier_ice__leq-volume = "wflow_glacierstore"
 
 soil_water__brooks-corey_epsilon_parameter = "c"
 soil_surface_water__infiltration_reduction_parameter = "cf_soil"
@@ -72,21 +102,36 @@ land_surface__slope = "Slope"
 
 subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio = "KsatHorFrac"
 
+land_water_allocation_area__number = "allocation_areas"
+"land~domestic__gross_water_demand_flux" = "dom_gross"
+"land~domestic__net_water_demand_flux" = "dom_net"
+"land~industry__gross_water_demand_flux" = "ind_gross"
+"land~industry__net_water_demand_flux" = "ind_net"
+"land~livestock__gross_water_demand_flux" = "lsk_gross"
+"land~livestock__net_water_demand_flux" = "lsk_net"
+land_surface_water__withdrawal_fraction = "SurfaceWaterFrac"
+"land~irrigated-paddy_area__number" = "paddy_irrigation_areas"
+"land~irrigated-non-paddy_area__number" = "nonpaddy_irrigation_areas"
+"land~irrigated-paddy__irrigation_trigger_flag" = "irrigation_trigger"
+"land~irrigated-non-paddy__irrigation_trigger_flag" = "irrigation_trigger"
+
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
+
 cyclic = [ 
     "vegetation__leaf-area_index",
-    "vertical.demand.domestic.demand.demand_gross", 
-    "vertical.demand.domestic.demand.demand_net", 
-    "vertical.demand.industry.demand.demand_gross", 
-    "vertical.demand.industry.demand.demand_net", 
-    "vertical.demand.livestock.demand.demand_gross", 
-    "vertical.demand.livestock.demand.demand_net", 
-    "vertical.demand.paddy.parameters.irrigation_trigger", 
-    "vertical.demand.nonpaddy.parameters.irrigation_trigger",
+    "land~domestic__gross_water_demand_flux",
+    "land~domestic__net_water_demand_flux", 
+    "land~industry__gross_water_demand_flux", 
+    "land~industry__net_water_demand_flux", 
+    "land~livestock__gross_water_demand_flux", 
+    "land~livestock__net_water_demand_flux", 
+    "land~irrigated-paddy__irrigation_trigger_flag", 
+    "land~irrigated-non-paddy__irrigation_trigger_flag",
 ]
 
 [model]
@@ -103,32 +148,6 @@ kw_land_tstep = 3600
 thicknesslayers = [ 50, 100, 50, 200, 800,]
 river_routing = "kinematic-wave"
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[state.vertical.glacier.variables]
-glacier_store = "glacierstore"
-
-[state.vertical.demand.paddy.variables]
-h = "h_paddy"
-
-[input.vertical.glacier.variables]
-glacier_store = "wflow_glacierstore"
-
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
 [model.water_demand]
 domestic = true
 industry = true
@@ -136,51 +155,12 @@ livestock = true
 paddy = true
 nonpaddy = true
 
-[state.lateral.river.variables]
-q = "q_river"
-h = "h_river"
-h_av = "h_av_river"
-
-[state.lateral.subsurface.variables]
-ssf = "ssf"
-
-[state.lateral.land.variables]
-q = "q_land"
-h = "h_land"
-h_av = "h_av_land"
-
-[input.vertical.allocation.parameters]
-areas = "allocation_areas"
-frac_sw_used = "SurfaceWaterFrac"
-
-[input.vertical.demand.domestic.demand]
-demand_gross = "dom_gross"
-demand_net = "dom_net"
-
-[input.vertical.demand.industry.demand]
-demand_gross = "ind_gross"
-demand_net = "ind_net"
-
-[input.vertical.demand.livestock.demand]
-demand_gross = "lsk_gross"
-demand_net = "lsk_net"
-
-[input.vertical.demand.paddy.parameters]
-irrigation_areas = "paddy_irrigation_areas"
-irrigation_trigger = "irrigation_trigger"
-
-[input.vertical.demand.nonpaddy.parameters]
-irrigation_areas = "nonpaddy_irrigation_areas"
-irrigation_trigger = "irrigation_trigger"
-
 [output]
 path = "output-piave-demand.nc"
 
-[output.lateral.river.variables]
-q_av = "q_river"
-
-[output.vertical.soil.variables]
-zi = "zi"
+[output.variables]
+river_water__time_average_of_volume_flow_rate = "q_river"
+soil_water_sat-zone_top__depth = "zi"
 
 [csv]
 path = "output-piave-demand.csv"
@@ -188,21 +168,21 @@ path = "output-piave-demand.csv"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"
 
 [[csv.column]]
 header = "Q"
 map = "gauges_grdc"
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"
 
 [[csv.column]]
 coordinate.x = 12.7243
 coordinate.y = 45.5851
 header = "paddy_h_bycoord"
-parameter = "vertical.demand.paddy.variables.h"
+parameter = "land_surface_water~paddy__depth"
 
 [[csv.column]]
 coordinate.x = 12.7243
 coordinate.y = 45.5851
 header = "irri_bycoord"
-parameter = "vertical.allocation.variables.irri_alloc"
\ No newline at end of file
+parameter = "land~irrigated__allocated_water_volume_flux"
\ No newline at end of file
diff --git a/test/sbm_simple.toml b/test/sbm_simple.toml
index e4b796450..06552aaec 100644
--- a/test/sbm_simple.toml
+++ b/test/sbm_simple.toml
@@ -17,6 +17,12 @@ ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -53,8 +59,8 @@ river__width = "wflow_riverwidth"
 land_surface_water_flow__manning_n_parameter = "N"
 land_surface__slope = "Slope"
 
-reservoir_areas__number = "wflow_reservoirareas"
-reservoir_locations__number = "wflow_reservoirlocs"
+reservoir_area__number = "wflow_reservoirareas"
+reservoir_location__number = "wflow_reservoirlocs"
 reservoir_surface__area = "ResSimpleArea"
 "reservoir_water_demand~required~downstream__volume_flow_rate" = "ResDemand"
 reservoir_water_release-below-spillway__max_volume_flow_rate = "ResMaxRelease"
@@ -66,26 +72,15 @@ subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
 
-[input.vertical.vegetation_parameter_set]
-leaf_area_index = "LAI"
-kext = "Kext"
-storage_specific_leaf = "Sl"
-storage_wood = "Swood"
-rootingdepth = "RootingDepth"
-
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
 [model]
 thicknesslayers = [100, 300, 800]
 type = "sbm"
@@ -97,9 +92,9 @@ path = "output_moselle_simple.csv"
 coordinate.x = 7.378
 coordinate.y = 50.204
 header = "Q"
-parameter = "lateral.river.variables.q"
+parameter = "river_water__volume_flow_rate"
 
 [[csv.column]]
 header = "recharge"
-parameter = "vertical.soil.variables.recharge"
+parameter = "soil_water_sat-zone_top__recharge_volume_flux"
 reducer = "mean"
diff --git a/test/sbm_swf_config.toml b/test/sbm_swf_config.toml
index 7e479ea0d..8cf912277 100644
--- a/test/sbm_swf_config.toml
+++ b/test/sbm_swf_config.toml
@@ -14,34 +14,28 @@ dir_output = "data/output"
 # if listed, the variable must be present in the NetCDF or error
 # if not listed, the variable can get a default value if it has one
 
-[state.vertical.interception.variables]
-canopy_storage = "canopystorage"
+[state.variables]
+vegetation_canopy_water__storage = "canopystorage"
 
-[state.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
 
-[state.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
 
-[state.lateral.river.variables]
-h = "h_river"
-h_av = "h_av_river"
-q = "q_river"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "h_av_river"
+river_water__volume_flow_rate = "q_river"
 
-[state.lateral.river.reservoir]
-volume = "volume_reservoir"
+reservoir_water__volume = "volume_reservoir"
 
-[state.lateral.subsurface.variables]
-ssf = "ssf"
+subsurface_water__volume_flow_rate  = "ssf"
 
-[state.lateral.land.variables]
-h = "h_land"
-h_av = "h_av_land"
-qx = "qx_land"
-qy = "qy_land"
+land_surface_water__depth = "h_land"
+land_surface_water__time_average_of_depth = "h_av_land"
+land_surface_water__x_component_of_volume_flow_rate = "qx_land"
+land_surface_water__y_component_of_volume_flow_rate = "qy_land"
 
 [input]
 path_forcing = "forcing-moselle.nc"
@@ -53,6 +47,12 @@ ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
 
+[input.forcing]
+atmosphere_water__precipitation_volume_flux = "precip"
+land_surface_water__potential_evaporation_volume_flux = "pet"
+atmosphere_air__temperature = "temp"
+
+[input.parameters]
 atmosphere_air__snowfall_temperature_threshold = "TT"
 atmosphere_air__snowfall_temperature_interval = "TTI"
 
@@ -92,8 +92,8 @@ land_surface_water_flow__manning_n_parameter = "N"
 land_surface_water_flow__ground_elevation = "FloodplainZ"
 land_surface__slope = "Slope"
 
-reservoir_areas__number = "wflow_reservoirareas"
-reservoir_locations__number = "wflow_reservoirlocs"
+reservoir_area__number = "wflow_reservoirareas"
+reservoir_location__number = "wflow_reservoirlocs"
 reservoir_surface__area = "ResSimpleArea"
 "reservoir_water_demand~required~downstream__volume_flow_rate" = "ResDemand"
 reservoir_water_release-below-spillway__max_volume_flow_rate = "ResMaxRelease"
@@ -105,23 +105,20 @@ subsurface_water__horizontal-to-vertical_saturated_hydraulic_conductivity_ratio
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
+[input.dynamic]
 forcing = [
-  "vertical.atmospheric_forcing.precipitation",
-  "vertical.atmospheric_forcing.temperature",
-  "vertical.atmospheric_forcing.potential_evaporation",
+  "atmosphere_water__precipitation_volume_flux",
+  "atmosphere_air__temperature",
+  "land_surface_water__potential_evaporation_volume_flux",
 ]
 
 cyclic = ["vegetation__leaf-area_index"]
 
-[input.vertical.atmospheric_forcing]
-potential_evaporation = "pet"
-precipitation = "precip"
-temperature = "temp"
-
 [model]
 kin_wave_iteration = true
 masswasting = true
 reinit = true
+reservoirs = true
 snow = true
 thicknesslayers = [100, 300, 800]
 min_streamorder = 3
@@ -132,31 +129,22 @@ type = "sbm"
 [output]
 path = "output_moselle_swf.nc"
 
-[output.vertical.interception.variables]
-canopy_storage = "canopystorage"
-
-[output.vertical.soil.variables]
-satwaterdepth = "satwaterdepth"
-tsoil = "tsoil"
-ustorelayerdepth = "ustorelayerdepth"
-
-[output.vertical.snow.variables]
-snow_storage = "snow"
-snow_water = "snowwater"
-
-[output.lateral.river.variables]
-h = "h_river"
-h_av = "hav_river"
-q = "q_river"
-q_av = "qav_river"
-
-[output.lateral.subsurface.variables]
-ssf = "ssf"
-
-[output.lateral.land.variables]
-h = "h_land"
-qx = "qx_land"
-qy = "qy_land"
+[output.variables]
+soil_water_sat-zone__depth = "satwaterdepth"
+soil_surface__temperature = "tsoil"
+soil_water_unsat-zone__depth-per-soil_layer = "ustorelayerdepth"
+"snowpack~dry__leq-depth" = "snow"
+"snowpack~liquid__depth" = "snowwater"
+river_water__depth = "h_river"
+river_water__time_average_of_depth = "hav_river"
+river_water__volume_flow_rate = "q_river"
+river_water__time_average_of_volume_flow_rate = "qav_river"
+reservoir_water__volume = "volume_reservoir"
+subsurface_water__volume_flow_rate  = "ssf"
+land_surface_water__depth = "h_land"
+land_surface_water__x_component_of_volume_flow_rate = "qx_land"
+land_surface_water__y_component_of_volume_flow_rate = "qy_land"
+vertical.interception.variables.canopy_storage = "canopystorage"
 
 [csv]
 path = "output_moselle_swf.csv"
@@ -164,5 +152,5 @@ path = "output_moselle_swf.csv"
 [[csv.column]]
 header = "Q"
 map = "gauges"
-parameter = "lateral.river.variables.q_av"
+parameter = "river_water__time_average_of_volume_flow_rate"