diff --git a/Project.toml b/Project.toml
index 0cbb5bd0c41e3fc4301c3f41a98b1bfad7d0d6b7..735281df498288ca8fca52cea46db76de523cada 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.5.0"
+version = "0.5.1"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 5677de6caab1b0252772ba174b895e1b0619d775..e189b75767e0829071736474160de81595c18332 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -27,15 +27,15 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
 [[ArnoldiMethod]]
 deps = ["LinearAlgebra", "Random", "StaticArrays"]
-git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
+git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
 uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
-version = "0.1.0"
+version = "0.2.0"
 
 [[ArrayInterface]]
 deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
-git-tree-sha1 = "e527b258413e0c6d4f66ade574744c94edef81f8"
+git-tree-sha1 = "265b06e2b1f6a216e0e8f183d28e4d354eab3220"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "3.1.40"
+version = "3.2.1"
 
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -67,9 +67,9 @@ version = "0.1.1"
 
 [[BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "61adeb0823084487000600ef8b1c00cc2474cd47"
+git-tree-sha1 = "365c0ea9a8d256686e97736d6b7fb0c880261a7a"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.2.0"
+version = "1.2.1"
 
 [[Bijections]]
 git-tree-sha1 = "705e7822597b432ebe152baa844b49f8026df090"
@@ -90,9 +90,9 @@ version = "0.1.6"
 
 [[ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "f885e7e7c124f8c92650d61b9477b9ac2ee607dd"
+git-tree-sha1 = "4c26b4e9e91ca528ea212927326ece5918a04b47"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.11.1"
+version = "1.11.2"
 
 [[ChangesOfVariables]]
 deps = ["LinearAlgebra", "Test"]
@@ -168,7 +168,7 @@ version = "4.0.4"
 [[CryoGrid]]
 deps = ["BenchmarkTools", "ComponentArrays", "ConstructionBase", "DataStructures", "Dates", "DiffEqBase", "DiffEqCallbacks", "DimensionalData", "Downloads", "ExprTools", "Flatten", "ForwardDiff", "IfElse", "InteractiveUtils", "Interpolations", "IntervalSets", "JSON3", "Lazy", "LinearAlgebra", "LoopVectorization", "ModelParameters", "OrdinaryDiffEq", "Parameters", "PreallocationTools", "ProgressMeter", "RecursiveArrayTools", "Reexport", "ReverseDiff", "RuntimeGeneratedFunctions", "Setfield", "SimulationLogs", "Statistics", "StructTypes", "SymbolicUtils", "Symbolics", "Test", "TimeSeries", "Unitful"]
 git-tree-sha1 = "45878a6ffaaf743c53141d40171a4efd675174dd"
-repo-rev = "HEAD"
+repo-rev = "refactor/setup-and-drivers"
 repo-url = ".."
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 version = "0.4.0"
@@ -210,9 +210,9 @@ uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
 [[DensityInterface]]
 deps = ["InverseFunctions", "Test"]
-git-tree-sha1 = "794daf62dce7df839b8ed446fc59c68db4b5182f"
+git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
 uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
-version = "0.3.3"
+version = "0.4.0"
 
 [[DiffEqBase]]
 deps = ["ArrayInterface", "ChainRulesCore", "DEDataArrays", "DataStructures", "Distributions", "DocStringExtensions", "FastBroadcast", "ForwardDiff", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "PreallocationTools", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"]
@@ -234,21 +234,21 @@ version = "1.0.3"
 
 [[DiffRules]]
 deps = ["LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "3287dacf67c3652d3fed09f4c12c187ae4dbb89a"
+git-tree-sha1 = "d8f468c5cd4d94e86816603f7d18ece910b4aaf1"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.4.0"
+version = "1.5.0"
 
 [[DimensionalData]]
 deps = ["Adapt", "ArrayInterface", "ConstructionBase", "Dates", "LinearAlgebra", "Random", "RecipesBase", "SparseArrays", "Statistics", "Tables"]
-git-tree-sha1 = "5184a92e568637ab0e985cdb76f426e97ff52126"
+git-tree-sha1 = "d6027664587200d6af27ea3fc3c18dca90dd983e"
 uuid = "0703355e-b756-11e9-17c0-8b28908087d0"
-version = "0.19.2"
+version = "0.19.5"
 
 [[Distances]]
-deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
-git-tree-sha1 = "837c83e5574582e07662bbbba733964ff7c26b9d"
+deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.10.6"
+version = "0.10.7"
 
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
@@ -256,9 +256,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 
 [[Distributions]]
 deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
-git-tree-sha1 = "ed2cf1e5388569cd00fa9cd28519434a5571633f"
+git-tree-sha1 = "7f3bec11f4bcd01bc1f507ebce5eadf1b0a78f47"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.26"
+version = "0.25.34"
 
 [[DocStringExtensions]]
 deps = ["LibGit2", "Markdown", "Pkg", "Test"]
@@ -290,15 +290,15 @@ version = "0.3.21"
 
 [[EllipsisNotation]]
 deps = ["ArrayInterface"]
-git-tree-sha1 = "9aad812fb7c4c038da7cab5a069f502e6e3ae030"
+git-tree-sha1 = "3fe985505b4b667e1ae303c9ca64d181f09d5c05"
 uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
-version = "1.1.1"
+version = "1.1.3"
 
 [[ExponentialUtilities]]
 deps = ["ArrayInterface", "LinearAlgebra", "Printf", "Requires", "SparseArrays"]
-git-tree-sha1 = "cb39752c2a1f83bbe0fda393c51c480a296042ad"
+git-tree-sha1 = "1b873816d2cfc8c0fcb1edcb08e67fdf630a70b7"
 uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
-version = "1.10.1"
+version = "1.10.2"
 
 [[ExprTools]]
 git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
@@ -360,6 +360,12 @@ version = "1.1.2"
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
+[[Graphs]]
+deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
+git-tree-sha1 = "92243c07e786ea3458532e199eb3feee0e7e08eb"
+uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
+version = "1.4.1"
+
 [[HostCPUFeatures]]
 deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
 git-tree-sha1 = "8f0dc80088981ab55702b04bba38097a44a1a3a9"
@@ -461,10 +467,10 @@ uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 version = "1.3.0"
 
 [[LabelledArrays]]
-deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
-git-tree-sha1 = "fa07d4ee13edf79a6ac2575ad28d9f43694e1190"
+deps = ["ArrayInterface", "ChainRulesCore", "LinearAlgebra", "MacroTools", "StaticArrays"]
+git-tree-sha1 = "3609bbf5feba7b22fb35fe7cb207c8c8d2e2fc5b"
 uuid = "2ee39098-c373-598a-b85f-a56591580800"
-version = "1.6.6"
+version = "1.6.7"
 
 [[Latexify]]
 deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
@@ -503,12 +509,6 @@ uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
-[[LightGraphs]]
-deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
-uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
-version = "1.3.5"
-
 [[LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
 git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd"
@@ -530,9 +530,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[LoopVectorization]]
 deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "Requires", "SIMDDualNumbers", "SLEEFPirates", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"]
-git-tree-sha1 = "2f8621145e87ce637a2f7ac701353c74d974455e"
+git-tree-sha1 = "9e10579c154f785b911d9ceb96c33fcc1a661171"
 uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
-version = "0.12.96"
+version = "0.12.99"
 
 [[MacroTools]]
 deps = ["Markdown", "Random"]
@@ -649,15 +649,15 @@ version = "1.4.1"
 
 [[OrdinaryDiffEq]]
 deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "LoopVectorization", "MacroTools", "MuladdMacro", "NLsolve", "Polyester", "PreallocationTools", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
-git-tree-sha1 = "6f76c887ddfd3f2a018ef1ee00a17b46bcf4886e"
+git-tree-sha1 = "758a3cc612baa4ec194dd2e71625d091edd8675a"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-version = "5.67.0"
+version = "5.68.0"
 
 [[PDMats]]
 deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "c8b8775b2f242c80ea85c83714c64ecfa3c53355"
+git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.3"
+version = "0.11.5"
 
 [[Parameters]]
 deps = ["OrderedCollections", "UnPack"]
@@ -740,9 +740,9 @@ uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
 version = "0.4.2"
 
 [[RecipesBase]]
-git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
+git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
 uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "1.1.2"
+version = "1.2.1"
 
 [[RecursiveArrayTools]]
 deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
@@ -774,10 +774,10 @@ uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.1.3"
 
 [[ReverseDiff]]
-deps = ["ChainRulesCore", "DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"]
-git-tree-sha1 = "14b09b7d28378e6e8be334edc0ade781b5c65e85"
+deps = ["ChainRulesCore", "DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "LogExpFunctions", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"]
+git-tree-sha1 = "2c7abf1c7b4f9bbfb0e755cf74d4bf0e43431b9e"
 uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-version = "1.10.0"
+version = "1.11.0"
 
 [[Rmath]]
 deps = ["Random", "Rmath_jll"]
@@ -819,9 +819,9 @@ version = "0.6.28"
 
 [[SciMLBase]]
 deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
-git-tree-sha1 = "ad2c7f08e332cc3bb05d33026b71fa0ef66c009a"
+git-tree-sha1 = "b3d23aa4e5f621b574b3b0d41c62c8624d27192a"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "1.19.4"
+version = "1.19.5"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -862,10 +862,10 @@ deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [[SparseDiffTools]]
-deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"]
-git-tree-sha1 = "0b0bd4086536520cfc452ba51b177be8cf634d91"
+deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "Graphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"]
+git-tree-sha1 = "f87076b43379cb0bd9f421cfe7c649fb510d8e4e"
 uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
-version = "1.17.2"
+version = "1.18.1"
 
 [[SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
@@ -896,21 +896,21 @@ deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[StatsAPI]]
-git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
+git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
 uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
-version = "1.0.0"
+version = "1.1.0"
 
 [[StatsBase]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
-git-tree-sha1 = "eb35dcc66558b2dda84079b9a1be17557d32091a"
+git-tree-sha1 = "2bb0cb32026a66037360606510fca5984ccc6b75"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.33.12"
+version = "0.33.13"
 
 [[StatsFuns]]
-deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
-git-tree-sha1 = "95072ef1a22b057b1e80f73c2a89ad238ae4cfff"
+deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
+git-tree-sha1 = "bedb3e17cc1d94ce0e6e66d3afa47157978ba404"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "0.9.12"
+version = "0.9.14"
 
 [[StrideArraysCore]]
 deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "Requires", "SIMDTypes", "Static", "ThreadingUtilities"]
@@ -930,15 +930,15 @@ uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[SymbolicUtils]]
 deps = ["AbstractTrees", "Bijections", "ChainRulesCore", "Combinatorics", "ConstructionBase", "DataStructures", "DocStringExtensions", "DynamicPolynomials", "IfElse", "LabelledArrays", "LinearAlgebra", "Metatheory", "MultivariatePolynomials", "NaNMath", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "TermInterface", "TimerOutputs"]
-git-tree-sha1 = "5a25c36acf4f288c88b4a83e3b87f0adc5e1f9f3"
+git-tree-sha1 = "5255e65d129c8edbde92fd2ede515e61098f93df"
 uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
-version = "0.18.0"
+version = "0.18.1"
 
 [[Symbolics]]
-deps = ["ConstructionBase", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "Metatheory", "NaNMath", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "TermInterface", "TreeViews"]
-git-tree-sha1 = "671779c01f26efbfc671edb924e085793ce301e6"
+deps = ["ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "IfElse", "Latexify", "Libdl", "LinearAlgebra", "MacroTools", "Metatheory", "NaNMath", "RecipesBase", "Reexport", "Requires", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArrays", "SymbolicUtils", "TermInterface", "TreeViews"]
+git-tree-sha1 = "1f738ebade1567d461c4db9ef15470a1fdf0b9b5"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "4.0.0"
+version = "4.1.1"
 
 [[TOML]]
 deps = ["Dates"]
@@ -1025,21 +1025,21 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [[Unitful]]
 deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
-git-tree-sha1 = "880f77d2cd4c6948e6bd55425b7b52f34dcd7f4b"
+git-tree-sha1 = "0992ed0c3ef66b0390e5752fe60054e5ff93b908"
 uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
-version = "1.9.1"
+version = "1.9.2"
 
 [[VectorizationBase]]
 deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "Hwloc", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"]
-git-tree-sha1 = "5239606cf3552aff43d79ecc75b1af1ce4625109"
+git-tree-sha1 = "17e5847bb36730d90801170ecd0ce4041a3dde86"
 uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
-version = "0.21.21"
+version = "0.21.22"
 
 [[VertexSafeGraphs]]
-deps = ["LightGraphs"]
-git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
+deps = ["Graphs"]
+git-tree-sha1 = "8351f8d73d7e880bfc042a8b6922684ebeafb35c"
 uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
-version = "0.1.2"
+version = "0.2.0"
 
 [[WoodburyMatrices]]
 deps = ["LinearAlgebra", "SparseArrays"]
diff --git a/src/Physics/Boundaries/Boundaries.jl b/src/Physics/Boundaries/Boundaries.jl
index ea2dfced2c15066a7522e73066f1aa53e495f066..0d768b2fc569bb4ab0438361e347413d85810451 100644
--- a/src/Physics/Boundaries/Boundaries.jl
+++ b/src/Physics/Boundaries/Boundaries.jl
@@ -33,7 +33,7 @@ struct Constant{S,T} <: BoundaryProcess
     Constant(::Type{S}, value::T) where {S<:BoundaryStyle,T} = new{S,T}(value)
 end
 ConstructionBase.constructorof(::Type{<:Constant{S}}) where {S} = value -> Constant(S,value)
-boundaryvalue(bc::Constant{S,T},l2,p2,s1,s2) where {S,T} = bc.value
+boundaryvalue(bc::Constant{S,T},l1,p2,l2,s1,s2) where {S,T} = bc.value
 
 BoundaryStyle(::Type{<:Constant{S}}) where {S} = S()
 
@@ -51,14 +51,14 @@ struct Periodic{S,T} <: BoundaryProcess
         new{S,T}(uconvert(u"s",period) |> dustrip, amplitude, phaseshift)
 end
 
-@inline boundaryvalue(bc::Periodic,l2,p2,s1,s2) = bc.amplitude*sin(Ï€*(1/bc.period)*t + bc.phaseshift)
+@inline boundaryvalue(bc::Periodic,l1,p2,l2,s1,s2) = bc.amplitude*sin(Ï€*(1/bc.period)*t + bc.phaseshift)
 
 BoundaryStyle(::Type{<:Periodic{S}}) where {S} = S()
 
 @with_kw struct Bias{P} <: BoundaryProcess
     bias::P = Param(0.0)
 end
-@inline boundaryvalue(bc::Bias,l2,p2,s1,s2) = bc.bias
+@inline boundaryvalue(bc::Bias,l1,p2,l2,s1,s2) = bc.bias
 
 BoundaryStyle(::Type{<:Bias}) = Dirichlet()
 
@@ -78,7 +78,7 @@ struct CombinedBoundaryProcess{B1,B2,F,S} <: BoundaryProcess
         new{B1,B2,F,typeof(BoundaryStyle(bc1))}(op,bc1,bc2)
     end
 end
-@inline boundaryvalue(cbc::CombinedBoundaryProcess,l2,p2,s1,s2) = cbc.op(cbc.bc1(l1,l2,p2,s1,s2), cbc.bc2(l1,l2,p2,s1,s2))
+@inline boundaryvalue(cbc::CombinedBoundaryProcess,l1,p2,l2,s1,s2) = cbc.op(cbc.bc1(l1,l2,p2,s1,s2), cbc.bc2(l1,l2,p2,s1,s2))
 variables(top::Top, cbc::CombinedBoundaryProcess) = tuplejoin(variables(top, cbc.bc1), variables(top, cbc.bc2))
 BoundaryStyle(::Type{CombinedBoundaryProcess{B1,B2,F,S}}) where {F,B1,B2,S} = S()
 # Overload arithmetic operators on boundary processes.
diff --git a/src/Physics/Boundaries/effects.jl b/src/Physics/Boundaries/effects.jl
index faccfb49d14ae9ef89992bb7171e12d9725b2fe6..af66fbbbc6f86f5352106541566993694e3aee76 100644
--- a/src/Physics/Boundaries/effects.jl
+++ b/src/Physics/Boundaries/effects.jl
@@ -15,6 +15,6 @@ function (damp::Damping)(u_top, u_sub, k_sub, Δsub, t)
         k_med = damp.k, # conductivity of damping medium (assumed uniform)
         k_sub = k_sub, # conductivity at grid center
         k = Numerics.harmonicmean(k_med, k_sub, d_med, d_sub);
-        -k*(u_sub - u_top)/d # flux (unscaled)
+        -k*(u_sub - u_top)/d/Δsub # flux per unit volume
     end
 end
diff --git a/src/Physics/HeatConduction/HeatConduction.jl b/src/Physics/HeatConduction/HeatConduction.jl
index d3ae3dc6cf430abb638379455629cd3c761378f9..ee4a2684fcb9bd3f17eeeefd46ec50fb2866d94d 100644
--- a/src/Physics/HeatConduction/HeatConduction.jl
+++ b/src/Physics/HeatConduction/HeatConduction.jl
@@ -49,7 +49,7 @@ Heat(::Val{(:Hâ‚›,:Hâ‚—)}; kwargs...) = Heat(;sp=PartitionedEnthalpy(), kwargs..
 Heat(::Val{:T}; kwargs...) = Heat(;sp=Temperature(), kwargs...)
 
 export enthalpy, heatcapacity, heatcapacity!, thermalconductivity, thermalconductivity!
-export heatconduction!, boundaryflux
+export heatconduction!
 include("heat.jl")
 export ConstantTemp, GeothermalHeatFlux, TemperatureGradient, NFactor, Damping
 include("heat_bc.jl")
diff --git a/src/Physics/HeatConduction/heat.jl b/src/Physics/HeatConduction/heat.jl
index 14e9b672b95fee7ce3dcff3c061d2ec14f4e3769..9b9830fad77bcbe8459412a616b21e42e7a58707 100644
--- a/src/Physics/HeatConduction/heat.jl
+++ b/src/Physics/HeatConduction/heat.jl
@@ -137,17 +137,11 @@ function prognosticstep!(::SubSurface, ::Heat{<:FreezeCurve,Temperature}, state)
     @inbounds @. state.dT = state.dH / state.Ceff
     return nothing
 end
-boundaryflux(::Neumann, top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, stop, ssub) =
-    @inbounds let δ₀ = Δ(ssub.grids.k)[1]
-        boundaryvalue(bc,sub,heat,stop,ssub)/δ₀
-    end
-boundaryflux(::Neumann, bot::Bottom, bc::BoundaryProcess, sub::SubSurface, heat::Heat, sbot, ssub) =
-    @inbounds let δ₀ = Δ(ssub.grids.k)[end]
-        boundaryvalue(bc,sub,heat,sbot,ssub)/δ₀
-    end
-function boundaryflux(::Dirichlet, top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, stop, ssub)
+boundaryflux(::Neumann, bc::BoundaryProcess, top::Top, heat::Heat, sub::SubSurface, stop, ssub) = boundaryvalue(bc,top,heat,sub,stop,ssub)
+boundaryflux(::Neumann, bc::BoundaryProcess, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub) = boundaryvalue(bc,bot,heat,sub,sbot,ssub)
+function boundaryflux(::Dirichlet, bc::BoundaryProcess, top::Top, heat::Heat, sub::SubSurface, stop, ssub)
     Δk = Δ(ssub.grids.k)
-    @inbounds let Tupper=boundaryvalue(bc,sub,heat,stop,ssub),
+    @inbounds let Tupper=boundaryvalue(bc,top,heat,sub,stop,ssub),
         Tsub=ssub.T[1],
         k=ssub.k[1],
         d=Δk[1],
@@ -155,9 +149,9 @@ function boundaryflux(::Dirichlet, top::Top, bc::BoundaryProcess, sub::SubSurfac
         -k*(Tsub-Tupper)/δ₀/d
     end
 end
-function boundaryflux(::Dirichlet, bot::Bottom, bc::BoundaryProcess, sub::SubSurface, heat::Heat, sbot, ssub)
+function boundaryflux(::Dirichlet, bc::BoundaryProcess, bot::Bottom, heat::Heat, sub::SubSurface, sbot, ssub)
     Δk = Δ(ssub.grids.k)
-    @inbounds let Tlower=boundaryvalue(bc,sub,heat,sbot,ssub),
+    @inbounds let Tlower=boundaryvalue(bc,bot,heat,sub,sbot,ssub),
         Tsub=ssub.T[end],
         k=ssub.k[end],
         d=Δk[1],
@@ -173,7 +167,7 @@ function interact!(top::Top, bc::BoundaryProcess, sub::SubSurface, heat::Heat, s
     # assumes (1) k has already been computed, (2) surface conductivity = cell conductivity
     @inbounds ssub.k[1] = ssub.k[2]
     # boundary flux
-    @inbounds ssub.dH[1] += boundaryflux(top, bc, sub, heat, stop, ssub)
+    @inbounds ssub.dH[1] += boundaryflux(bc, top, heat, sub, stop, ssub)
     return nothing # ensure no allocation
 end
 """
@@ -184,7 +178,7 @@ function interact!(sub::SubSurface, heat::Heat, bot::Bottom, bc::BoundaryProcess
     # assumes (1) k has already been computed, (2) bottom conductivity = cell conductivity
     @inbounds ssub.k[end] = ssub.k[end-1]
     # boundary flux
-    @inbounds ssub.dH[end] += boundaryflux(bot, bc, sub, heat, sbot, ssub)
+    @inbounds ssub.dH[end] += boundaryflux(bc, bot, heat, sub, sbot, ssub)
     return nothing # ensure no allocation
 end
 """
diff --git a/src/Physics/HeatConduction/heat_bc.jl b/src/Physics/HeatConduction/heat_bc.jl
index ad2eee8e70c7661520ab6e4626c90dd483de9183..35185cd4823c9f8f7da8ab6f5c61747891549380 100644
--- a/src/Physics/HeatConduction/heat_bc.jl
+++ b/src/Physics/HeatConduction/heat_bc.jl
@@ -18,7 +18,7 @@ BoundaryStyle(::Type{<:TemperatureGradient{<:Damping}}) = Neumann()
     winterfactor::W = Param(1.0, bounds=(0.0,1.0)) # applied when Tair <= 0
     summerfactor::S = Param(1.0, bounds=(0.0,1.0)) # applied when Tair > 0
 end
-@inline function boundaryvalue(bc::TemperatureGradient{<:NFactor},l2,p2::Heat,s1,s2)
+@inline function boundaryvalue(bc::TemperatureGradient{<:NFactor}, l1, ::Heat, l2, s1, s2)
     let nfw = bc.effect.winterfactor,
         nfs = bc.effect.summerfactor,
         Tair = bc.T(s1.t),
@@ -27,7 +27,7 @@ end
     end
 end
 # damped temperature gradient
-@inline function boundaryvalue(bc::TemperatureGradient{<:Damping},l2,p2::Heat,s1,s2)
+@inline function boundaryvalue(bc::TemperatureGradient{<:Damping}, l1, ::Heat, l2, s1,s2)
     let Ttop = bc.T(s1.t),
         Tsub = s2.T[1],
         δ = Δ(s2.grids.k)[1], # grid cell size
diff --git a/src/methods.jl b/src/methods.jl
index 7bbd7853bd3ffe0d2b185057f9427b1e50d26ac7..c403ed9c0b40ddf7dced8a0a4c6f7e1ebe0e80a3 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -34,27 +34,25 @@ prognosticstep!(l::Layer, p::Process, state) = error("no prognostic step defined
 """
     interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2)
 
-Defines a boundary interaction between two processes on adjacent layers.
+Defines a boundary interaction between two processes on adjacent layers. For any interaction, the order of the arguments
+follows decreasing depth, i.e. the first layer/process is always on top of the second layer/process. This ordering matters
+and separate dispatches must be provided for interactions in reverse order.
 """
 interact!(::Layer, ::Process, ::Layer, ::Process, state1, state2) = nothing
 """
-    boundaryflux(top::Top, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2)
-    boundaryflux(bot::Bottom, bc::BoundaryProcess, sub::SubSurface, bc1::SubSurfaceProcess, s1, s2)
-    boundaryflux(s::BoundaryStyle, ::Top, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2)
-    boundaryflux(s::BoundaryStyle, sub::SubSurface, p::SubSurfaceProcess, ::Bottom, bc::BoundaryProcess, s1, s2)
+    boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
+    boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2)
 
 Computes the flux dH/dt at the boundary layer. Calls boundaryflux(BoundaryStyle(B),...) to allow for generic implementations by boundary condition type.
 """
-boundaryflux(top::Top, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2) = boundaryflux(BoundaryStyle(bc), top, bc, sub, p, s1, s2)
-boundaryflux(bot::Bottom, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2) = boundaryflux(BoundaryStyle(bc), bot, bc, sub, p, s1, s2)
-boundaryflux(s::BoundaryStyle, ::Top, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2) = error("missing implementation of $(typeof(s)) top boundaryflux for $(typeof(bc)) + $(typeof(p)) on $(typeof(sub))")
-boundaryflux(s::BoundaryStyle, ::Bottom, bc::BoundaryProcess, sub::SubSurface, p::SubSurfaceProcess, s1, s2) = error("missing implementation of $(typeof(s)) bottom boundaryflux for $(typeof(bc)) + $(typeof(p)) on $(typeof(sub))")
+boundaryflux(bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = boundaryflux(BoundaryStyle(bc), bc, b, p, sub, s1, s2)
+boundaryflux(s::BoundaryStyle, bc::BoundaryProcess, b::Union{Top,Bottom}, p::SubSurfaceProcess, sub::SubSurface, s1, s2) = error("missing implementation of $(typeof(b)) $(typeof(s)) boundaryflux for $(typeof(bc)) + $(typeof(p)) on $(typeof(sub))")
 """
-    boundaryvalue(bc::BoundaryProcess, layer::SubSurface, p::SubSurfaceProcess, sbc, ssub)
+    boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub)
 
 Computes the value of the boundary condition specified by `bc` for the given layer/process combinations.
 """
-boundaryvalue(bc::BoundaryProcess, layer, p, sbc, ssub) = error("missing implementation of boundaryvalue for $(typeof(bc)) on $(typeof(layer)) with $(typeof(p))")
+boundaryvalue(bc::BoundaryProcess, b, p, layer, sbc, ssub) = error("missing implementation of boundaryvalue for $(typeof(b)) $(typeof(bc)) on $(typeof(layer)) with $(typeof(p))")
 """
     observe(::Val{name}, ::Layer, ::Process, state1)
 
diff --git a/src/types.jl b/src/types.jl
index 671103d3ea67d66bb52d276c21adb8fdd069da79..50b63f92e5ed60e88b200e04f223b80a7d2d9a33 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -77,7 +77,7 @@ Base.Broadcast.broadcastable(p::Process) = Ref(p)
 """
 Trait that specifies the "style" or kind of boundary condition. This can be used to write generic
 implementations of `interact!` that are (relatively) agnostic to specific implementations of
-`BoundaryProcess`. A good example of this can be found in `HeatConduction.boundaryflux`.
+`BoundaryProcess`. A good example of this can be found in the `boundaryflux` method interface.
 """
 abstract type BoundaryStyle end
 """
diff --git a/test/Physics/HeatConduction/heat_conduction_tests.jl b/test/Physics/HeatConduction/heat_conduction_tests.jl
index ad593dfd9fc56be3076affec247c00d068667af9..50cffb0650f03f2dd541e521522b5e54c7a08c7f 100644
--- a/test/Physics/HeatConduction/heat_conduction_tests.jl
+++ b/test/Physics/HeatConduction/heat_conduction_tests.jl
@@ -37,15 +37,15 @@ include("../../types.jl")
 			T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
 			∂H = zeros(length(T₀))u"J/s/m^3"
 			state = (T=T₀,k=k,dH=∂H,grids=(T=xc,k=x),t=0.0)
-			@test boundaryflux(Top(),bc,sub,heat,state,state) > 0.0u"J/s/m^3"
-			@test boundaryflux(Bottom(),bc,sub,heat,state,state) < 0.0u"J/s/m^3"
+			@test boundaryflux(bc,Top(),heat,sub,state,state) > 0.0u"J/s/m^3"
+			@test boundaryflux(bc,Bottom(),heat,sub,state,state) < 0.0u"J/s/m^3"
 		end
 		@testset "top: -, bot: +" begin
 			T₀ = Vector(LinRange(27,-23,length(xc)))u"°C"
 			∂H = zeros(length(T₀))u"J/s/m^3"
 			state = (T=T₀,k=k,dH=∂H,grids=(T=xc,k=x),t=0.0)
-			@test boundaryflux(Top(),bc,sub,heat,state,state) < 0.0u"J/s/m^3"
-			@test boundaryflux(Bottom(),bc,sub,heat,state,state) > 0.0u"J/s/m^3"
+			@test boundaryflux(bc,Top(),heat,sub,state,state) < 0.0u"J/s/m^3"
+			@test boundaryflux(bc,Bottom(),heat,sub,state,state) > 0.0u"J/s/m^3"
 		end
 		@testset "inner edge boundary (positive)" begin
 			T₀ = Vector(sin.(ustrip.(xc).*π))u"°C"
@@ -61,6 +61,13 @@ include("../../types.jl")
 			@test ∂H[1] < 0.0u"J/s/m^3"
 			@test ∂H[end] < 0.0u"J/s/m^3"
 		end
+		@testset "Neumann boundary" begin
+			bc = Constant(Neumann, -1.0u"J/m^3")
+			T₀ = Vector(LinRange(-23,27,length(xc)))u"°C"
+			∂H = zeros(length(T₀))u"J/s/m^3"
+			state = (T=T₀,k=k,dH=∂H,grids=(T=xc,k=x),t=0.0)
+			@test boundaryflux(bc,Top(),heat,sub,state,state) == -1.0u"J/m^3"
+		end
 	end
 end
 @testset "Boundary conditions" begin
@@ -70,7 +77,7 @@ end
 		tgrad = TemperatureGradient(forcing, NFactor(winterfactor=0.5, summerfactor=1.0))
 		sub = TestGroundLayer()
 		heat = Heat()
-		f1(t) = boundaryvalue(tgrad,sub,heat,(t=t,),(t=t,))
+		f1(t) = boundaryvalue(tgrad,Top(),heat,sub,(t=t,),(t=t,))
 		Tres = f1.(Dates.datetime2epochms.(ts)./1000.0)
 		@test all(Tres .≈ [1.0,0.5,-0.25,-0.5,0.1])
 	end
@@ -95,8 +102,8 @@ end
 		# compute boundary fluxes;
 		# while not correct in general, for this test case we can just re-use state for both layers.
 		state = (T=T_K,k=k,dH=dT,grids=(T=xc,k=x),t=t)
-		dT[1] += boundaryflux(Top(),bc,sub,heat,state,state)
-		dT[end] += boundaryflux(Bottom(),bc,sub,heat,state,state)
+		dT[1] += boundaryflux(bc,Top(),heat,sub,state,state)
+		dT[end] += boundaryflux(bc,Bottom(),heat,sub,state,state)
 		# strip units from dT before returning dT to the solver
 		return ustrip.(dT)
 	end