diff --git a/README.md b/README.md
index 7d41667a47a7faf284043a246fdfef29421a265b..2127659debce0107841bfd6d9bc5a06c39ef59e0 100755
--- a/README.md
+++ b/README.md
@@ -51,12 +51,12 @@ using Plots: plot
 # The forcing file will be automatically downloaded to the input/ folder if not already present.
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
 # get preset soil and initial temperature profile for Samoylov
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
 initT = initializer(:T, tempprofile)
 # choose grid with 5cm spacing
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
 # build a default Tile with heat conduction on the given soil profile
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
     TemperatureBC(forcings.Tair),
     GeothermalHeatFlux(0.053u"W/m^2"),
     soilprofile,
diff --git a/docs/src/manual/architecture.md b/docs/src/manual/architecture.md
index eb0720e2b34c9a70854733c0882261822714f621..5cadab0602ac9a2db344c2128fb1e4f708326333 100644
--- a/docs/src/manual/architecture.md
+++ b/docs/src/manual/architecture.md
@@ -65,7 +65,7 @@ using CryoGrid
 using CryoGrid.Diagnostics
 
 soil = Ground()
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
 state = Diagnostics.build_dummy_state(grid, soil)
 
 @which CryoGrid.computediagnostic!(soil, state)
diff --git a/docs/src/manual/overview.md b/docs/src/manual/overview.md
index 03357710a8154c4c8cd6f447ced8154f8f156534..414f950a831bd3945b6014483e4099e7e641bd1b 100644
--- a/docs/src/manual/overview.md
+++ b/docs/src/manual/overview.md
@@ -18,7 +18,7 @@ strat = Stratigraphy(
     0.0u"m" => Ground(soilprofile, HeatBalance(:H; freezecurve=DallAmico())),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
 );
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
 # define initial conditions for temperature using a given profile;
 # The default initializer linearly interpolates between profile points.
 initT = initializer(:T, tempprofile)
diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md
index 2d887e003a9e1d0a1ee98299b674b802e6694fa9..75c28cb83857ba96b4ceb0ec547dd514a05d3b38 100644
--- a/docs/src/quickstart.md
+++ b/docs/src/quickstart.md
@@ -1,6 +1,6 @@
 ## Quick start
 
-After [installing](installation.md) `CryoGrid.jl`, you can get started right away with a simple soil heat model. The [`Presets`](@ref) module (aliased `CryoGrid.Presets`) provides pre-specified model configurations that can be obtained with a single function call. It is also possible to provide a custom soil and temperature profile to `SoilHeatTile`; here the values given in `Presets.SamoylovDefault` is used.
+After [installing](installation.md) `CryoGrid.jl`, you can get started right away with a simple soil heat model. The [`Presets`](@ref) module (aliased `CryoGrid.Presets`) provides pre-specified model configurations that can be obtained with a single function call. It is also possible to provide a custom soil and temperature profile to `SoilHeatTile`; here the values given in `SamoylovDefault` is used.
 
 ```julia
 using CryoGrid
@@ -10,12 +10,12 @@ using Plots
 # The forcing file will be automatically downloaded to the input/ folder if not already present.
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044, :Tair => u"°C");
 # get preset soil and initial temperature profile for Samoylov
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
 initT = initializer(:T, tempprofile)
 # choose grid with 5cm spacing
-grid = CryoGrid.Presets.DefaultGrid_5cm
+grid = CryoGrid.DefaultGrid_5cm
 # build Tile from the given soil and temperature profiles
-tile = CryoGrid.Presets.SoilHeatTile(TemperatureBC(forcings.Tair), GeothermalHeatFlux(0.053u"W/m^2"), soilprofile, initT, grid=grid)
+tile = CryoGrid.SoilHeatTile(TemperatureBC(forcings.Tair), GeothermalHeatFlux(0.053u"W/m^2"), soilprofile, initT, grid=grid)
 # define time span (1 year)
 tspan = (DateTime(2010,11,30),DateTime(2011,11,30))
 u0, du0 = initialcondition!(tile, tspan)
diff --git a/examples/cglite_parameter_ensembles.jl b/examples/cglite_parameter_ensembles.jl
index 91bb16c8143935a083c52e3a29916cba06831d0c..4c6c40e53418828e36ff15dae11ee1ffe42779f2 100644
--- a/examples/cglite_parameter_ensembles.jl
+++ b/examples/cglite_parameter_ensembles.jl
@@ -42,7 +42,7 @@ strat = Stratigraphy(
     soil_layers,
     z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+modelgrid = CryoGrid.DefaultGrid_2cm
 tile = Tile(strat, modelgrid, ssinit);
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
 # Here we specify a time span of 10 years.
diff --git a/examples/heat_freeW_bucketW_samoylov.jl b/examples/heat_freeW_bucketW_samoylov.jl
index 96d6762546aca339005bf94d9726ef2b6b8ba548..355b734e208a175d9dfc6dcf88198eed22599802 100644
--- a/examples/heat_freeW_bucketW_samoylov.jl
+++ b/examples/heat_freeW_bucketW_samoylov.jl
@@ -6,7 +6,7 @@
 # Frist, load forcings and define boundary conditions.
 using CryoGrid
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
-_, tempprofile = CryoGrid.Presets.SamoylovDefault;
+_, tempprofile = CryoGrid.SamoylovDefault;
 initT = initializer(:T, tempprofile)
 initsat = initializer(:sat, (l,state) -> state.sat .= l.para.sat);
 upperbc = WaterHeatBC(
@@ -23,7 +23,7 @@ strat = @Stratigraphy(
     2.0u"m" => Ground(SimpleSoil(por=0.10,sat=1.0,org=0.0), heat=HeatBalance(), water=WaterBalance(BucketScheme())),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2")),
 );
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+modelgrid = CryoGrid.DefaultGrid_2cm
 tile = Tile(strat, modelgrid, forcings, initT, initsat);
 
 # Now we set up the problem and solve using the integrator interface.
diff --git a/examples/heat_freeW_lake_lite_implicit.jl b/examples/heat_freeW_lake_lite_implicit.jl
index c06904729a386e354d93503f19d70dcf37d8d2d5..dc4d2da8df057789a3e340f0be33ab03579a4254 100644
--- a/examples/heat_freeW_lake_lite_implicit.jl
+++ b/examples/heat_freeW_lake_lite_implicit.jl
@@ -20,7 +20,7 @@ tempprofile_linear = TemperatureProfile(
     10.0u"m" => -10.0u"°C", 
     1000.0u"m" => 10.2u"°C"
 )
-modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm))
+modelgrid = Grid(vcat(-1.0u"m":0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
 z_top = -1.0u"m"
 z_sub = keys(soilprofile)
 z_bot = modelgrid[end]
diff --git a/examples/heat_freeW_lite_implicit.jl b/examples/heat_freeW_lite_implicit.jl
index eb5b993126799fab3fe36b4c50aa12d5338be341..aa88e5574dd0859035b4ec239ea82f74d14e46ad 100644
--- a/examples/heat_freeW_lite_implicit.jl
+++ b/examples/heat_freeW_lite_implicit.jl
@@ -36,7 +36,7 @@ strat = @Stratigraphy(
     soil_layers...,
     z_bot => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
-modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.Presets.DefaultGrid_2cm))
+modelgrid = Grid(vcat(z_top:0.02u"m":-0.02u"m", CryoGrid.DefaultGrid_2cm))
 tile = Tile(strat, modelgrid, forcings, ssinit);
 
 # Since the solver can take daily timesteps, we can easily specify longer simulation time spans at minimal cost.
diff --git a/examples/heat_freeW_samoylov.jl b/examples/heat_freeW_samoylov.jl
index e2eb192b1713ba630d7b8b914b3dbe0bd683b936..a54e08068be210bfbd23febabb54be83d0e00c5b 100644
--- a/examples/heat_freeW_samoylov.jl
+++ b/examples/heat_freeW_samoylov.jl
@@ -23,14 +23,14 @@ soilprofile = SoilProfile(
 );
 
 # We construct a state variable initializer for temperature `T` from the temperature profile preset for Samoylov.
-initT = initializer(:T, CryoGrid.Presets.SamoylovDefault.tempprofile);
+initT = initializer(:T, CryoGrid.SamoylovDefault.tempprofile);
 
 # We choose the default grid discretization with 5 cm spacing at the surface.
 grid = CryoGrid.DefaultGrid_5cm;
 
 # Now we construct the Tile using the built-in model configuration `SoilHeatTile` which defines a
 # standalone soil straigraphy with only heat conduction and no water flow.
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
     :H,
     TemperatureBC(Input(:Tair), NFactor(nf=Param(0.6), nt=Param(0.9))),
     GeothermalHeatFlux(0.053u"W/m^2"),
@@ -59,7 +59,4 @@ out = CryoGridOutput(sol)
 import Plots
 zs = [1,10,20,30,50,100,200,500,1000]u"cm"
 cg = Plots.cgrad(:copper,rev=true);
-Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
-
-using BenchmarkTools
-@profview @btime $tile($du0, $u0, $prob.p, $prob.tspan[1])
+Plots.plot(out.T[Z(Near(zs))], color=cg[LinRange(0.0,1.0,length(zs))]', ylabel="Temperature", leg=false, size=(800,500), dpi=150)
\ No newline at end of file
diff --git a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
index 4248173e186576545ba3f53daf17cc8f10d13d80..2f47c9d5f74688b974b886933a78f7d3eb960991 100644
--- a/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
+++ b/examples/heat_freeW_seb_snow_bucketW_samoylov.jl
@@ -8,7 +8,7 @@ using CryoGrid
 using OrdinaryDiffEq
 
 # First, load the forcings and construct the Tile.
-modelgrid = CryoGrid.Presets.DefaultGrid_2cm;
+modelgrid = CryoGrid.DefaultGrid_2cm;
 soilprofile = SoilProfile(
     0.0u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.75),
     0.1u"m" => SimpleSoil(por=0.80,sat=1.0,org=0.25),
@@ -17,9 +17,9 @@ soilprofile = SoilProfile(
     10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0),
 );
 ## mid-winter temperature profile
-tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
+tempprofile = CryoGrid.SamoylovDefault.tempprofile
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
-tempprofile = CryoGrid.Presets.SamoylovDefault.tempprofile
+tempprofile = CryoGrid.SamoylovDefault.tempprofile
 initT = initializer(:T, tempprofile)
 seb = SurfaceEnergyBalance()
 swb = SurfaceWaterBalance()
diff --git a/examples/heat_freeW_snow_samoylov.jl b/examples/heat_freeW_snow_samoylov.jl
index 2a42322224ff9297e48c922294b4bd661c5632d1..2bdc5d3e603f3bf227e5fa890f68378ec9f6c9f9 100644
--- a/examples/heat_freeW_snow_samoylov.jl
+++ b/examples/heat_freeW_snow_samoylov.jl
@@ -15,7 +15,7 @@ soilprofile = SoilProfile(
     3.0u"m" => SimpleSoil(por=0.50,sat=1.0,org=0.0),
     10.0u"m" => SimpleSoil(por=0.30,sat=1.0,org=0.0),
 );
-initT = initializer(:T, CryoGrid.Presets.SamoylovDefault.tempprofile)
+initT = initializer(:T, CryoGrid.SamoylovDefault.tempprofile)
 initsat = initializer(:sat, 1.0)
 z_top = -2.0u"m"
 z_sub = keys(soilprofile)
@@ -39,7 +39,7 @@ strat = @Stratigraphy(
     ground_layers...,
     z_bot => Bottom(GeothermalHeatFlux(0.053u"J/s/m^2"))
 );
-modelgrid = CryoGrid.Presets.DefaultGrid_5cm
+modelgrid = CryoGrid.DefaultGrid_5cm
 tile = Tile(strat, modelgrid, forcings, initT, initsat)
 # define time span, 2 years + 3 months
 tspan = (DateTime(2010,9,30), DateTime(2012,9,30))
diff --git a/examples/heat_sfcc_constantbc.jl b/examples/heat_sfcc_constantbc.jl
index 11cee1160fc93ed5e5008d37d4444574b8f78c9a..6cc3f12c378449ba837f2ee01c268ac4a4ea3288 100644
--- a/examples/heat_sfcc_constantbc.jl
+++ b/examples/heat_sfcc_constantbc.jl
@@ -8,7 +8,7 @@
 using CryoGrid
 
 # Select default grid and initial temperature profile.
-grid = CryoGrid.Presets.DefaultGrid_10cm
+grid = CryoGrid.DefaultGrid_10cm
 tempprofile = TemperatureProfile(
     0.0u"m" => -20.0u"°C",
     1000.0u"m" => 10.0u"°C",
@@ -31,7 +31,7 @@ Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc)
 # specified.
 heatop = Heat.Diffusion1D(:H)
 initT = initializer(:T, tempprofile)
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
     heatop,
     ## 10 W/m^2 in and out, i.e. net zero flux
     ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
diff --git a/examples/heat_sfcc_richardseq_samoylov.jl b/examples/heat_sfcc_richardseq_samoylov.jl
index 36aada3402318d30617c5748d24f5b211be30c5f..877d07f80adc593d20e8f6247ad07a6fe0fc1214 100644
--- a/examples/heat_sfcc_richardseq_samoylov.jl
+++ b/examples/heat_sfcc_richardseq_samoylov.jl
@@ -38,7 +38,7 @@ strat = @Stratigraphy(
     2.0u"m" => Ground(SimpleSoil(por=0.10,sat=1.0,org=0.0,freezecurve=sfcc); heat, water),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
-grid = CryoGrid.Presets.DefaultGrid_2cm
+grid = CryoGrid.DefaultGrid_2cm
 tile = Tile(strat, grid, forcings, initT);
 u0, du0 = @time initialcondition!(tile, tspan)
 prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600, savevars=(:T,:θw,:θwi,:kw));
diff --git a/examples/heat_sfcc_salt_constantbc.jl b/examples/heat_sfcc_salt_constantbc.jl
index bfc22d37b9970456d4b82f8eeef5ffbd46353250..e686b01359d1c4c45c69d2ea091d8a91d3f43c52 100644
--- a/examples/heat_sfcc_salt_constantbc.jl
+++ b/examples/heat_sfcc_salt_constantbc.jl
@@ -14,7 +14,7 @@ strat = @Stratigraphy(
     0.0u"m" => SalineGround(SimpleSoil(freezecurve=sfcc), heat=HeatBalance(:T), salt=SaltMassBalance()),
     1000.0u"m" => Bottom(GeothermalHeatFlux(0.053u"W/m^2"))
 );
-modelgrid = CryoGrid.Presets.DefaultGrid_10cm
+modelgrid = CryoGrid.DefaultGrid_10cm
 tile = Tile(strat, modelgrid, initT, initsalt, initpor)
 tspan = (DateTime(1990,1,1),DateTime(2000,12,31))
 u0, du0 = initialcondition!(tile, tspan)
diff --git a/examples/heat_sfcc_samoylov.jl b/examples/heat_sfcc_samoylov.jl
index 3b6b81192cb4a7f9177701624e1851fcc295dd73..add6edc0b85956aa87824933299e223342503a08 100644
--- a/examples/heat_sfcc_samoylov.jl
+++ b/examples/heat_sfcc_samoylov.jl
@@ -5,14 +5,15 @@
 
 using CryoGrid
 
-# First we set up the soil heat model:
+# First we set up the soil heat model. Note that the default soil profile
+# for Samoylov already includes the appropriate freeze curves.
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA5_fitted_daily_1979_2020);
-grid = CryoGrid.Presets.DefaultGrid_5cm
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
+grid = CryoGrid.DefaultGrid_5cm
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
 initT = initializer(:T, tempprofile)
 upperbc = TemperatureBC(Input(:Tair), NFactor(nf=0.6, nt=0.9))
 lowerbc = GeothermalHeatFlux(0.053u"W/m^2")
-tile = CryoGrid.Presets.SoilHeatTile(upperbc, lowerbc, soilprofile, forcings, initT; grid=grid)
+tile = CryoGrid.SoilHeatTile(upperbc, lowerbc, soilprofile, forcings, initT; grid=grid)
 tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
 u0, du0 = @time initialcondition!(tile, tspan)
 prob = CryoGridProblem(tile, u0, tspan, saveat=3*3600.0, savevars=(:T,));
diff --git a/examples/heat_simple_autodiff_grad.jl b/examples/heat_simple_autodiff_grad.jl
index 5b88521e0d56bae9f67a88819fa8cc649ab7f689..99f74eedc7339b226b8549a1c2da431fe9ec1722 100644
--- a/examples/heat_simple_autodiff_grad.jl
+++ b/examples/heat_simple_autodiff_grad.jl
@@ -7,10 +7,10 @@
 # Set up forcings and boundary conditions similarly to other examples:
 using CryoGrid
 forcings = loadforcings(CryoGrid.Forcings.Samoylov_ERA_obs_fitted_1979_2014_spinup_extended_2044);
-soilprofile, tempprofile = CryoGrid.Presets.SamoylovDefault
-grid = CryoGrid.Presets.DefaultGrid_5cm
+soilprofile, tempprofile = CryoGrid.SamoylovDefault
+grid = CryoGrid.DefaultGrid_5cm
 initT = initializer(:T, tempprofile)
-tile = CryoGrid.Presets.SoilHeatTile(
+tile = CryoGrid.SoilHeatTile(
     :T,
     TemperatureBC(Input(:Tair), NFactor(nf=Param(0.5), nt=Param(0.9))),
     GeothermalHeatFlux(0.053u"W/m^2"),
diff --git a/src/CryoGrid.jl b/src/CryoGrid.jl
index bad4afaa280860cd99f62a43c8eed3a23314875d..1f13ec07ca63536f95f1356b6119ec308bab866d 100755
--- a/src/CryoGrid.jl
+++ b/src/CryoGrid.jl
@@ -36,6 +36,8 @@ import Interpolations
 @reexport using Unitful
 @reexport using UnPack
 
+using SciMLBase: intervals # resolve name collision
+
 export Interpolations
 
 # Common types and methods
@@ -79,6 +81,7 @@ include("IO/InputOutput.jl")
 
 include("Tiles/Tiles.jl")
 @reexport using .Tiles
+using .Tiles: layers # resolve name collision
 parameters = Tiles.parameters
 
 export ConstantBC, PeriodicBC, ConstantValue, PeriodicValue, ConstantFlux, PeriodicFlux
@@ -99,8 +102,9 @@ include("problem.jl")
 export CGEuler, CryoGridIntegrator, CryoGridSolution
 include("Solvers/Solvers.jl")
 
-# Presets
-include("Presets/Presets.jl")
+# Built-in model definitions
+export SoilHeatTile, SamoylovDefault
+include("models.jl")
 
 using .InputOutput: Resource
 
diff --git a/src/IO/forcings/forcings_json.jl b/src/IO/forcings/forcings_json.jl
index 3c56b0d8593fab94e1629d3aa9cd032b73992de3..0a19ff877e2597d3dcd632f24d2bc73a06223015 100644
--- a/src/IO/forcings/forcings_json.jl
+++ b/src/IO/forcings/forcings_json.jl
@@ -41,7 +41,7 @@ function loadforcings(format::ForcingFormatJSON{1}, filename::String, kwargs...)
     end
     return Interpolated1D(DimStack((; forcings...)); kwargs...)
 end
-function loadforcings(format::ForcingFormatJSON{2}, filename::String)
+function loadforcings(format::ForcingFormatJSON{2}, filename::String; kwargs...)
     dict = open(filename, "r") do file; JSON3.read(file) end
     # convert JSON3 dict for data field to Julia dict
     data = Dict(dict[:data]...)
diff --git a/src/Physics/Soils/para/simple.jl b/src/Physics/Soils/para/simple.jl
index 213b902b0a277d1675f925a17a2afb111d3f6ce5..13c7b82543e16210fd9e9ebd46019978aa0fcf78 100644
--- a/src/Physics/Soils/para/simple.jl
+++ b/src/Physics/Soils/para/simple.jl
@@ -35,12 +35,12 @@ SoilThermalProperties(
     kh_i = ThermalProperties().kh_i,
     kh_a = ThermalProperties().kh_a,
     kh_o=Param(0.25, units=u"W/m/K", domain=0..Inf), # organic [Hillel (1982)]
-    kh_m=Param(3.8, untis=u"W/m/K", domain=0..Inf), # mineral [Hillel (1982)]
+    kh_m=Param(3.8, units=u"W/m/K", domain=0..Inf), # mineral [Hillel (1982)]
     ch_w = ThermalProperties().ch_w,
     ch_i = ThermalProperties().ch_i,
     ch_a = ThermalProperties().ch_a,
     ch_o=Param(2.5e6, units=u"J/K/m^3", domain=0..Inf), # heat capacity organic
-    ch_m=Param(2.0e6, untis=u"J/K/m^3", domain=0..Inf), # heat capacity mineral
+    ch_m=Param(2.0e6, units=u"J/K/m^3", domain=0..Inf), # heat capacity mineral
     kwargs...,
 ) = ThermalProperties(; kh_w, kh_i, kh_a, kh_m, kh_o, ch_w, ch_i, ch_a, ch_m, ch_o, kwargs...)
 
diff --git a/src/Physics/Surface/SEB/seb.jl b/src/Physics/Surface/SEB/seb.jl
index 4a98aa7b9b0fa9e8fb37fd1b1c9fd9f003f20524..b5245cdb2e20edef8f3b91fd519ebd8f31f49e2f 100644
--- a/src/Physics/Surface/SEB/seb.jl
+++ b/src/Physics/Surface/SEB/seb.jl
@@ -96,7 +96,7 @@ function SurfaceEnergyBalance(;
     stabfun::StabilityFunctions = HøgstrømSHEBA(),
 )
     inputs = (; Tair, Lin, Sin, pr, qh, wind)
-    SurfaceEnergyBalance(; inputs, para, solscheme, stabfun)
+    SurfaceEnergyBalance(inputs, para, solscheme, stabfun)
 end
 
 """
diff --git a/src/Physics/Surface/SEB/seb_state.jl b/src/Physics/Surface/SEB/seb_state.jl
index d2703a4887ed8bcfdb12040756ea06976a66aa20..972179688cfe6f8cccbfc129ea896815eff46a8c 100644
--- a/src/Physics/Surface/SEB/seb_state.jl
+++ b/src/Physics/Surface/SEB/seb_state.jl
@@ -1,10 +1,10 @@
 """
-    SEBInputs{TT1,TT2,TR,TP,TQ,TW,TZ,Tsurf}
+    SEBInputs{TT1,TT2,TR,TP,TQ,TW,Tsurf}
 
 Non-prognostic input variables to the SEB that are assumed to be given.
 """
-Base.@kwdef struct SEBInputs{TT1,TT2,TR,TP,TQ,TW,TZ,Tsurf}
-    Ts::TT1 = nothing
+Base.@kwdef struct SEBInputs{TT1,TT2,TR,TP,TQ,TW,Tsurf}
+    Ts::TT1 = 0.0u"°C"
     Tair::TT2 = 0.0u"°C"
     Lin::TR = 0.0u"W/m^2"
     Sin::TR = 0.0u"W/m^2"
diff --git a/src/Presets/Presets.jl b/src/models.jl
similarity index 86%
rename from src/Presets/Presets.jl
rename to src/models.jl
index 2978ec7e9cf19b78341bab86a69e68df7660dc65..ce5ebbcef783ebb1a05f06f9a96c275ec1c99135 100644
--- a/src/Presets/Presets.jl
+++ b/src/models.jl
@@ -1,22 +1,3 @@
-"""
-Pre-built CryoGrid configurations for rapid prototyping.
-"""
-module Presets
-
-using CryoGrid
-using CryoGrid.InputOutput: Resource
-using CryoGrid.Numerics
-using CryoGrid.Utils
-
-# physics modules
-using CryoGrid.Heat
-using CryoGrid.Hydrology
-using CryoGrid.Soils
-
-using Statistics
-
-export SoilHeatTile, SamoylovDefault
-
 """
     SoilHeatTile([heatop=:H], upperbc::BoundaryProcess, soilprofile::Profile, init::CryoGrid.Initializer...; grid::Grid=DefaultGrid_10cm, tile_kwargs...) where {F<:FreezeCurve}
 
@@ -61,5 +42,3 @@ const SamoylovDefault = (
         1000.0u"m" => 10.2u"°C",
     )
 )
-
-end
\ No newline at end of file
diff --git a/test/CryoGridSBIExt/runtests.jl b/test/CryoGridSBIExt/runtests.jl
index 5ddb3e0bbbc20c32b76f87b190f354a6a01131ee..aa2f09df67ee49c4c48760531a26ce590232975a 100644
--- a/test/CryoGridSBIExt/runtests.jl
+++ b/test/CryoGridSBIExt/runtests.jl
@@ -15,14 +15,15 @@ function cryogrid_test_problem(tspan)
     );
     heatop = Heat.Diffusion1D(:H)
     initT = initializer(:T, tempprofile)
-    tile = CryoGrid.Presets.SoilHeatTile(
+    tile = CryoGrid.SoilHeatTile(
         heatop,
         ## 10 W/m^2 in and out, i.e. net zero flux
         PeriodicBC(HeatBalance, CryoGrid.Dirichlet, 365.0u"d", 15.0u"°C", 0.0, -5.0u"°C"),
         ConstantBC(HeatBalance, CryoGrid.Neumann, -0.1u"W/m^2"),
         soilprofile,
+        inputs(),
         initT;
-        grid=CryoGrid.Presets.DefaultGrid_10cm,
+        grid=CryoGrid.DefaultGrid_10cm,
     );
     u0, _ = initialcondition!(tile, tspan);
     prob = CryoGridProblem(tile, u0, tspan)
diff --git a/test/IO/forcing_tests.jl b/test/IO/forcing_tests.jl
index 47be0699ea8110d9f182517f1757e0e4c3732c45..60f2dfdaf90ac19b85edfbbd1cd0abd86c734b22 100644
--- a/test/IO/forcing_tests.jl
+++ b/test/IO/forcing_tests.jl
@@ -2,23 +2,24 @@ using CryoGrid
 using Dates
 using Test, BenchmarkTools
 
-@testset "InterpolatedForcing" begin
+@testset "Interpolated1D" begin
     ts = DateTime(2020,1,1):Day(1):DateTime(2020,3,1)
     ys = randn(length(ts))
-    forcing = CryoGrid.InterpolatedForcing(ts, ys, :test)
-    @test all(forcing.interpolant.knots[1] .≈ convert_t.(ts))
-    @test forcing(Dates.datetime2epochms(ts[1])/1000.0) ≈ ys[1]
-    @test forcing(Dates.datetime2epochms(ts[end])/1000.0) ≈ ys[end]
-    @test_throws BoundsError forcing(-1.0)
-    @test_throws BoundsError forcing(Dates.datetime2epochms(ts[end])/1000.0+1.0)
+    arr = DimStack((test=DimArray(ys, (Ti(ts),)),))
+    interp = InputOutput.Interpolated1D(arr)
+    f = interp.test
+    @test f(Dates.datetime2epochms(ts[1])/1000.0) ≈ ys[1]
+    @test f(Dates.datetime2epochms(ts[end])/1000.0) ≈ ys[end]
+    @test_throws BoundsError f(-1.0)
+    @test_throws BoundsError f(Dates.datetime2epochms(ts[end])/1000.0+1.0)
     t1,y1 = ts[1], ys[1]
     t2,y2 = ts[2], ys[2]
-    @test forcing((Dates.datetime2epochms(t1) + Dates.datetime2epochms(t2))/2000.0) ≈ (y1+y2)/2
+    @test f((Dates.datetime2epochms(t1) + Dates.datetime2epochms(t2))/2000.0) ≈ (y1+y2)/2
     t = Dates.datetime2epochms(t1)/1000.0
-    benchres = @benchmark $forcing($t)
+    benchres = @benchmark $f($t)
     @test benchres.allocs == 0
     out = zeros(100)
     queries = t .+ (1:100);
-    benchres = @benchmark $out .= $forcing.($queries)
+    benchres = @benchmark $out .= $f.($queries)
     @test benchres.allocs == 0
 end
\ No newline at end of file
diff --git a/test/Physics/Heat/heat_conduction_tests.jl b/test/Physics/Heat/heat_conduction_tests.jl
index 692cf824dfe689a8d26e7b50402962f120be151c..78ed483181d7a2d09623112a0fcda6e0722c35df 100644
--- a/test/Physics/Heat/heat_conduction_tests.jl
+++ b/test/Physics/Heat/heat_conduction_tests.jl
@@ -79,19 +79,22 @@ using Test
 end
 @testset "Boundary conditions" begin
 	@testset "n-factors" begin
-		ts = DateTime(2010,1,1):Hour(1):DateTime(2010,1,1,4)
-		forcing = InterpolatedForcing(ts, [1.0,0.5,-0.5,-1.0,0.1]u"°C", :Tair)
-		tgrad = TemperatureBC(forcing, NFactor(nf=0.5, nt=1.0))
+		T_pos = TemperatureBC(10.0, NFactor(nf=0.5, nt=0.75))
+		T_neg = TemperatureBC(-10.0, NFactor(nf=0.5, nt=0.75))
 		heat = HeatBalance()
 		sub = TestGroundLayer(heat)
 		zerobc = ConstantBC(HeatBalance, CryoGrid.Dirichlet, 0.0u"°C")
 		function f1(t)
 			state = (T_ub=[Inf], nfactor=[Inf], t=t)
-			computediagnostic!(Top(zerobc), tgrad, state)
-			return boundaryvalue(tgrad, state)
+			computediagnostic!(Top(T_pos), T_pos, state)
+			T_ub = boundaryvalue(T_pos, state)
+			computediagnostic!(Top(T_neg), T_neg, state)
+			T_lb = boundaryvalue(T_neg, state)
+			return T_ub, T_lb
 		end
-		Tres = f1.(Dates.datetime2epochms.(ts)./1000.0)
-		@test all(Tres .≈ [1.0,0.5,-0.25,-0.5,0.1])
+		T1, T2 = f1(0.0)
+		@test T1 == 7.5
+		@test T2 == -5.0
 	end
 end
 @testset "Fourier solution" begin
diff --git a/test/Physics/Hydrology/water_balance_tests.jl b/test/Physics/Hydrology/water_balance_tests.jl
index ac819ccd68c33f4b23fc9c577d62644548124063..b4f0873064469310c2253384728f7208bd17bcd0 100644
--- a/test/Physics/Hydrology/water_balance_tests.jl
+++ b/test/Physics/Hydrology/water_balance_tests.jl
@@ -6,9 +6,9 @@ using CryoGrid.Utils
 using Test
 
 @testset "Bucket scheme" begin
-    testgrid = CryoGrid.Presets.DefaultGrid_2cm
+    testgrid = CryoGrid.DefaultGrid_2cm
     @testset "variables" begin
-        sub = TestGroundLayer(WaterBalance(BucketScheme()))
+        sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
         vars = CryoGrid.variables(sub)
         prognostic_vars = filter(CryoGrid.isprognostic, vars)
         diagnostic_vars = filter(CryoGrid.isdiagnostic, vars)
@@ -23,7 +23,7 @@ using Test
         @test isa(proc, WaterBalance)
     end
     @testset "computediagnostic!" begin
-        sub = TestGroundLayer(WaterBalance(BucketScheme()))
+        sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
         state = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
         # initialize fully saturated
         state.sat .= 1.0
@@ -33,7 +33,7 @@ using Test
         @test allfinite(state.θw)
     end
     @testset "computeprognostic!" begin
-        sub = TestGroundLayer(WaterBalance(BucketScheme()))
+        sub = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
         state = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
         # initialize fully saturated
         state.sat .= 1.0
@@ -53,8 +53,8 @@ using Test
         @test all(state.jw[2:end-1] .> zero(eltype(state.jw)))
     end
     @testset "interact!" begin
-        sub1 = TestGroundLayer(WaterBalance(BucketScheme()))
-        sub2 = TestGroundLayer(WaterBalance(BucketScheme()))
+        sub1 = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
+        sub2 = TestGroundLayer(WaterBalance(BucketScheme())) |> stripparams
         state1 = Diagnostics.build_dummy_state(testgrid[0.0u"m"..10.0u"m"], sub1, with_units=true)
         state2 = Diagnostics.build_dummy_state(testgrid[10.0u"m"..1000.0u"m"], sub2, with_units=true)
         # case 1: both saturated
diff --git a/test/Physics/Salt/runtests.jl b/test/Physics/Salt/runtests.jl
index d9dc484526d0bfc2b434dbe8a266f81a52cbd855..70f5154a42c50f92e3e5901d980c07dd764c6dd7 100644
--- a/test/Physics/Salt/runtests.jl
+++ b/test/Physics/Salt/runtests.jl
@@ -5,7 +5,7 @@ using CryoGrid.Utils
 using Test
 
 @testset "Salt" begin
-    testgrid = CryoGrid.Presets.DefaultGrid_2cm
+    testgrid = CryoGrid.DefaultGrid_2cm
     @testset "variables" begin
         soil = SalineGround()
         vars = CryoGrid.variables(soil)
diff --git a/test/Physics/Surface/seb_tests.jl b/test/Physics/Surface/seb_tests.jl
index c77626706b4246b3a4f4a2a961fbac46dd1446a5..d4d96beebdb5c9261795ba5a9a4e83848ada29b8 100644
--- a/test/Physics/Surface/seb_tests.jl
+++ b/test/Physics/Surface/seb_tests.jl
@@ -9,8 +9,8 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
     Sin = upreferred(200.0u"W/m^2")
     forcings = (; Tair, pr, qh, wind, Lin, Sin)
     @testset "Iterative" begin
-        seb = SurfaceEnergyBalance(
-            forcings;
+        seb = SurfaceEnergyBalance(;
+            forcings...,
             solscheme = Surface.Iterative(),
         )
         # TODO: would be good to have SEB work correctly with units to verify correctness;
@@ -29,8 +29,8 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
         @test seb_output.Qg > zero(seb_output.Qg)
     end
     @testset "Analytical" begin
-        seb = SurfaceEnergyBalance(
-            forcings;
+        seb = SurfaceEnergyBalance(;
+            forcings...,
             solscheme = Surface.Analytical(),
         )
         seb = pstrip(seb)
@@ -46,8 +46,8 @@ Surface.surfaceproperties(seb::SurfaceEnergyBalance, sub::TestGroundLayer) = seb
         @test seb_output.Qg > zero(seb_output.Qg)
     end
     @testset "Numerical" begin
-        seb = SurfaceEnergyBalance(
-            forcings;
+        seb = SurfaceEnergyBalance(;
+            forcings...,
             solscheme=Surface.Numerical(),
         )
         seb = pstrip(seb)
diff --git a/test/Physics/Surface/swb_tests.jl b/test/Physics/Surface/swb_tests.jl
index 6a52498477fa5069a69fb840d0cf33fca8ef900d..96cf4dbd64e2ba98b40dd8b198b6b0001960720b 100644
--- a/test/Physics/Surface/swb_tests.jl
+++ b/test/Physics/Surface/swb_tests.jl
@@ -1,13 +1,8 @@
 @testset "Surface Water Balance" begin
-    testgrid = CryoGrid.Presets.DefaultGrid_2cm
+    testgrid = CryoGrid.DefaultGrid_2cm
     @testset "interact!" begin
-        sub = TestGroundLayer(WaterBalance(BucketScheme()))
-        top = Top(
-            SurfaceWaterBalance(
-                ConstantForcing(1e-6u"m/s", :rainfall),
-                ConstantForcing(0.0u"m/s", :snowfall),
-            )
-        )
+        sub = stripparams(TestGroundLayer(WaterBalance(BucketScheme())))
+        top = Top(SurfaceWaterBalance(1e-6u"m/s", 0.0u"m/s"))
         stop = Diagnostics.build_dummy_state(testgrid, top, with_units=true)
         ssub = Diagnostics.build_dummy_state(testgrid, sub, with_units=true)
         # initialize fully saturated
diff --git a/test/Solvers/cglite_implicit_tests.jl b/test/Solvers/cglite_implicit_tests.jl
index bc91216c520c2d2ce0b19c672dab7edfba6809e7..5c437e86586f25f1c5f4b5030d71c0debac9aa7c 100644
--- a/test/Solvers/cglite_implicit_tests.jl
+++ b/test/Solvers/cglite_implicit_tests.jl
@@ -27,7 +27,7 @@ end
         α = upreferred(strat.soil.para.heat.kh_m) / upreferred(strat.soil.para.heat.ch_m)
         T_analytic = Heat.heat_conduction_linear_periodic_ub(T₀, A, P, ustrip(α))
         initT = initializer(:T, (layer, state) -> state.T .= T_analytic.(cells(state.grid), 0.0))
-        modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+        modelgrid = CryoGrid.DefaultGrid_2cm
         # modelgrid = Grid(z_top:0.02u"m":z_bot)
         # modelgrid = Grid(vcat(0.0:0.02:1.0, 1.05:0.05:5.0, 5.1:0.1:10.0)*u"m")
         tile = Tile(strat, modelgrid, initT)
@@ -58,7 +58,7 @@ end
             z_bot => Bottom(ConstantFlux(HeatBalance, 0.0))
         );
         initT = initializer(:T, -1.0)
-        modelgrid = CryoGrid.Presets.DefaultGrid_2cm
+        modelgrid = CryoGrid.DefaultGrid_2cm
         tile = Tile(strat, modelgrid, initT)
         # define time span, 5 years
         tspan = (0.0, 5*365*24*3600.0)
@@ -68,7 +68,6 @@ end
         sol = solve(prob, LiteImplicitEuler(), dt=24*3600)
         out = CryoGridOutput(sol)
 
-        kh_w, kh_i, kh_a, kh_m, kh_o = Heat.thermalconductivities(strat.soil)
         θm = mineral(strat.soil)
         θo = organic(strat.soil)
         θ_s = (θw=0.0, θi=porosity(strat.soil), θa=0.0, θm=θm, θo=θo)
diff --git a/test/benchmarks.jl b/test/benchmarks.jl
index 87c01f4f0c7a6146d96635ec8a73988dcb0077d3..19e67e87c353b50b959b296f0fa067d7c9410b3a 100644
--- a/test/benchmarks.jl
+++ b/test/benchmarks.jl
@@ -15,7 +15,7 @@ function run_solver_benchmark(solvers, grids, dts, tspan, upperbc;
         # high accuracy reference solution, 1 minute time steps
         println("Computing reference solution for grid spacing $min_dx with forward Euler @ $dt_ref s time steps ...")
         # basic 1-layer heat conduction model (defaults to free water freezing scheme)
-        model = Presets.SoilHeat(:H, upperbc, profile; grid=grid, freezecurve=freezecurve)
+        model = SoilHeat(:H, upperbc, profile; grid=grid, freezecurve=freezecurve)
         p = copy(model.pproto)
         p .= params
         prob = CryoGridProblem(model,tspan,p)
@@ -55,7 +55,7 @@ forcings = loadforcings(filename);
 # use air temperature as upper boundary forcing
 tair = TemperatureBC(forcings.Tair)
 solvers = [Euler, DP5, ROCK2, ROCK4, Trapezoid, ROS3P]
-grids = [Presets.DefaultGrid_2cm, Presets.DefaultGrid_5cm, Presets.DefaultGrid_10cm, Presets.DefaultGrid_20cm]
+grids = [DefaultGrid_2cm, DefaultGrid_5cm, DefaultGrid_10cm, DefaultGrid_20cm]
 dts = [2*60.0, 10*60.0, 30*60.0, 3600.0]
 tspan = (DateTime(2010,10,30),DateTime(2011,10,30))
 results_freeW_tair = run_solver_benchmark(solvers, grids, dts, tspan, tair)
diff --git a/test/test_problems.jl b/test/test_problems.jl
index d12036f201ce6d9fe4bbae0f764efd2c0e35d547..93c81569e42a8b93acceaa2e3f8601472cf7ad07 100644
--- a/test/test_problems.jl
+++ b/test/test_problems.jl
@@ -13,9 +13,9 @@ function test_heat_conduction_freeW_periodic_bc()
     initT = initializer(:T, tempprofile)
     # Define the simulation time span.
     tspan = (0.0,2*24*3600.0)
-    discretization = CryoGrid.Presets.DefaultGrid_10cm
+    discretization = CryoGrid.DefaultGrid_10cm
     Qbot = -0.01u"W/m^2"
-    tile = CryoGrid.Presets.SoilHeatTile(
+    tile = CryoGrid.SoilHeatTile(
         heatop,
         PeriodicBC(HeatBalance, CryoGrid.Neumann, 24*3600.0, 20.0, 0.0, -5.0),
         ConstantBC(HeatBalance, CryoGrid.Neumann, Qbot),
@@ -32,7 +32,7 @@ function test_water_flow_bucket_scheme(topflux=0.0u"m/s")
     initsat = initializer(:sat, 0.5)
     # Define the simulation time span.
     tspan = (0.0, 24*3600.0)
-    discretization = CryoGrid.Presets.DefaultGrid_10cm
+    discretization = CryoGrid.DefaultGrid_10cm
     strat = Stratigraphy(
         0.0u"m" => Top(ConstantBC(WaterBalance, CryoGrid.Neumann, topflux),),
         0.0u"m" => Ground(SimpleSoil(por=0.5, sat=0.5), heat=nothing, water=WaterBalance(BucketScheme())),