Repository: Deltares/Wflow.jl (language: python)

Issue:
Remove vertical concepts `HBV` and `FLEXTopo`
### Feature type

Removing existing functionality

### Improvement Description

See Discussion https://github.com/Deltares/Wflow.jl/discussions/342.

### Implementation Description

_No response_

### Additional Context

_No response_


Discussion / hints:

Fix patch:
diff --git a/docs/make.jl b/docs/make.jl
index a710b5fb3..22dd23a5d 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -23,8 +23,6 @@ pages = [
         "model_docs/model_configurations.md",
         "Vertical concepts" => [
             "model_docs/vertical/sbm.md",
-            "model_docs/vertical/hbv.md",
-            "model_docs/vertical/flextopo.md",
             "model_docs/vertical/sediment.md",
             "model_docs/shared_concepts.md",
         ],
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index a0bd2349c..6133e9600 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## v1.0.0
+
+### Fixed
+
+### Changed
+- Removed vertical concepts `HBV` and `FLEXTopo`.
+
+### Added
+
 ## [unreleased]
 
 ### Fixed
diff --git a/docs/src/images/flextopo_julia_1class.png b/docs/src/images/flextopo_julia_1class.png
deleted file mode 100644
index cb662a324..000000000
Binary files a/docs/src/images/flextopo_julia_1class.png and /dev/null differ
diff --git a/docs/src/images/flextopo_julia_3class.png b/docs/src/images/flextopo_julia_3class.png
deleted file mode 100644
index dd8175800..000000000
Binary files a/docs/src/images/flextopo_julia_3class.png and /dev/null differ
diff --git a/docs/src/images/hbv-soilmoist.png b/docs/src/images/hbv-soilmoist.png
deleted file mode 100644
index 1e6e3ed7d..000000000
Binary files a/docs/src/images/hbv-soilmoist.png and /dev/null differ
diff --git a/docs/src/images/hbv-upper.png b/docs/src/images/hbv-upper.png
deleted file mode 100644
index 4aae2eb12..000000000
Binary files a/docs/src/images/hbv-upper.png and /dev/null differ
diff --git a/docs/src/images/hbv96.png b/docs/src/images/hbv96.png
deleted file mode 100644
index aa95d9f6e..000000000
Binary files a/docs/src/images/hbv96.png and /dev/null differ
diff --git a/docs/src/model_docs/model_configurations.md b/docs/src/model_docs/model_configurations.md
index f2872153c..6baf6a2ec 100644
--- a/docs/src/model_docs/model_configurations.md
+++ b/docs/src/model_docs/model_configurations.md
@@ -123,128 +123,6 @@ lateral.land => struct ShallowWaterLand{T}
 The local inertial approach is described in more detail in the section [Local inertial
 model](@ref local_inertial).
 
-## [wflow\_hbv](@id config_hbv)
-The Hydrologiska Byrans Vattenbalansavdelning (HBV) model was introduced back in 1972 by the
-Swedish Meteological and Hydrological Institute (SMHI).  The HBV model is mainly used for
-runoff simulation and hydrological forecasting. The model is particularly useful for
-catchments where snow fall and snow melt are dominant factors, but application of the model
-is by no means restricted to these type of catchments.
-
-The model is based on the HBV-96 model. However, the hydrological routing represented in HBV
-by a triangular function controlled by the MAXBAS parameter has been removed. Instead, the
-kinematic wave function is used to route the water downstream. All runoff that is generated
-in a cell in one of the HBV reservoirs is added to the kinematic wave reservoir at the end
-of a timestep. There is no connection between the different HBV cells within the model.
-
-A catchment is divided into a number of grid cells. For each of the cells individually,
-daily runoff is computed through application of the HBV-96 of the HBV model. The use of the
-grid cells offers the possibility to turn the HBV modelling concept, which is originally
-lumped, into a distributed model.
-
-![wflow_hbv model](../images/hbv96.png)
-
-The figure above shows a schematic view of hydrological response simulation with the
-HBV-modelling concept. The land-phase of the hydrological cycle is represented by three
-different components: a snow routine, a soil routine and a runoff response routine. Each
-component is discussed in more detail below.
-
-The vertical HBV concept is described in section [HBV vertical concept](@ref vert_hbv). The
-routing for river and overland flow is described in the section [Kinematic wave](@ref
-kin_wave).
-
-Below the mapping for wflow\_hbv (type `hbv`) to the vertical HBV concept (instance of
-`struct HBV`) and the different lateral concepts. For an explanation about the type
-parameters between curly braces after the `struct` name see the section on model parameters.
-
-```julia
-vertical => struct HBV{T}
-lateral.subsurface => struct LateralSSF{T}
-lateral.land => struct SurfaceFlow{T,R,L}
-lateral.river => struct SurfaceFlow{T,R,L}
-lateral.river.lake => struct NaturalLake{T} # optional
-lateral.river.reservoir => struct SimpleReservoir{T} # optional
-```
-
-## [wflow\_flextopo](@id config_flextopo)
-The FLEXTopo model is a process-based model, which consists of different parallel classes
-connected through their groundwater storage. These classes are usually delineated from
-topographical data to represent the variability in hydrological processes across
-user-defined Hydrological Response Units (HRU). The main assumption underlying the concept,
-which was first introduced by Savenije (2010), is that different parts of the landscape
-fulfill different tasks in runoff generation and, hence, can be represented by different
-model structures. The strength of the concept is that the definition of classes and
-associated model structures is modular and flexible and not fixed to a predefined model
-structure. The flexible approach allows to develop process-based models for different
-topographic, climatic, geologic and land use conditions, making use of the available data
-and expert knowledge.
-
-The kinematic wave function is used to route the water downstream. In a similar way as for
-HBV, all runoff that is generated in a cell in one of the FLEXTopo storages is added to the
-kinematic wave reservoir at the end of a timestep. There is no connection between the
-different vertical FLEXTopo cells within the model. The FLEXTopo model is implemented in a
-fully distributed way in the wflow Julia framework.
-
-In wflow\_flextopo, the user is free to determine the number of classes and which model
-components to include or exclude for each class, this is done in the TOML file. Currently,
-for each storage, it is possible to bypass the storage and pass on the fluxes to the next
-model component. Interested users can contribute to the code by adding other
-conceptualizations for each storage components.
-
-```julia
-[model]
-type = "flextopo"
-classes = ["h", "p", "w"]     #user can set the number and name of each class.
-
-# for each component which is class specific, the user can select which conceptualization
-# to apply for each class as defined above in classes = ["h", "p", "w"]
-select_snow = ["common_snow_hbv"]
-# available options are ["common_snow_hbv", "common_snow_no_storage"]
-select_interception = ["interception_overflow", "interception_overflow", "interception_overflow"]
-# available options are ["interception_overflow", "interception_no_storage"]
-select_hortonponding = ["hortonponding_no_storage", "hortonponding_no_storage", "hortonponding_no_storage"]
-# available options are ["hortonponding", "hortonponding_no_storage"]
-select_hortonrunoff = ["hortonrunoff_no_storage", "hortonrunoff_no_storage", "hortonrunoff_no_storage"]
-# available options are ["hortonrunoff", "hortonrunoff_no_storage"]
-select_rootzone = ["rootzone_storage", "rootzone_storage", "rootzone_storage"]
-# available options are ["rootzone_storage", "rootzone_no_storage"]
-select_fast = ["fast_storage", "fast_storage", "fast_storage"]
-# available options are ["fast_storage", "fast_no_storage"]
-select_slow = ["common_slow_storage"]
-# available options are ["common_slow_storage", "slow_no_storage"]
-```
-
-A schematic representation of the most complete model structure including all storage
-components, as currently implemented in the code, is shown in the Figure below. When setting
-up the model with multiple classes, model structures can be adapted by bypassing storages or
-turning parameter values on or off (e.g.: percolation or capillary rise, non-linear versus
-linear outflow of the fast runoff etc.), an example of a three class model is shown in
-[FLEXTopo vertical concept](@ref vert_flextopo).
-
-![flextopo_julia_1class.png](../images/flextopo_julia_1class.png)
-*Schematic representation of the FLEXTopo model for a single class model including all
-storages and fluxes. Main parameters are denoted in red.*
-
-In the staticmaps, the user needs to provide maps of the fraction of each class within each
-cell, as shown below with `hrufrac`. For each model parameter which is class specific, an
-extra dimension `classes` is required in the staticmaps netcdf. For an example model, see
-[FLEXTopo example model](@ref wflow_flextopo_data).
-
-```julia
-[input.vertical]
-hrufrac = "hrufrac_lu"
-```
-
-Parameter multiplication of model parameters which are defined for several classes is
-possible through the TOML file:
-
-```julia
-[input.vertical.kf]
-netcdf.variable.name = "kf"
-scale = [1.0, 3.0, 4.0]
-offset = [0.0, 0.0, 0.0]
-class = ["h", "p", "w"]
-```
-
 ## [wflow\_sediment](@id config_sediment)
 The processes and fate of many particles and pollutants impacting water quality at the
 catchment level are intricately linked to the processes governing sediment dynamics. Both
@@ -284,8 +162,7 @@ required dynamic inputs to run wflow\_sediment are:
 -  River water level in the kinematic wave,
 -  Rainfall interception by the vegetation.
 
-These inputs can be obtained from other wflow models such as wflow\_sbm, wflow\_hbv or from
-other sources.
+These inputs can be obtained from wflow\_sbm or from other sources.
 
 Model outputs can be saved for both the inland and the instream part of the model. Some
 examples are listed below.
@@ -323,7 +200,4 @@ Bedconc = "Bedconc"
 ## References
 + Köhler, L., Mulligan, M., Schellekens, J., Schmid, S., Tobón, C., 2006, Hydrological
   impacts of converting tropical montane cloud forest to pasture, with initial reference to
-  northern Costa Rica. Final Technical Report DFID‐FRP Project No. R799.
-
-+ Savenije, H. H. G. (2010). HESS opinions “topography driven conceptual modelling (FLEX-Topo).”
-  Hydrology and Earth System Sciences, 14(12), 2681–2692. https://doi.org/10.5194/hess-14-2681-2010
+  northern Costa Rica. Final Technical Report DFID‐FRP Project No. R799.
\ No newline at end of file
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 02477283c..1edea2f83 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -135,164 +135,6 @@ profile `kv` is used and `z_layered` is required as input.
 | `waterlevel_river` | water level river | mm | - |
 | `total_storage` | total water storage (excluding floodplains, lakes and reservoirs) | mm | - |
 
-
-## [HBV](@id params_hbv)
-The Table below shows the parameters (fields) of struct `HBV`, including a description of
-these parameters, the unit, and default value if applicable. The parameters in bold
-represent model parameters that can be set through static and forcing input data (netCDF),
-and can be listed in the TOML configuration file under `[input.vertical]`, to map the
-internal model parameter to the external netCDF variable.
-
-|  parameter | description    | unit | default |
-|:---------------| --------------- | ---------------------- | ----- |
-| **`cfmax`**  | degree-day factor | mm ᵒC``^{-1}`` Δt``^{-1}`` | 3.75653 mm ᵒC``^{-1}`` day``^{-1}``  |
-| **`tt`**  | threshold temperature for snowfall| ᵒC | -1.41934  |
-| **`tti`**  | threshold temperature interval length | ᵒC | 1.0  |
-| **`ttm`**  | threshold temperature for snowmelt  | ᵒC | -1.41934  |
-| **`whc`**  | water holding capacity as fraction of current snow pack  | - | 0.1  |
-| **`g_tt`**  | threshold temperature for snowfall above glacier  | ᵒC| 0.0  |
-| **`g_cfmax`**  | Degree-day factor for glacier  | mm ᵒC``^{-1}`` Δt``^{-1}``| 3.0 mm ᵒC``^{-1}`` day``^{-1}`` |
-| **`g_sifrac`**  | fraction of the snowpack on top of the glacier converted into ice  | Δt``^{-1}``| 0.001 day``^{-1}`` |
-| **`glacierfrac`**  | fraction covered by a glacier | - | 0.0  |
-| **`glacierstore`**  | water within the glacier  | mm | 5500.0  |
-| **`fc`**  | field capacity  | mm | 260.0  |
-| **`betaseepage`**  | exponent in soil runoff generation equation  | - | 1.8  |
-| **`lp`**  | fraction of field capacity below which actual evaporation=potential evaporation | - | 0.53 |
-| **`k4`**  | recession constant baseflow | Δt``^-1`` | 0.02307 day``^{-1}`` |
-| **`kquickflow`**  | recession constant upper reservoir | Δt``^-1`` |  0.09880 day``^{-1}`` |
-| **`suz`**  | Level over which `k0` is used | mm |  100.0 |
-| **`k0`**  | recession constant upper reservoir | Δt``^-1`` | 0.30 day``^{-1}`` |
-| **`khq`**  | recession rate at flow `hq`| Δt``^-1`` | 0.09880 day``^{-1}`` |
-| **`hq`**  | high flow rate hq for which recession rate of upper reservoir is known | mm Δt``^-1`` | 3.27 mm day``^{-1}`` |
-| **`alphanl`**  | measure of non-linearity of upper reservoir | - | 1.1 |
-| **`perc`**  | percolation from upper to lower zone |  mm Δt``^-1`` | 0.4 mm day``^{-1}`` |
-| **`cfr`**  | refreezing efficiency constant in refreezing of freewater in snow | - | 0.05  |
-| **`pcorr`**  | correction factor for precipitation  | - | 1.0  |
-| **`rfcf`**  | correction factor for rainfall   | - | 1.0  |
-| **`sfcf`**  | correction factor for snowfall | - | 1.0  |
-| **`cflux`**  | maximum capillary rise from runoff response routine to soil moisture routine | mm Δt``^-1`` | 2.0 mm day``^{-1}``  |
-| **`icf`**  | maximum interception storage (in forested and non-forested areas) | mm | 2.0  |
-| **`cevpf`**  | correction factor for potential evaporation | - | 1.0  |
-| **`epf`**  | exponent of correction factor for evaporation on days with precipitation | mm``^{-1}`` | 1.0  |
-| **`ecorr`**  | evaporation correction | - | 1.0  |
-| **`precipitation`**           | precipitation     | mm Δt``^-1`` | - |
-| **`temperature`**             | temperature    | ᵒC | - |
-| **`potential_evaporation`**            | potential evapotranspiration     | mm Δt``^-1`` | - |
-| `potsoilevap`             | potential soil evaporation      | mm Δt``^-1`` | - |
-| `soilevap`             | soil evaporation | mm Δt``^-1`` | - |
-| `intevap`             | evaporation from interception storage     | mm Δt``^-1`` | - |
-| `actevap`             | actual evapotranspiration (intevap + soilevap) | mm Δt``^-1`` | - |
-| `interceptionstorage`             | actual interception storage     | mm | - |
-| `snowwater`             | available free water in snow    | mm | - |
-| `snow`             | snow pack    | mm | - |
-| `rainfallplusmelt`    | snow melt + precipitation as rainfall    | mm Δt``^-1`` | - |
-| `soilmoisture`        | actual soil moisture | mm | - |
-| `directrunoff`        | direct runoff to upper zone    | mm Δt``^-1`` | - |
-| `hbv_seepage`         | recharge to upper zone  | mm Δt``^-1``  | - |
-| `in_upperzone`         | water inflow into upper zone  | mm Δt``^-1``  | - |
-| `upperzonestorage`     | water content of the upper zone  | mm  | - |
-| `quickflow`         |  specific runoff (quickflow part)  | mm Δt``^-1``  | - |
-| `real_quickflow`         | specific runoff (quickflow), if K upper zone is precalculated  | mm Δt``^-1``  | - |
-| `percolation`         | actual percolation to the lower zone | mm Δt``^-1``  | - |
-| `capflux`         | capillary rise | mm Δt``^-1``  | - |
-| `lowerzonestorage`         | water content of the lower zone | mm  | - |
-| `baseflow`         | specific runoff (baseflow part) per cell  | mm Δt``^-1``  | - |
-| `runoff`         | total specific runoff per cell  | mm Δt``^-1``  | - |
-
-
-## [FLEXtopo](@id params_flextopo)
-The Table below shows the parameters (fields) of struct `FLEXTOPO`, including a description
-of these parameters, the unit, and default value if applicable. The parameters in bold
-represent model parameters that can be set through static and forcing input data (netCDF),
-and can be listed in the TOML configuration file under `[input.vertical]`, to map the
-internal model parameter to the external netCDF variable.
-
-|  parameter | description    | unit | default |
-|:---------------| --------------- | ---------------------- | ----- |
-| **`cfmax`**  | degree-day factor | mm ᵒC``^{-1}`` Δt``^{-1}`` | 3.75653 mm ᵒC``^{-1}`` day``^{-1}``  |
-| **`tt`**  | threshold temperature for snowfall| ᵒC | -1.41934  |
-| **`tti`**  | threshold temperature interval length | ᵒC | 1.0  |
-| **`ttm`**  | threshold temperature for snowmelt  | ᵒC | -1.41934  |
-| **`whc`**  | water holding capacity as fraction of current snow pack  | - | 0.1  |
-| **`cfr`**  | refreezing efficiency constant in refreezing of freewater in snow | - | 0.05  |
-| **`g_tt`**  | threshold temperature for snowfall above glacier  | ᵒC| 0.0  |
-| **`g_cfmax`**  | Degree-day factor for glacier  | mm ᵒC``^{-1}`` Δt``^{-1}``| 3.0 mm ᵒC``^{-1}`` day``^{-1}`` |
-| **`g_sifrac`**  | fraction of the snowpack on top of the glacier converted into ice  | Δt``^{-1}`` | 0.001 day``^{-1}`` |
-| **`glacierfrac`**  | fraction covered by a glacier | - | 0.0  |
-| **`glacierstore`**  | water within the glacier  | mm | 5500.0  |
-| **`ecorr`**  | evaporation correction | - | 1.0  |
-| **`pcorr`**  | correction factor for precipitation  | - | 1.0  |
-| **`rfcf`**  | correction factor for rainfall   | - | 1.0  |
-| **`sfcf`**  | correction factor for snowfall | - | 1.0  |
-| **`imax`**  | maximum interception storage (``I_\mathrm{max}``) | mm | 3.0  |
-| **`shmax`**  | maximum horton ponding storage capacity (``S_\mathrm{Hmax}``)  | mm | 30.0  |
-| **`srmax`**  | maximum root zone storage capacity (``S_\mathrm{Rmax}``)  | mm | 260.0  |
-| **`beta`**  | exponent in soil runoff generation equation  | - | 0.3  |
-| **`lp`**  | fraction of root zone capacity below which actual evaporation=potential evaporation (``L_\mathrm{P}``) | - | 0.3 |
-| **`ks`**  | recession constant slow groundwater storage (``K_\mathrm{S}``) | Δt``^-1`` | 0.006 day``^{-1}`` |
-| **`kf`**  | recession constant fast storage (``K_\mathrm{F}``) | Δt``^-1`` |  0.1 day``^{-1}`` |
-| **`khf`**  | recession constant horton runoff storage (``K_\mathrm{Hf}``) | Δt``^-1`` |  0.5 day``^{-1}`` |
-| **`alfa`**  | measure of non-linearity of upper reservoir (``\alpha``) | - | 1.0 |
-| **`perc`**  | maximum percolation flux from root zone to slow storage (``Q_\mathrm{perc,max}``) |  mm Δt``^-1``| 0.30 mm day``^{-1}`` |
-| **`cap`**  | maximum capillary rise from slow storage to root zone (``Q_\mathrm{cap,max}``)| mm Δt``^-1`` | 0.20 mm day``^{-1}``  |
-| **`ds`**  | splitter parameter determining fraction of root zone outflow to slow storage (``d_\mathrm{s}``) | - | 0.2 |
-| **`shmin`**  | minimum storage capacity in horton ponding (relative to ``S_\mathrm{Hmax}``) (``S_\mathrm{Hmin}``) | [-] |  0.2 |
-| **`facc0`**  | maximum modelled accumulated frost resulting in `shmin` (``F_\mathrm{acc,fr0}``) | [ᵒC Δt] |  -3.0 |
-| **`facc1`**  | minimum modelled accumulated frost resulting in `shmin` (``F_\mathrm{acc,fr1}``) | [ᵒC Δt] |  0.0 |
-| **`fdec`**  | exponent for the decline of infiltration capacity (``F_\mathrm{dec}``) | [-] |  0.2 |
-| **`fmax`**  | maximum infiltration capacity from horton ponding (``F_\mathrm{max}``) | [mm Δt``^-1``] |  2.0 |
-| **`kmf`**  | melt coefficient of frozen topsoil (``K_\mathrm{mf}``) | [-] |  1.0 |
-| **`hrufrac`**  | fraction of class within cell (``F_\mathrm{hrufrac}``)| - | 1/length(classes) |
-| **`precipitation`**           | precipitation     | mm Δt``^-1`` | - |
-| **`temperature`**             | temperature    | ᵒC | - |
-| **`potential_evaporation`**            | potential evapotranspiration     | mm Δt``^-1`` | - |
-| `precipcorr`           | corrected precipitation     | mm Δt``^-1`` | - |
-| `epotcorr`           | corrected potential evaporation     | mm Δt``^-1`` | - |
-| `snow`             | snow water (``S_\mathrm{W}``)   | mm | - |
-| `snowwater`             | available free water in snow    | mm | - |
-| `interceptionstorage`             | interception storage (``S_\mathrm{I}``)    | mm | - |
-| `interceptionstorage_m`             | average interception storage over classes (``S_\mathrm{I}``)    | mm | - |
-| `hortonpondingstorage`             | horton ponding storage (``S_\mathrm{Hp}``)    | mm | - |
-| `hortonpondingstorage_m`             | average horton ponding storage over classes (``S_\mathrm{Hp}``)    | mm | - |
-| `hortonrunoffstorage`             | horton runoff storage (``S_\mathrm{Hf}``)    | mm | - |
-| `hortonrunoffstorage_m`             | average horton runoff storage over classes (``S_\mathrm{Hf}``)    | mm | - |
-| `rootzonestorage`             | root zone storage (``S_\mathrm{R}``)  | mm | - |
-| `rootzonestorage_m`             | average root zone storage over classes (``S_\mathrm{R}``)  | mm | - |
-| `faststorage`             | fast storage (``S_\mathrm{F}``)  | mm | - |
-| `faststorage_m`             | average fast storage over classes (``S_\mathrm{F}``)  | mm | - |
-| `slowstorage`             | slow storage (``S_\mathrm{S}``)  | mm | - |
-| `potsoilevap`             | potential soil evaporation (``E_\mathrm{P}``)     | mm Δt``^-1`` | - |
-| `soilevap`             | soil evaporation | mm Δt``^-1`` | - |
-| `intevap`             | evaporation from interception storage (``E_\mathrm{I}``)    | mm Δt``^-1`` | - |
-| `intevap_m`             | average evaporation from interception storage over classes  (``E_\mathrm{I}``)   | mm Δt``^-1`` | - |
-| `hortonevap`             | evaporation from horton ponding storage (``E_\mathrm{H}``)    | mm Δt``^-1`` | - |
-| `hortonevap_m`             | average evaporation from horton ponding storage over classes  (``E_\mathrm{H}``)   | mm Δt``^-1`` | - |
-| `rootevap`             | evaporation from root zone storage (``E_\mathrm{R}``)    | mm Δt``^-1`` | - |
-| `rootevap_m`             | average evaporation from root zone storage over classes  (``E_\mathrm{R}``)   | mm Δt``^-1`` | - |
-| `actevap`             | actual evapotranspiration (intevap + hortonevap + rootevap) (``E_\mathrm{A}``) | mm Δt``^-1`` | - |
-| `actevap_m`             | average actual evapotranspiration (intevap + hortonevap + rootevap) over classes (``E_\mathrm{A}``) | mm Δt``^-1`` | - |
-| `precipeffective`    | Effective precipitation (``P_\mathrm{E}``)   | mm Δt``^-1`` | - |
-| `rainfallplusmelt`    | snow melt + precipitation as rainfall (``P_\mathrm{M} + P_\mathrm{R}``)   | mm Δt``^-1`` | - |
-| `snowmelt`    | snowfall   | mm Δt``^-1`` | - |
-| `snowfall`    | snowfall   | mm Δt``^-1`` | - |
-| `facc`    | modeled accumulated frost   | ᵒC Δt | - |
-| `qhortonpond`        | Flux from the hortonian ponding storage to the hortonian runoff storage (``Q_\mathrm{H}``)   | mm Δt``^-1`` | - |
-| `qhortonrootzone`        | Flux from the hortonian ponding storage to the root zone storage (``Q_\mathrm{HR}``)   | mm Δt``^-1`` | - |
-| `qhortonrun`        | Flux from the hortonian runoff storage (``Q_\mathrm{Hf}``)   | mm Δt``^-1`` | - |
-| `qrootzone`        | Flux from the root zone storage (``Q_\mathrm{R}``)  | mm Δt``^-1`` | - |
-| `qrootzonefast`        | Pref. recharge to fast storage (``Q_\mathrm{RF}``)  | mm Δt``^-1`` | - |
-| `qrootzoneslow_m`        | Pref. recharge to slow storage sum classes (``Q_\mathrm{RS}``)   | mm Δt``^-1`` | - |
-| `qcapillary`        | Capillary flux from the slow to the root-zone storage (``Q_\mathrm{cap}``)  | mm Δt``^-1`` | - |
-| `qcapillary_m`        | Capillary flux from the slow to the root-zone storage sum classes (``Q_\mathrm{cap}``)   | mm Δt``^-1`` | - |
-| `qpercolation`        | Percolation flux from the root-zone to the slow storage (``Q_\mathrm{perc}``)  | mm Δt``^-1`` | - |
-| `qpercolation_m`        | Percolation flux from the root-zone to the slow storage sum classes (``Q_\mathrm{perc}``)   | mm Δt``^-1`` | - |
-| `qfast`        | runoff from fast storage (``Q_\mathrm{F}``)  | mm Δt``^-1`` | - |
-| `qfast_tot`        | sum of fast runoff (from fast and horton runoff storages) over classes  | mm Δt``^-1`` | - |
-| `qslow`        | runoff from slow storage (``Q_\mathrm{S}``)  | mm Δt``^-1`` | - |
-| `runoff`         | total specific runoff per cell (qslow + qfast_tot) (``Q``) | mm Δt``^-1``  | - |
-| `wb_tot`             | total water balance | mm Δt``^-1`` | - |
-
-
 ## [Sediment](@id params_sediment)
 
 The Table below shows external parameters that can be set through static input data
diff --git a/docs/src/model_docs/shared_concepts.md b/docs/src/model_docs/shared_concepts.md
index ebda6144b..cbb8f1940 100644
--- a/docs/src/model_docs/shared_concepts.md
+++ b/docs/src/model_docs/shared_concepts.md
@@ -48,10 +48,10 @@ exceeds `whc`, either through snow melt or incoming rainfall, the surplus water
 
 ### Glacier modelling
 
-Glacier processes can be modelled if the snow model is enabled. For the vertical HBV concept
-snow modelling is not optional. Glacier modelling is very close to snow modelling and
-considers two main processes: glacier build-up from snow turning into firn/ice (using the
-HBV-light model) and glacier melt (using a temperature degree-day model).
+Glacier processes can be modelled if the snow model is enabled. Glacier modelling is very
+similar to snow modelling and considers two main processes: glacier build-up from snow turning
+into firn/ice (using the HBV-light model) and glacier melt (using a temperature degree-day
+model).
 
 The definition of glacier boundaries and initial volume is defined in three parameters.
 `glacierfrac` is a parameter that gives the fraction of each grid cell covered by a glacier
diff --git a/docs/src/model_docs/structures.md b/docs/src/model_docs/structures.md
index 78ca20a85..5f25e6f3f 100644
--- a/docs/src/model_docs/structures.md
+++ b/docs/src/model_docs/structures.md
@@ -21,7 +21,7 @@ end
 The `lateral` field of the `struct Model` can contain different lateral concepts. For each
 wflow model these different lateral concepts are mapped through the use of a `NamedTuple`.
 The `vertical` field of the `struct Model` always contains one vertical concept, for example
-the SBM or HBV vertical concept.
+the SBM vertical concept.
 
 Below an example how lateral concepts are mapped for the SBM model through a `NamedTuple`:
 
diff --git a/docs/src/model_docs/vertical/flextopo.md b/docs/src/model_docs/vertical/flextopo.md
deleted file mode 100644
index 114b54920..000000000
--- a/docs/src/model_docs/vertical/flextopo.md
+++ /dev/null
@@ -1,357 +0,0 @@
-# [FLEXTopo](@id vert_flextopo)
-
-## Introduction
-This section describes the different vertical processes available as part of the vertical
-FLEXTopo concept. This concept is part of the wflow\_flextopo model. The FLEXTopo model is a
-process-based model, which consists of different parallel classes connected through their
-groundwater storage. These classes are usually delineated from topographical data to
-represent the variability in hydrological processes across user-defined Hydrological
-Response Units (HRU). The main assumption underlying the concept, which was first introduced
-by Savenije (2010), is that different parts of the landscape fulfill different tasks in
-runoff generation and, hence, can be represented by different model structures. Commonly
-used classes include hillslopes, plateau and wetlands. Hillslopes are steep areas in a
-catchment and are generally forested. The dominant runoff process in hillslopes is assumed
-to be characterized by subsurface flow. Plateaus are defined as relatively flat and are
-relatively high above the stream, with deep groundwater levels. Depending on the specific
-conditions, the dominant runoff processes are groundwater recharge, quick subsurface flow
-and hortonian overland flow, which is especially important in agricultural areas. Saturation
-overland flow and capillary rise are the dominant processes on the riparian wetland class,
-where groundwater levels are shallow and assumed to rise quickly during an event. The
-strength of the concept is that the definition of classes and associated model structures is
-modular and flexible and not constrained by a predefined fixed model structure. The flexible
-approach allows to develop process-based models for different topographic, climatic,
-geologic and land use conditions, making use of the available data and expert knowledge. The
-FLEXTopo modeling approach has been applied in a lumped and distributed way in various
-applications across the globe (Gao et al., 2014; Euser et al., 2015; Hanus et al., 2021;
-Hrachowitz et al. 2021; Hulsman et al. 2021; Bouaziz et al., 2022).
-
-The wflow\_flextopo model is set-up in a modular way, implying that the user is free to
-determine the number of classes and which processes and/or parameters to include or exclude
-for each class. For each cell in the model domain, the percentage of each class in the cell
-is provided by the user in the staticmaps. The most complete model structure implemented in
-wflow\_flextopo is shown in the figure below. However, it is also possible to bypass each
-bucket or to deactivate processes through parameter values (see [FLEXTopo
-configuration](@ref config_flextopo)). It is also possible for users to contribute to the
-code by adding different conceptualizations for the different storages. When defining
-several classes, the model structures for each class are implemented in parallel, except for
-the common glacier, snow and groundwater processes, which are not class specific. An example
-of a three classes model is also shown in the Figure below.
-
-![flextopo_julia_1class.png](../../images/flextopo_julia_1class.png)
-*Schematic representation of the FLEXTopo model for a single class model including all
-storages and fluxes. Main parameters are denoted in red.*
-
-![flextopo_julia_3class.png](../../images/flextopo_julia_3class.png)
-*Example of a three class model with different model structure configurations per class*
-
-The descriptions below of each of the FLEXTopo model components are given for the most
-complete model structure with the symbols as shown in the schematic representation of the
-one class model in the Figure above. By bypassing storages and/or setting parameters to
-specific values, the model structure can be adapted to a user-defined model structure.
-
-## Snow
-The snow model is described in [Snow and glaciers](@ref snow_and_glac).
-
-## Glaciers
-Glacier processes are described in [Snow and glaciers](@ref snow_and_glac). Glacier
-modelling is enabled by specifying the following in the TOML file:
-
-```toml
-[model]
-glacier = true
-```
-
-## Correction factors for forcing data
-The ``e_\mathrm{corr}`` and ``p_\mathrm{corr}`` model parameters can be used to adjust the
-potential evaporation and precipitation, respectively.
-
-## Interception
-After the snow module, rainfall ``P_\mathrm{R}`` [mm t``^{-1}``] and snowmelt
-``P_\mathrm{M}`` [mm t``^{-1}``] enter the interception storage ``S_\mathrm{I}`` \[mm\]. The
-maximum interception storage is defined by the ``I_\mathrm{max}`` \[mm\] parameter for each
-class. Interception evaporation ``E_\mathrm{I}`` [mm t``^{-1}``] occurs at potential rate
-``E_\mathrm{P}`` [mm t``^{-1}``] as long as there is enough water in the interception
-storage. Effective precipitation ``P_\mathrm{E}`` [mm t``^{-1}``] flows out of the
-interception store when the storage capacity is exceeded.  Interception evaporation is
-subtracted from potential evaporation to ensure total evaporation does not exceed potential
-evaporation.
-
-The following equations apply:
-
-```math
-    \mathrm{d}S_\mathrm{I}/\mathrm{d}t = (P_\mathrm{R} + P_\mathrm{M}) - E_\mathrm{I} - P_\mathrm{E}
-```
-
-```math
-    P_\mathrm{E} = \mathrm{max}(0, (S_\mathrm{I} - I_\mathrm{max})/\mathrm{d}t)
-```
-
-```math
-    E_\mathrm{I} = \mathrm{min}(E_\mathrm{P}, S_\mathrm{I}/\mathrm{d}t)
-```
-
-## Hortonion ponding and runoff
-Hortonian overland flow processes are represented by a combination of two storages: a horton
-ponding storage ``S_\mathrm{Hp}`` [mm] and a runoff generating storage ``S_\mathrm{Hf}``
-[mm]. This conceptualization was introduced by de Boer-Euser (2017) and included in the
-plateau class to represent hortonian overland flow (infiltration excess overland flow) in
-agricultural fields. When the storage capacity of the ponding storage is exceeded, runoff
-``Q_\mathrm{H}`` [mm t``^{-1}``] is generated that enters the horton runoff storage. The
-horton runoff generating storage is included to slightly smooth the precipitation signal.
-However, the response time of the runoff generation storage ``K_\mathrm{Hf}`` [t``^{-1}``]
-is very short to generate fast runoff ``Q_\mathrm{Hf}`` [mm t``^{-1}``].
-
-Effective precipitation ``P_\mathrm{E}`` [mm t``^{-1}``] from the interception module enters
-the horton ponding storage, which has a maximum storage capacity ``S_\mathrm{Hmax}`` [mm].
-When the inflow exceeds the storage capacity, direct runoff is generated
-``Q_\mathrm{H,direct}`` [mm t``^{-1}``] and net infiltration in the horton ponding storage
-is denoted as ``Q_\mathrm{H,in,net}`` [mm t``^{-1}``].
-
-Evaporation from the horton ponding storage ``E_\mathrm{H}`` [mm t``^{-1}``]  is based on a
-simple formulation to express water stress. The equation describes how actual evaporation is
-linearly reduced when the relative horton ponding storage ``\overline{S_\mathrm{Hp}}`` [-]
-is below a certain threshold ``L_\mathrm{P}`` [-] parameter.
-
-A beta function with parameter ``\beta`` [-] is used to split the net infiltrating water to
-storage and to runoff ``Q_\mathrm{H}`` [mm t``^{-1}``], to which is added the direct runoff
-``Q_\mathrm{H,direct}`` [mm t``^{-1}``].
-
-The shape of the beta function for various values of ``\beta`` [-] is shown below:
-
-```@setup plot
-    using Printf
-    using CairoMakie
-```
-
-```@example plot
-    let                                                                                     # hide
-        fig = Figure(resolution = (800, 400))                                               # hide
-        ax = Axis(fig[1, 1], xlabel = "S/Smax [-]", ylabel = "Fraction of runoff [-]")      # hide
-        x = 0:0.01:1                                                                        # hide
-        betas = [0.3, 1.0, 3.0]                                                             # hide
-        for beta in betas                                                                   # hide
-            lines!(ax, x, (1 .- (1 .- x).^beta), label = @sprintf("beta = %.1f", beta))     # hide
-        end                                                                                 # hide
-        Legend(fig[1, 2], ax, "beta")                                                       # hide
-        fig                                                                                 # hide
-    end                                                                                     # hide
-```
-
-Additionally, water infiltrates from the horton ponding storage to the root zone storage
-(``Q_\mathrm{HR}`` [mm t``^{-1}``]), based on a formulation using a maximum infiltration
-capacity ``F_\mathrm{max}`` [mm t``^{-1}``] and a decay coefficient ``F_\mathrm{dec}`` [-].
-
-The maximum storage capacity of the horton ponding storage is reduced when soils are likely
-to be frozen ``S_\mathrm{Hmax,frost}`` [mm], i.e. during periods when the temperature is
-below zero for several consecutive days.
-
-The following equations apply for the Horton ponding storage ``S_\mathrm{Hp}`` [mm]:
-
-```math
-    \mathrm{d}S_\mathrm{Hp}/\mathrm{d}t = P_\mathrm{E} - E_\mathrm{H}  - Q_\mathrm{H} - Q_\mathrm{HR}
-```
-
-The constitutive equations are:
-
-```math
-    Q_\mathrm{H,direct}=\mathrm{max}((S_\mathrm{Hp}+P_\mathrm{E}−S_\mathrm{Hmax});0.0)
-```
-
-```math
-    Q_\mathrm{H,in,net} = P_\mathrm{E} − Q_\mathrm{H,direct}
-```
-
-```math
-    \overline{S_\mathrm{Hp}} = S_\mathrm{Hp}/S_\mathrm{H,max}
-```
-
-```math
-    E_\mathrm{H} = \mathrm{min} ( (E_\mathrm{P} - E_\mathrm{I}) \cdot \mathrm{min}(\overline{S_\mathrm{Hp}}/L_\mathrm{P},1), S_\mathrm{Hp}/\mathrm{d}t  )
-```
-
-```math
-    Q_\mathrm{H} = Q_\mathrm{H,in,net} \cdot (1-(1-\overline{S_\mathrm{Hp}})^\beta)
-```
-
-```math
-    Q_\mathrm{HR} = F_\mathrm{max} \cdot \mathrm{exp}(-F_\mathrm{dec} \cdot \overline{S_\mathrm{Hp}})
-```
-
-The reduction of the storage capacity of the horton ponding storage during frozen soil
-conditions is calculated following the equations provided by de Boer-Euser (2017):
-
-```math
-    S_\mathrm{Hmax,frost} = F_\mathrm{T} \cdot S_\mathrm{Hmax}
-```
-
-with:
-
-```math
-F_\mathrm{T} =
-    \begin{cases}
-      S_\mathrm{Hmin} & \text{if $F_\mathrm{acc,fr} < F_\mathrm{acc,fr0}$}\\
-      \frac{F_\mathrm{acc}}{F_\mathrm{acc,fr1} - F_\mathrm{acc,fr0}} - \frac{F_\mathrm{acc,fr0}}{F_\mathrm{acc,fr1} - F_\mathrm{acc,fr0}}  & \text{if $ F_\mathrm{acc,fr0} \le F_\mathrm{acc,fr} \le F_\mathrm{acc,fr1}$}\\
-      1 & \text{if $F_\mathrm{acc,fr} > F_\mathrm{acc,fr,1}$}
-    \end{cases}
-```
-
-where ``S_\mathrm{Hmin}`` [-], ``F_\mathrm{acc,fr0}`` [degree t], ``F_\mathrm{acc,fr1}``
-[degree t] and ``K_\mathrm{mf}`` [-] are all model parameters to describe: a coefficient to
-reduce ``S_\mathrm{Hmax}`` to a minimum storage capacity, the minimum and maximum modelled
-accumulated frost and a melt coefficient for the frozen topsoil, respectively.
-
-The following equations apply for the Horton fast runoff storage ``S_\mathrm{Hf}`` [mm]:
-
-```math
-    \mathrm{d}S_\mathrm{Hf}/\mathrm{d}t = Q_\mathrm{H} - Q_\mathrm{Hf}
-```
-
-```math
-    Q_\mathrm{Hf} = K_\mathrm{Hf}^{-1} \cdot S_\mathrm{Hf}
-```
-
-## Root zone soil moisture
-The incoming water from the interception and hortonian routines ``Q_\mathrm{HR}`` [mm
-t``^{-1}``] enters the root zone storage ``S_\mathrm{R}`` [mm]. The root zone storage has a
-maximum capacity ``S_\mathrm{Rmax}`` [mm], which represents the volume of water in the
-unsaturated root zone, which is available to the roots of vegetation for transpiration.
-Abundant water which exceeds the capacity of the root zone storage cannot infiltrate and
-becomes directly available for runoff ``Q_\mathrm{R,direct}`` [mm t``^{-1}``]. The net
-infiltration in the root zone storage is denoted as ``Q_\mathrm{R,in,net}`` [mm t``^{-1}``].
-
-A simple formulation to express water stress is used to calculate evaporation
-``E_\mathrm{R}`` [mm t``^{-1}``] from the root zone storage. The equation describes how
-actual evaporation is linearly reduced when the relative root zone storage
-``\overline{S_\mathrm{R}}`` [-] is below a certain threshold ``L_\mathrm{P}`` [-] parameter.
-
-Next, a beta function with parameter ``\beta`` [-] describes the partitioning of incoming
-water to the root zone storage and to runoff ``Q_\mathrm{R}`` [mm t``^{-1}``]. The water
-that leaves the root zone storage is partitioned into the fast storage and through
-preferential recharge to the slow storage, based on a splitter parameter ``d_\mathrm{s}``
-[-].
-
-Water may also leave the root zone storage through percolation to the slow groundwater
-``Q_\mathrm{perc}`` [mm t``^{-1}``] or enter the root zone storage from the slow groundwater
-through capillary rise ``Q_\mathrm{cap}`` [mm t``^{-1}``], based on a maximum percolation
-parameter ``Q_\mathrm{perc,max}`` [mm t``^{-1}``] and a maximum capillary rise flux
-parameter ``Q_\mathrm{cap,max}`` [mm t``^{-1}``].
-
-
-The water balance equation for the root zone storage is:
-
-```math
-    \mathrm{d}S_\mathrm{R}/\mathrm{d}t = Q_\mathrm{HR} - E_\mathrm{R}  - Q_\mathrm{R} - Q_\mathrm{perc} + Q_\mathrm{cap}
-```
-
-The constitutive equations are:
-
-```math
-    Q_\mathrm{R,direct}=\mathrm{max}((S_\mathrm{R}+Q_\mathrm{HR}−S_\mathrm{Rmax});0.0)
-```
-
-```math
-    Q_\mathrm{R,in,net} = Q_\mathrm{HR} − Q_\mathrm{R,direct}
-```
-
-```math
-    \overline{S_\mathrm{R}} = S_\mathrm{R}/S_\mathrm{R,max}
-```
-
-```math
-    E_\mathrm{R} = \mathrm{min} ( (E_\mathrm{P} - E_\mathrm{I} - E_\mathrm{H}) \cdot \mathrm{min}(\overline{S_\mathrm{R}}/L_\mathrm{P},1), S_\mathrm{R}/\mathrm{d}t  )
-```
-
-```math
-    Q_\mathrm{R} = Q_\mathrm{R,in,net} \cdot (1-(1-\overline{S_\mathrm{R}})^\beta)
-```
-
-```math
-    Q_\mathrm{perc} = Q_\mathrm{perc,max} \cdot \overline{S_\mathrm{R}}
-```
-
-```math
-    Q_\mathrm{cap} = Q_\mathrm{cap,max} \cdot (1 - \overline{S_\mathrm{R}})
-```
-
-## Fast storage and runoff
-The outflow from the root zone storage ``Q_\mathrm{R}`` [mm t``^{-1}``] is split with the
-splitter parameter ``d_\mathrm{s}`` [-] into inflow in the fast storage ``Q_\mathrm{RF}``
-[mm t``^{-1}``] and inflow in the slow storage ``Q_\mathrm{RS}`` [mm t``^{-1}``] to
-represent preferential recharge. The fast runoff storage ``S_\mathrm{F}`` [mm] generates
-fast runoff ``Q_\mathrm{F}`` [mm t``^{-1}``] through a simple non-linear equation with a
-recession constant ``K_\mathrm{F}`` [t``^{-1}``] and an exponent ``\alpha`` [-].
-
-The following equations apply:
-
-```math
-    \mathrm{d}S_\mathrm{F}/\mathrm{d}t = Q_\mathrm{RF} - Q_\mathrm{F}
-```
-
-```math
-    Q_\mathrm{RF} = Q_\mathrm{R} \cdot (1-d_\mathrm{s})
-```
-
-```math
-    Q_\mathrm{F} = K_\mathrm{F}^{-1} \cdot S_\mathrm{F}^{\alpha}
-```
-
-
-## Common slow groundwater storage and runoff
-The slow groundwater storage ``S_\mathrm{S}`` [mm] is a shared storage for all the different
-classes. It is filled through preferential recharge from the outflow of the root zone
-storage ``Q_\mathrm{RS}`` [mm t``^{-1}``] and through percolation ``Q_\mathrm{perc}`` [mm
-t``^{-1}``]. It empties through capillary rise ``Q_\mathrm{cap}`` [mm t``^{-1}``] and
-through a linear outflow ``Q_\mathrm{S}`` [mm t``^{-1}``] with recession timescale
-coefficient ``K_\mathrm{S}`` [t``^{-1}``].
-
-Total streamflow ``Q_\mathrm{TOT}`` [mm t``^{-1}``] is the weighted sum of the horton fast
-runoff and the fast runoff from the different classes based on the fraction of each class in
-a cell ``F_\mathrm{hrufrac}`` [-] and the slow runoff, which is then routed downstream along
-the river network through the kinematic wave.
-
-```math
-    \mathrm{d}S_\mathrm{S}/\mathrm{d}t = Q_\mathrm{RS} + Q_\mathrm{perc} - Q_\mathrm{S} - Q_\mathrm{cap}
-```
-
-```math
-    Q_\mathrm{RS} = Q_\mathrm{R} \cdot d_\mathrm{s}
-```
-
-```math
-    Q_\mathrm{S} = K_\mathrm{S}^{-1} \cdot S_\mathrm{S}
-```
-
-```math
-Q_\mathrm{TOT} =  Q_\mathrm{S} + \sum_{class=1}^{n} (Q_\mathrm{F,class} + Q_\mathrm{Hf,class}) \cdot F_\mathrm{hrufrac,class}
-```
-
-## References
-+ de Boer-Euser, T. (2017). Added value of distribution in rainfall-runoff models for the Meuse basin
-   PhD thesis, Delft University of Technology. https://doi.org/10.4233/uuid:89a78ae9-7ffb-4260-b25d-698854210fa8
-
-+ Bouaziz, L. J. E., Aalbers, E. E., Weerts, A. H., Hegnauer, M., Buiteveld, H., Lammersen, R., Stam, J., Sprokkereef, E., Savenije, H. H. G., and Hrachowitz, M. (2022)
-    Ecosystem adaptation to climate change: the sensitivity of hydrological predictions to time-dynamic model parameters,
-    Hydrol. Earth Syst. Sci., 26, 1295–1318, https://doi.org/10.5194/hess-26-1295-2022
-
-+ Euser, T., Hrachowitz, M., Winsemius, H. C., & Savenije, H. H. G. (2015).
-    The effect of forcing and landscape distribution on performance and consistency of model structures.
-    Hydrological Processes, 29(17), 3727–3743. https://doi.org/10.1002/hyp.10445
-
-+ Gao, H., Hrachowitz, M., Fenicia, F., Gharari, S., & Savenije, H. H. G. (2014).
-    Testing the realism of a topography-driven model (FLEX-Topo) in the nested catchments of the Upper Heihe, China.
-    Hydrology and Earth System Sciences, 18(5), 1895–1915. https://doi.org/10.5194/hess-18-1895-2014
-
-+ Hanus, S., Hrachowitz, M., Zekollari, H., Schoups, G., Vizcaino, M., and Kaitna, R. (2021)
-    Future changes in annual, seasonal and monthly runoff signatures in contrasting Alpine catchments in Austria,
-    Hydrol. Earth Syst. Sci., 25, 3429–3453, https://doi.org/10.5194/hess-25-3429-2021
-
-+ Hrachowitz, M., Stockinger, M., Coenders-Gerrits, M., van der Ent, R., Bogena, H., Lücke, A., and Stumpp, C. (2021)
-    Reduction of vegetation-accessible water storage capacity after deforestation affects catchment travel time distributions
-    and increases young water fractions in a headwater catchment, Hydrol. Earth Syst. Sci., 25, 4887–4915, https://doi.org/10.5194/hess-25-4887-2021
-
-+ Hulsman, P., Savenije, H. H. G., & Hrachowitz, M. (2021). Learning from satellite observations:
-    Increased understanding of catchment processes through stepwise model improvement.
-    Hydrology and Earth System Sciences, 25(2), 957–982. https://doi.org/10.5194/hess-25-957-2021
-
-+ Savenije, H. H. G. (2010). HESS opinions “topography driven conceptual modelling (FLEX-Topo).”
-    Hydrology and Earth System Sciences, 14(12), 2681–2692. https://doi.org/10.5194/hess-14-2681-2010
diff --git a/docs/src/model_docs/vertical/hbv.md b/docs/src/model_docs/vertical/hbv.md
deleted file mode 100644
index a0c2ae92f..000000000
--- a/docs/src/model_docs/vertical/hbv.md
+++ /dev/null
@@ -1,195 +0,0 @@
-# [HBV](@id vert_hbv)
-
-## Introduction
-This section describes the different vertical processes available as part of the vertical
-HBV concept. This concept is part of the wflow\_hbv model.
-
-## Snow
-The snow model is described in [Snow and glaciers](@ref snow_and_glac).
-
-## Glaciers
-Glacier processes are described in [Snow and glaciers](@ref snow_and_glac). Glacier modelling is enabled
-by specifying the following in the TOML file:
-
-```toml
-[model]
-glacier = true
-```
-## Potential Evaporation
-The `cevpf` model parameter is used to adjust the potential evaporation based on land use.
-In the original HBV version `cevpfo` is used, a factor for forest land use only.
-
-## Interception
-For interception storage a single `icf` parameter is used according to the land use. In this
-implementation interception evaporation is subtracted to ensure total evaporation does not
-exceed potential evaporation. From this storage evaporation equal to the potential
-evaporation rate will occur as long as water is available, even if it is stored as snow. All
-water enters this store first, there is no concept of free throughfall (e.g. through gaps in
-the canopy). In the model a running water budget is kept of the interception store:
-
-+ The available storage (`icf`- actual storage) is filled with the water coming from the
-  snow routine (``Q_{in}``)
-+ Any surplus water now becomes the new ``Q_{in}``
-+ Interception evaporation is determined as the minimum of the current interception storage
-  and the potential evaporation
-
-## The soil routine
-The incoming water from the snow and interception routines, ``Q_{in}``, is available for
-infiltration in the soil routine. The soil layer has a limited capacity, `fc`, to hold soil
-water, which means if `fc` is exceeded the abundant water cannot infiltrate and,
-consequently, becomes directly available for runoff.
-
-```math
-    Q_{dr}=max((SM+Q_{in}−fc);0.0)
-```
-
-where ``Q_{dr}`` is the abundant soil water (also referred to as direct runoff) and ``SM``
-is the soil moisture content. Consequently, the net amount of water that infiltrates into
-the soil, ``I_{net}``, equals:
-
-```math
-I_{net} = Q_{in} − Q_{dr}
-```
-
-Part of the infiltrating water, ``I_{net}``, will runoff through the soil layer (seepage).
-This runoff volume, ``SP``, is related to the soil moisture content, ``SM``, through the
-following power relation:
-
-```math
-SP = \left(\frac{SM}{fc}\right)^\beta I_{net}
-```
-
-where ``\beta`` is an empirically based parameter. Application of this equation implies that
-the amount of seepage water increases with increasing soil moisture content. The fraction of
-the infiltrating water which does not runoff, ``I_{net}−SP``, is added to the available
-amount of soil moisture, ``SM``. The ``\beta`` parameter affects the amount of supply to the
-soil moisture reservoir that is transferred to the quick response reservoir. Values of
-``\beta`` vary generally between 1 and 3. Larger values of ``\beta`` reduce runoff and
-indicate a higher absorption capacity of the soil.
-
-![hbv-soilmoist.png](../../images/hbv-soilmoist.png)
-
-*Schematic view of the soil moisture routine*
-
-A percentage of the soil moisture will evaporate. This percentage is related to the
-potential evaporation and the available amount of soil moisture:
-
-```math
-    E_a = \frac{SM}{T_m} E_p \, ; \, SM<T_m \\~\\
-    E_a = E_p \, ; \, SM \geq T_m
-```
-
-where ``E_a`` is the actual evaporation, ``E_p`` is the potential evaporation and ``T_m``
-(``\leq fc``) is a user defined threshold, above which the actual evaporation equals the
-potential evaporation. ``T_m`` is defined as ``LP * fc`` in which ``LP`` is a soil dependent
-evaporation factor (``LP\leq 1``).
-
-In the original model (Bergström, 1992), a correction to ``E_a`` is applied in case of
-interception. If ``E_a`` from the soil moisture storage plus ``E_i`` exceeds ``ET_p−E_i``
-(``E_i`` = interception evaporation) then the exceeding part is multiplied by a factor
-(``1-e_{red}``), where the parameter ``e_{red}`` varies between 0 and 1. This correction is
-not present in the wflow\_hbv model.
-
-## The runoff response routine
-
-The volume of water which becomes available for runoff, ``S_{dr}+SP``, is transferred to the
-runoff response routine. In this routine the runoff delay is simulated through the use of a
-number of linear reservoirs.
-
-Two linear reservoirs are defined to simulate the different runoff processes: the upper zone
-(generating quick runoff and interflow) and the lower zone (generating slow runoff). The
-available runoff water from the soil routine (i.e. direct runoff, ``S_{dr}``, and seepage,
-``SP``) in principle ends up in the lower zone, unless the percolation threshold, ``PERC``,
-is exceeded, in which case the redundant water ends up in the upper zone:
-
-```math
-\Delta V_{LZ}=min(PERC;(S_{dr}+SP)) \\~\\
-\Delta V_{UZ}=max(0.0;(S_{dr}+SP−PERC))
-```
-
-where ``V_{UZ}`` is the content of the upper zone, ``V_{LZ}`` is the content of the lower
-zone.
-
-Capillary flow from the upper zone to the soil moisture reservoir is modeled according to:
-
-```math
-Q_{cf}=cflux \; (fc−SM)/fc
-```
-
-where ``cflux`` is the maximum capillary flux [mm day``^{-1}``].
-
-The upper zone generates quick runoff (``Q_q``) using:
-```math
-Q_q=K \, UZ^{(1+\alpha)}
-```
-where ``K`` is the upper zone recession coefficient, and ``\alpha`` determines the amount of
-non-linearity. Within HBV-96, the value of ``K`` is determined from three other parameters:
-``\alpha``, ``KHQ``, and ``HQ`` [mm day``^{-1}``]. The value of ``HQ`` represents an outflow
-rate of the upper zone for which the recession rate is equal to ``KHQ``. If we define
-``UZ_{HQ}`` to be the content of the upper zone at outflow rate ``HQ`` we can write the
-following equation:
-
-```math
-    HQ=K \cdot UZ^{(1+\alpha)} = KHQ \cdot UZ_{HQ}
-```
-
-If we eliminate ``UZ_{HQ}`` we obtain:
-
-```math
-    HQ = K \left(\frac{HQ}{KHQ}\right)^{(1+\alpha)}
-```
-
-Rewriting for ``K`` results in:
-```math
-    K = KQH^{(1−\alpha)} HQ^{−\alpha}
-```
-
-The lower zone is a linear reservoir, which means the rate of slow runoff, ``Q_{LZ}``, which
-leaves this zone during one time step equals:
-
-```math
-Q_{LZ} = K_{LZ} \, V_{LZ}
-```
-
-where ``K_{LZ}`` is the reservoir constant.
-
-The upper zone is also a linear reservoir, but it is slightly more complicated than the
-lower zone because it is divided into two zones: a lower part in which interflow is
-generated and an upper part in which quick flow is generated.
-
-![hbv-upper](../../images/hbv-upper.png)
-
-*Schematic view of the upper zone*
-
-If the total water content of the upper zone, ``V_{UZ}`` is lower than a threshold ``UZ1``,
-the upper zone only generates interflow. On the other hand, if ``V_{UZ}`` exceeds ``UZ1``,
-part of the upper zone water will runoff as quick flow:
-
-```math
-    Q_i = K_i min(UZ1;V_{uz})\\
-    Q_q = K_q max((V_{UZ}−UZ1);0.0)
-```
-
-Where ``Q_i`` is the amount of generated interflow in one time step, ``Q_q`` is the amount
-of generated quick flow in one time step and ``K_i`` and ``K_q`` are reservoir constants for
-interflow and quick flow respectively.
-
-The total runoff rate, ``Q``, is equal to the sum of the three different runoff components:
-
-```math
-    Q = Q_{LZ}+Q_i+Q_q
-```
-
-The runoff behavior in the runoff response routine is controlled by two threshold values
-``P_m`` and ``UZ1`` in combination with three reservoir parameters, ``K_{LZ}``, ``K_i`` and
-``K_q``.
-
-In order to represent the differences in delay times between the three runoff components,
-the reservoir constants have to meet the following requirement:
-```math
-    K_{LZ}<K_i<K_q
-```
-
-## References
- + Bergström, S., 1992, The HBV model–its structure and applications, SMHI Reports RH 4,
-   Norrköping.
\ No newline at end of file
diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md
index 740943c42..0d16a4c12 100644
--- a/docs/src/user_guide/additional_options.md
+++ b/docs/src/user_guide/additional_options.md
@@ -144,13 +144,12 @@ and the kinematic wave components of this model can run on multiple threads. The
 [local inertial model for river flow](@ref config_sbm_gwf_lie_river_land) and the optional
 [local inertial model for river (1D) and land (2D)](@ref config_sbm_gwf_lie_river_land),
 both part of wflow\_sbm, can also run on multiple threads. The threading functionality for
-the kinematic wave may also be useful for models that make (partly) use of this routing
-approach as the [wflow_hbv](@ref config_hbv) model and the wflow\_sbm model [SBM +
-Groundwater flow](@ref config_sbm_gwf). The multi-threading functionality in wflow is
-considered experimental, see also the following
-[issue](https://github.com/Deltares/Wflow.jl/issues/139), where an error was not thrown
-running code multi-threaded. Because of this we advise to start with running a wflow model
-single-threaded (for example during the testing phase of setting up an new wflow model).
+the kinematic wave may also be useful for the wflow\_sbm model [SBM + Groundwater flow](@ref
+config_sbm_gwf). The multi-threading functionality in wflow is considered experimental, see
+also the following [issue](https://github.com/Deltares/Wflow.jl/issues/139), where an error
+was not thrown running code multi-threaded. Because of this we advise to start with running
+a wflow model single-threaded (for example during the testing phase of setting up an new
+wflow model).
 
 For information on how to start Julia with multiple threads we refer to [How to start Julia
 with multiple
diff --git a/docs/src/user_guide/intro.md b/docs/src/user_guide/intro.md
index a42634ed6..55f77a860 100644
--- a/docs/src/user_guide/intro.md
+++ b/docs/src/user_guide/intro.md
@@ -6,7 +6,7 @@ configure your model using the TOML file (model settings) and the netCDF gridded
 
 Sample data and model set up is also included for the Moselle River Basin (a major tributary
 of the Rhine River), which can be used to explore the model software. The model is set up
-for the wflow\_sbm, wflow\_hbv and wflow\_sediment model concepts.
+for the wflow\_sbm and wflow\_sediment model concepts.
 
 Finally, information is provided on setting up your own model, including building a model
 from scratch or alternatively using Deltares [HydroMT](https://github.com/Deltares/hydromt)
diff --git a/docs/src/user_guide/model-setup.md b/docs/src/user_guide/model-setup.md
index 6b8598e62..fccb10638 100644
--- a/docs/src/user_guide/model-setup.md
+++ b/docs/src/user_guide/model-setup.md
@@ -77,17 +77,6 @@ Reservoirs or lakes can be part of the kinematic wave or local inertial model fo
 and input parameters are described in [Reservoirs](@ref reservoir_params) and [Lakes](@ref
 lake_params).
 
-The [wflow\_hbv](@ref config_hbv) model configuration consists besides the river and
-overland flow components of the [HBV](@ref vert_hbv) vertical concept. Input parameters for
-this component are described in the [HBV](@ref params_hbv) section of Model parameters. For
-the river and overland flow components the kinematic wave approach is used.
-
-The [wflow\_flextopo](@ref config_flextopo) model configuration consists besides the river
-and overland flow components of the [FLEXTopo](@ref vert_flextopo) vertical concept. Input
-parameters for this component are described in the [FLEXTopo](@ref params_flextopo) section
-of Model parameters. For the river and overland flow components the kinematic wave approach
-is used.
-
 The [wflow\_sediment](@ref config_sediment) model configuration consists of the vertical
 [Soil Erosion](@ref) concept and the input parameters for this concept are described in the
 [Sediment](@ref params_sediment) section of the Model parameters. The parameters of the
diff --git a/docs/src/user_guide/sample_data.md b/docs/src/user_guide/sample_data.md
index e35ded8b8..1994e15bc 100644
--- a/docs/src/user_guide/sample_data.md
+++ b/docs/src/user_guide/sample_data.md
@@ -8,8 +8,6 @@ are listed in the Table below:
 |:--------------- | ------------------|
 | wflow\_sbm + kinematic wave | [sbm_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sbm_config.toml) |
 | wflow\_sbm + groundwater flow | [sbm\_gwf\_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sbm_gwf_config.toml) |
-| wflow\_hbv | [hbv_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/hbv_config.toml) |
-| wflow\_flextopo | [flextopo_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/flextopo_config.toml) |
 | wflow_sediment | [sediment_config.toml](https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/sediment_config.toml) |
 
 The associated Model files (input static, forcing and state files) can easily be downloaded
@@ -55,42 +53,6 @@ download(forcing, joinpath(datadir, "forcing-sbm-groundwater.nc"))
 download(toml_url, toml_path)
 ```
 
-## wflow\_hbv
-```julia
-# urls to TOML and netCDF of the Moselle example model
-toml_url = "https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/hbv_config.toml"
-staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.1/staticmaps-lahn.nc"
-forcing = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.0/forcing-lahn.nc"
-
-# create a "data" directory in the current directory
-datadir = joinpath(@__DIR__, "data")
-mkpath(datadir)
-toml_path = joinpath(@__DIR__, "hbv_config.toml")
-
-# download resources to current and data dirs
-download(staticmaps, joinpath(datadir, "staticmaps-lahn.nc"))
-download(forcing, joinpath(datadir, "forcing-lahn.nc"))
-download(toml_url, toml_path)
-```
-
-## [wflow\_flextopo](@id wflow_flextopo_data)
-```julia
-# urls to TOML and netCDF of the Meuse example model
-toml_url = "https://raw.githubusercontent.com/Deltares/Wflow.jl/master/test/flextopo_config.toml"
-staticmaps = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.8/staticmaps_flex_meuse.nc"
-forcing = "https://github.com/visr/wflow-artifacts/releases/download/v0.2.8/forcing_meuse.nc"
-
-# create a "data" directory in the current directory
-datadir = joinpath(@__DIR__, "data")
-mkpath(datadir)
-toml_path = joinpath(@__DIR__, "flextopo_config.toml")
-
-# download resources to current and data dirs
-download(staticmaps, joinpath(datadir, "staticmaps_flex_meuse.nc"))
-download(forcing, joinpath(datadir, "forcing_meuse.nc"))
-download(toml_url, toml_path)
-```
-
 ## wflow\_sediment
 ```julia
 # urls to TOML and netCDF of the Moselle example model
diff --git a/docs/src/user_guide/step1_requirements.md b/docs/src/user_guide/step1_requirements.md
index 01c0d0ce5..d15e05a16 100644
--- a/docs/src/user_guide/step1_requirements.md
+++ b/docs/src/user_guide/step1_requirements.md
@@ -16,8 +16,7 @@ via CSV files. A brief overview of the different files:
 
 There are several model configurations supported by wflow. These model configurations
 require slightly different input requirements, yet the general structure is similar for each
-model. A wflow model configuration consists of a `vertical` concept like the [SBM](@ref
-vert_sbm), [HBV](@ref vert_hbv) or [FLEXTOPO](@ref vert_flextopo) in
+model. A wflow model configuration consists of a `vertical` [SBM](@ref vert_sbm) concept in
 combination with `lateral` concepts that control how water is routed for example over the
 land or river domain. For the wflow\_sbm model different model configurations are possible.
 The following model configurations are supported in wflow:
@@ -28,9 +27,7 @@ The following model configurations are supported in wflow:
       optional floodplain)
     - SBM + kinematic wave for subsurface flow + local inertial river (1D) and land (2D)
     - SBM + groundwater flow + kinematic wave for surface flow
- - wflow\_hbv: HBV + kinematic wave for surface routing
- - wflow\_flextopo: FLEXTOPO + kinematic wave for surface routing
- - wflow\_sediment as post processing of wflow\_sbm or wflow\_hbv output
+ - wflow\_sediment as post processing of wflow\_sbm output
 
 In the following pages, some examples will be given on how to prepare a basic wflow\_sbm
 model. Example data for other model configurations is provided in the section with [sample
diff --git a/docs/src/user_guide/step2_settings_file.md b/docs/src/user_guide/step2_settings_file.md
index 490679524..6081e45fb 100644
--- a/docs/src/user_guide/step2_settings_file.md
+++ b/docs/src/user_guide/step2_settings_file.md
@@ -56,7 +56,7 @@ Model specific settings can be included in the model section of the TOML file.
 
 ```toml
 [model]
-type = "sbm"                        # one of ("sbm", "sbm_gwf, "hbv")
+type = "sbm"                        # one of ("sbm" or "sbm_gwf)
 masswasting = false                 # include lateral snow transport in the model, default is false
 snow = false                        # include snow modelling, default is false
 reinit = true                       # cold (reinit = true) or warm state (reinit = false), default is true
@@ -243,16 +243,14 @@ extract data for certain locations (e.g. `gauges`) or areas (e.g. `subcatchment`
 netCDF location names are extracted from these maps. For a specific location (grid cell) a
 `location` is required. For layered model parameters and variables that have an extra
 dimension `layer` and are part of the vertical `sbm` concept it is possible to specify an
-internal layer index (see also example below). For model parameters and variables that have
-an extra dimension `classes` and are part of the vertical `FLEXTopo` concept it is possible
-to specify the class name. If multiple layers or classes are desired, this can be specified
-in separate `[[netcdf.variable]]` entries. Note that the specification of the extra
-dimension is not optional when wflow is integrated with Delft-FEWS, for netCDF scalar data
-an extra dimension is not allowed by the `importNetcdfActivity` of the Delft-FEWS General
-Adapter. In the section [Output CSV section](@ref), similar functionality is available for
-CSV. For integration with Delft-FEWS, see also [Run from Delft-FEWS](@ref run_fews), it is
-recommended to write scalar data to netCDF format since the General Adapter of Delft-FEWS
-can ingest this data format directly.
+internal layer index (see also example below). If multiple layers are desired, this can be
+specified in separate `[[netcdf.variable]]` entries. Note that the specification of the
+extra dimension is not optional when wflow is integrated with Delft-FEWS, for netCDF scalar
+data an extra dimension is not allowed by the `importNetcdfActivity` of the Delft-FEWS
+General Adapter. In the section [Output CSV section](@ref), similar functionality is
+available for CSV. For integration with Delft-FEWS, see also [Run from Delft-FEWS](@ref
+run_fews), it is recommended to write scalar data to netCDF format since the General Adapter
+of Delft-FEWS can ingest this data format directly.
 
 ```toml
 [netcdf]
@@ -301,10 +299,8 @@ extract data for certain locations (e.g. `gauges`) or areas (e.g. `subcatchment`
 case a single entry can lead to multiple columns in the CSV file, which will be of the form
 `header_id`, e.g. `Q_20`, for a gauge with integer ID 20. For layered model parameters and
 variables that have an extra dimension `layer` and are part of the vertical `sbm` concept an
-internal layer index (see also example below) should be specified. For model parameters and
-variables that have an extra dimension `classes` and are part of the vertical `FLEXTopo`
-concept it is possible to specify the class name. If multiple layers or classes are desired,
-this can be specified in separate `[[csv.column]]` entries.
+internal layer index (see also example below) should be specified. If multiple layers are
+desired, this can be specified in separate `[[csv.column]]` entries.
 
 The double brackets in `[[csv.column]]` is TOML syntax to indicate that it is part of a
 list. You may specify as many entries as you wish.
diff --git a/src/Wflow.jl b/src/Wflow.jl
index 15784641d..199d2d0a1 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -106,23 +106,17 @@ end
 # different model types (used for dispatch)
 struct SbmModel end         # "sbm" type / sbm_model.jl
 struct SbmGwfModel end      # "sbm_gwf" type / sbm_gwf_model.jl
-struct HbvModel end         # "hbv" type / hbv_model.jl
-struct FlextopoModel end    # "flextopo" type / flextopo_model.jl
 struct SedimentModel end    # "sediment" type / sediment_model.jl
 
 # prevent a large printout of model components and arrays
 Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))
 
 include("horizontal_process.jl")
-include("hbv.jl")
 include("sbm.jl")
-include("flextopo.jl")
 include("flow.jl")
 include("sediment.jl")
 include("reservoir_lake.jl")
-include("hbv_model.jl")
 include("sbm_model.jl")
-include("flextopo_model.jl")
 include("sediment_model.jl")
 include("vertical_process.jl")
 include("groundwater/connectivity.jl")
@@ -188,12 +182,8 @@ function run(config::Config)
         initialize_sbm_model(config)
     elseif modeltype == "sbm_gwf"
         initialize_sbm_gwf_model(config)
-    elseif modeltype == "hbv"
-        initialize_hbv_model(config)
     elseif modeltype == "sediment"
         initialize_sediment_model(config)
-    elseif modeltype == "flextopo"
-        initialize_flextopo_model(config)
     else
         error("unknown model type")
     end
diff --git a/src/bmi.jl b/src/bmi.jl
index a2a7b9899..040a474b0 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -27,12 +27,8 @@ function BMI.initialize(::Type{<:Model}, config_file)
         initialize_sbm_model(config)
     elseif modeltype == "sbm_gwf"
         initialize_sbm_gwf_model(config)
-    elseif modeltype == "hbv"
-        initialize_hbv_model(config)
     elseif modeltype == "sediment"
         initialize_sediment_model(config)
-    elseif modeltype == "flextopo"
-        initialize_flextopo_model(config)
     else
         error("unknown model type")
     end
diff --git a/src/flextopo.jl b/src/flextopo.jl
deleted file mode 100644
index d547a34b9..000000000
--- a/src/flextopo.jl
+++ /dev/null
@@ -1,882 +0,0 @@
-@get_units @exchange @grid_type @grid_location @with_kw struct FLEXTOPO{T,N}
-    # Model time step [s]
-    dt::T | "s" | 0 | "none" | "none"
-    # Number of classes
-    nclass::Int | "-" | 0 | "none" | "none"
-    # Number of cells
-    n::Int | "-" | 0 | "none" | "none"
-    #dictionary with all possible functions for each store
-    dic_function::Dict | "-" | 0 | "none" | "none"
-    #current class
-    kclass::Vector{Int64} | "-" | 0 | "none" | "none"
-    classes::Vector{String} | "-" | 0 | "none" | "none"
-    select_snow::Vector{String} | "-" | 0 | "none" | "none"
-    select_interception::Vector{String} | "-" | 0 | "none" | "none"
-    select_hortonponding::Vector{String} | "-" | 0 | "none" | "none"
-    select_hortonrunoff::Vector{String} | "-" | 0 | "none" | "none"
-    select_rootzone::Vector{String} | "-" | 0 | "none" | "none"
-    select_fast::Vector{String} | "-" | 0 | "none" | "none"
-    select_slow::Vector{String} | "-" | 0 | "none" | "none"
-
-    #fraction of each class
-    hrufrac::Vector{SVector{N,T}} | "-"
-    ## PARAMETERS
-    ##SNOW
-    # Correction factor for precipitation [-]
-    pcorr::Vector{T} | "-"
-    # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹]
-    cfmax::Vector{T} | "mm ᵒC-1 dt-1"
-    # Threshold temperature for snowfall [ᵒC]
-    tt::Vector{T} | "ᵒC"
-    # Threshold temperature interval length [ᵒC]
-    tti::Vector{T} | "ᵒC"
-    # Threshold temperature for snowmelt [ᵒC]
-    ttm::Vector{T} | "ᵒC"
-    # Water holding capacity as fraction of current snow pack [-]
-    whc::Vector{T} | "-"
-    # Refreezing efficiency constant in refreezing of freewater in snow [-]
-    cfr::Vector{T} | "-"
-    # Correction factor for precipitation [-]
-    rfcf::Vector{T} | "-"
-    # Correction factor for snowfall [-]
-    sfcf::Vector{T} | "-"
-    ## GLACIER
-    # Threshold temperature for snowfall above glacier [ᵒC]
-    g_tt::Vector{T} | "ᵒC"
-    # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
-    g_cfmax::Vector{T} | "mm ᵒC-1 dt-1"
-    # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
-    g_sifrac::Vector{T} | "dt-1"
-    # Water within the glacier [mm]
-    glacierstore::Vector{T} | "mm"
-    # Fraction covered by a glacier [-]
-    glacierfrac::Vector{T} | "-"
-    ##INTERCEPTION
-    # Maximum interception storage (in forested and non-forested areas) [mm]
-    imax::Vector{SVector{N,T}} | "mm"
-    # Evap correction [-]
-    ecorr::Vector{T} | "-"
-    ##HORTON
-    # Maximum storage capacity in the hortonian ponding storage [mm]
-    shmax::Vector{SVector{N,T}} | "mm"
-    #recession coefficient of the hortonian runoff storage [dt-1]
-    khf::Vector{SVector{N,T}} | "dt-1"
-    #maximum modelled accumulated frost resulting in shmin [ᵒC Δt]
-    facc0::Vector{SVector{N,T}} | "ᵒC"
-    #minimum modelled accumulated frost resulting in shmax [ᵒC Δt]
-    facc1::Vector{SVector{N,T}} | "ᵒC"
-    #exponent for the decline of infiltration capacity [-]
-    fdec::Vector{SVector{N,T}} | "-"
-    #maximum infiltration capacity from horton ponding [mm dt-1]
-    fmax::Vector{SVector{N,T}} | "mm dt-1"
-    #minimum storage capacity in horton ponding (relative to shmax) [-]
-    shmin::Vector{SVector{N,T}} | "-"
-    #melt coefficient for melt of frozen topsoil [-]
-    kmf::Vector{SVector{N,T}} | "-"
-    ##ROOTZONE
-    # maximum root-zone storage capacity [mm]
-    srmax::Vector{SVector{N,T}} | "mm"
-    # Fraction of root zone storage below which actual evaporation is potential evaporation [-]
-    lp::Vector{SVector{N,T}} | "-"
-    # Exponent in soil runoff generation equation [-]
-    beta::Vector{SVector{N,T}} | "-"
-    # maximum percolation rate [mm Δt⁻¹]
-    perc::Vector{SVector{N,T}} | "mm dt-1"
-    # maximum capillary rise rate [mm Δt⁻¹]
-    cap::Vector{SVector{N,T}} | "mm dt-1"
-    #FAST
-    # Exponent for non linear recession [-]
-    alfa::Vector{SVector{N,T}} | "-"
-    #recession coefficient of fast storage [dt-1]
-    kf::Vector{SVector{N,T}} | "dt-1"
-    # fraction of qrootzone to slowstorage (1-ds to faststorage)
-    ds::Vector{SVector{N,T}} | "-"
-    # SLOW
-    #recession coefficient of slow storage [dt-1]
-    ks::Vector{T} | "dt-1"
-
-    ## STATES
-    ##SNOW
-    # Snow water equivalent [mm]
-    snow::Vector{T} | "mm"
-    # Liquid water content in the snow pack [mm]
-    snowwater::Vector{T} | "mm"
-    # Interception storage [mm]
-    interceptionstorage::Vector{SVector{N,T}} | "mm"
-    # Storage in the hortonian ponding reservoir [mm]
-    hortonpondingstorage::Vector{SVector{N,T}} | "mm"
-    # Storage in the hortonian runoff generation [mm]
-    hortonrunoffstorage::Vector{SVector{N,T}} | "mm"
-    # Storage in the root-zone [mm]
-    rootzonestorage::Vector{SVector{N,T}} | "mm"
-    # Storage in the root-zone relative to maximum root-zone storage capacity [mm]
-    srootzone_over_srmax::Vector{SVector{N,T}} | "mm"
-    # Storage in the fast store [mm]
-    faststorage::Vector{SVector{N,T}} | "mm"
-    # Storage in the slow reservoir (for qcapillary calc) [mm]
-    slowstorage::Vector{T} | "mm"
-    #states previous time step to calc water balance [mm]
-    states_::Vector{SVector{N,T}} | "mm"
-    #states previous time step to calc water balance combined based on perc class. [mm]
-    states_m::Vector{T} | "mm"
-    #states averaged over classes
-    interceptionstorage_m::Vector{T} | "mm"
-    hortonpondingstorage_m::Vector{T} | "mm"
-    hortonrunoffstorage_m::Vector{T} | "mm"
-    srootzone_m::Vector{T} | "mm"
-    faststorage_m::Vector{T} | "mm"
-    srootzone_over_srmax_m::Vector{T} | "mm"
-
-    ## FLUXES
-    #SNOW
-    # Precipitation [mm Δt⁻¹]
-    precipitation::Vector{T} | "mm dt-1"
-    # Temperature [ᵒC]
-    temperature::Vector{T} | "ᵒC"
-    # Potential evapotranspiration [mm Δt⁻¹]
-    potential_evaporation::Vector{T} | "mm dt-1"
-    # Potential evapotranspiration corrected [mm Δt⁻¹]
-    epotcorr::Vector{T} | "mm dt-1"
-    # Precipitation corrected [mm Δt⁻¹]
-    precipcorr::Vector{T} | "mm dt-1"
-    # Snow melt + precipitation as rainfall [mm]
-    rainfallplusmelt::Vector{T} | "mm dt-1"
-    # Snowfall [mm]
-    snowfall::Vector{T} | "mm dt-1"
-    # Snowmelt [mm]
-    snowmelt::Vector{T} | "mm dt-1"
-    #INTERCEPTION
-    # Potential soil evaporation [mm Δt⁻¹]
-    potsoilevap::Vector{SVector{N,T}} | "mm dt-1"
-    # Evaporation from interception storage [mm Δt⁻¹]
-    intevap::Vector{SVector{N,T}} | "mm dt-1"
-    #effective precipitation [mm Δt⁻¹]
-    precipeffective::Vector{SVector{N,T}} | "mm dt-1"
-    # Evaporation from interception sum classes [mm Δt⁻¹]
-    intevap_m::Vector{T} | "mm dt-1"
-    #HORTONPONDING
-    # Evaporation from the hortonion ponding storage [-]
-    hortonevap::Vector{SVector{N,T}} | "mm dt-1"
-    # Flux from the hortonian ponding storage to the hortonian runoff storage [mm Δt⁻¹]
-    qhortonpond::Vector{SVector{N,T}} | "mm dt-1"
-    # Flux from the hortonian ponding storage to the root zone storage [mm Δt⁻¹]
-    qhortonrootzone::Vector{SVector{N,T}} | "mm dt-1"
-    # modeled accumulated frost [ᵒC Δt]
-    facc::Vector{SVector{N,T}} | "ᵒC Δt"
-    # Evaporation from the hortonian sum classes [mm Δt⁻¹]
-    hortonevap_m::Vector{T} | "mm dt-1"
-    #HORTONRUNOFF
-    # Flux from the hortonian runoff storage [mm Δt⁻¹]
-    qhortonrun::Vector{SVector{N,T}} | "mm dt-1"
-    #ROOTZONE
-    # Evaporation from the root-zone storage [mm Δt⁻¹]
-    rootevap::Vector{SVector{N,T}} | "mm dt-1"
-    # Flux from the root-zone storage [mm Δt⁻¹]
-    qrootzone::Vector{SVector{N,T}} | "mm dt-1"
-    # Pref. recharge to fast storage [mm Δt⁻¹]
-    qrootzonefast::Vector{SVector{N,T}} | "mm dt-1"
-    # Pref. recharge to slow storage sum classes [mm Δt⁻¹]
-    qrootzoneslow_m::Vector{T} | "mm dt-1"
-    # Capillary flux from the slow to the root-zone storage [mm Δt⁻¹]
-    qcapillary::Vector{SVector{N,T}} | "mm dt-1"
-    # Capillary flux from the slow to the root-zone storage sum classes [mm Δt⁻¹]
-    qcapillary_m::Vector{T} | "mm dt-1"
-    # Percolation flux from the root-zone to the slow storage [mm Δt⁻¹]
-    qpercolation::Vector{SVector{N,T}} | "mm dt-1"
-    # Percolation flux from the root-zone to the slow storage sum classes [mm Δt⁻¹]
-    qpercolation_m::Vector{T} | "mm dt-1"
-    # Evaporation from the root-zone storage, interception and hortonian [mm Δt⁻¹]
-    actevap::Vector{SVector{N,T}} | "mm dt-1"
-    # Evaporation from the root-zone storage, interception and hortonian sum classes [mm Δt⁻¹]
-    actevap_m::Vector{T} | "mm dt-1"
-    # Evaporation from the root-zone storage sum classes [mm Δt⁻¹]
-    rootevap_m::Vector{T} | "mm dt-1"
-    #FAST
-    # runoff from fast reservoir [mm Δt⁻¹]
-    qfast::Vector{SVector{N,T}} | "mm dt-1"
-    #SLOW
-    # runoff from slow reservoir [mm Δt⁻¹]
-    qslow::Vector{T} | "mm dt-1"
-    #Total [mm Δt⁻¹]
-    runoff::Vector{T} | "mm dt-1"
-    # fast runoff sum classes [mm Δt⁻¹]
-    qfast_tot = Vector{T} | "mm dt-1"
-
-
-    ## WATERBALANCES
-    #water balance snow store
-    wb_snow::Vector{T} | "mm dt-1"
-    #water balance interception storage [mm Δt⁻¹]
-    wb_interception::Vector{SVector{N,T}} | "mm dt-1"
-    #water balance hortonian ponding storage [mm Δt⁻¹]
-    wb_hortonponding::Vector{SVector{N,T}} | "mm dt-1"
-    #water balance hortonian runoff storage [mm Δt⁻¹]
-    wb_hortonrunoff::Vector{SVector{N,T}} | "mm dt-1"
-    #water balance root-zone storage [mm Δt⁻¹]
-    wb_rootzone::Vector{SVector{N,T}} | "mm dt-1"
-    #water balance fast storage [mm Δt⁻¹]
-    wb_fast::Vector{SVector{N,T}} | "mm dt-1"
-    #water balance slow storage [mm Δt⁻¹]
-    wb_slow::Vector{T} | "mm dt-1"
-    #total water balance [mm Δt⁻¹]
-    wb_tot::Vector{T} | "mm dt-1"
-
-end
-
-
-function common_snow_hbv(flextopo::FLEXTOPO)
-    for i = 1:flextopo.n
-        # precip correction
-        precipcorr = flextopo.precipitation[i] * flextopo.pcorr[i]
-
-        #hbv snow
-        snow, snowwater, snowmelt, rainfallplusmelt, snowfall = snowpack_hbv(
-            flextopo.snow[i],
-            flextopo.snowwater[i],
-            precipcorr,
-            flextopo.temperature[i],
-            flextopo.tti[i],
-            flextopo.tt[i],
-            flextopo.ttm[i],
-            flextopo.cfmax[i],
-            flextopo.whc[i],
-        )
-
-        #wb
-        wb_snow =
-            precipcorr - rainfallplusmelt - snow + flextopo.snow[i] - snowwater +
-            flextopo.snowwater[i]
-
-        #update stores
-        # states_ is sum of states of previous time steps to compute WB at the end (per class); states_m is for combined stores
-        flextopo.states_m[i] = flextopo.snow[i] + flextopo.snowwater[i]
-        flextopo.snowfall[i] = snowfall
-        flextopo.snowmelt[i] = snowmelt
-        flextopo.snow[i] = snow
-        flextopo.snowwater[i] = snowwater
-        flextopo.rainfallplusmelt[i] = rainfallplusmelt
-        flextopo.precipcorr[i] = precipcorr
-        flextopo.wb_snow[i] = wb_snow
-    end
-end
-
-function common_snow_no_storage(flextopo::FLEXTOPO)
-    for i = 1:flextopo.n
-        snow = 0.0
-        snowwater = 0.0
-        snowfall = 0.0
-        snowmelt = 0.0
-        # precip correction
-        precipcorr = flextopo.precipitation[i] * flextopo.pcorr[i]
-        rainfallplusmelt = precipcorr
-
-        wb_snow = precipcorr - rainfallplusmelt - snow + flextopo.snow[i]
-
-        #update stores
-        # states_ is sum of states of previous time steps to compute WB at the end (per class); states_m is for combined stores
-        flextopo.states_m[i] = flextopo.snow[i] + flextopo.snowwater[i]
-        flextopo.snowfall[i] = snowfall
-        flextopo.snowmelt[i] = snowmelt
-        flextopo.snow[i] = snow
-        flextopo.snowwater[i] = snowwater
-        flextopo.rainfallplusmelt[i] = rainfallplusmelt
-        flextopo.wb_snow[i] = wb_snow
-        flextopo.precipcorr[i] = precipcorr
-    end
-end
-
-function common_glaciers(flextopo::FLEXTOPO, config)
-    modelglacier = get(config.model, "glacier", false)::Bool
-    for i = 1:flextopo.n
-        if modelglacier
-            # Run Glacier module and add the snowpack on-top of it.
-            # Estimate the fraction of snow turned into ice (HBV-light).
-            # Estimate glacier melt.
-
-            flextopo.snow[i], snow2glacier, flextopo.glacierstore[i], glaciermelt =
-                glacier_hbv(
-                    flextopo.glacierfrac[i],
-                    flextopo.glacierstore[i],
-                    flextopo.snow[i],
-                    flextopo.temperature[i],
-                    flextopo.g_tt[i],
-                    flextopo.g_cfmax[i],
-                    flextopo.g_sifrac[i],
-                    Second(flextopo.dt),
-                )
-            # Convert to mm per grid cell and add to snowmelt
-            glaciermelt = glaciermelt * flextopo.glacierfrac[i]
-            rainfallplusmelt = flextopo.rainfallplusmelt[i] + glaciermelt
-        else
-            rainfallplusmelt = flextopo.rainfallplusmelt[i]
-        end
-    end
-end
-
-
-function interception_no_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        intevap = 0.0
-        precipeffective = max(flextopo.rainfallplusmelt[i], 0.0)
-        interceptionstorage = 0.0
-
-        # correction for potential evaporation
-        epotcorr = flextopo.ecorr[i] * flextopo.potential_evaporation[i]
-        restevap = max(0.0, epotcorr - intevap)
-
-        #wb interceptionstore
-        wb_interception =
-            flextopo.rainfallplusmelt[i] - intevap - precipeffective - interceptionstorage +
-            flextopo.interceptionstorage[i][k]
-
-        #update stores
-        flextopo.states_[i] =
-            setindex(flextopo.states_[i], flextopo.interceptionstorage[i][k], k)
-        flextopo.interceptionstorage[i] =
-            setindex(flextopo.interceptionstorage[i], interceptionstorage, k)
-        flextopo.intevap[i] = setindex(flextopo.intevap[i], intevap, k)
-        flextopo.precipeffective[i] =
-            setindex(flextopo.precipeffective[i], precipeffective, k)
-        flextopo.potsoilevap[i] = setindex(flextopo.potsoilevap[i], restevap, k)
-        flextopo.epotcorr[i] = epotcorr
-        flextopo.wb_interception[i] =
-            setindex(flextopo.wb_interception[i], wb_interception, k)
-
-        #average storage over classes
-        flextopo.interceptionstorage_m[i] =
-            sum(flextopo.interceptionstorage[i] .* flextopo.hrufrac[i])
-        flextopo.intevap_m[i] = sum(flextopo.intevap[i] .* flextopo.hrufrac[i])
-
-    end
-end
-
-function interception_overflow(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        # rainfall added to interception store interceptionstorage
-        interception = min(
-            flextopo.rainfallplusmelt[i],
-            flextopo.imax[i][k] - flextopo.interceptionstorage[i][k],
-        )
-        # current interception storage
-        interceptionstorage = flextopo.interceptionstorage[i][k] + interception
-        precipeffective = flextopo.rainfallplusmelt[i] - interception
-
-        # correction for potential evaporation
-        epotcorr = flextopo.ecorr[i] * flextopo.potential_evaporation[i]
-
-        # evaporation from interception storage
-        intevap = min(interceptionstorage, epotcorr)
-        interceptionstorage = interceptionstorage - intevap
-        restevap = max(0.0, epotcorr - intevap)
-
-        #wb interceptionstore
-        wb_interception =
-            flextopo.rainfallplusmelt[i] - intevap - precipeffective - interceptionstorage +
-            flextopo.interceptionstorage[i][k]
-
-        #update stores
-        flextopo.states_[i] =
-            setindex(flextopo.states_[i], flextopo.interceptionstorage[i][k], k)
-        flextopo.interceptionstorage[i] =
-            setindex(flextopo.interceptionstorage[i], interceptionstorage, k)
-        flextopo.intevap[i] = setindex(flextopo.intevap[i], intevap, k)
-        flextopo.precipeffective[i] =
-            setindex(flextopo.precipeffective[i], precipeffective, k)
-        flextopo.potsoilevap[i] = setindex(flextopo.potsoilevap[i], restevap, k)
-        flextopo.epotcorr[i] = epotcorr
-        flextopo.wb_interception[i] =
-            setindex(flextopo.wb_interception[i], wb_interception, k)
-
-        #average storage over classes
-        flextopo.interceptionstorage_m[i] =
-            sum(flextopo.interceptionstorage[i] .* flextopo.hrufrac[i])
-        flextopo.intevap_m[i] = sum(flextopo.intevap[i] .* flextopo.hrufrac[i])
-
-    end
-end
-
-
-function hortonponding_no_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        hortonevap = 0.0
-        qhortonpond = 0.0
-        qhortonrootzone = max(flextopo.precipeffective[i][k], 0.0)
-        hortonpondingstorage = 0.0
-
-        wb_hortonponding =
-            flextopo.precipeffective[i][k] - hortonevap - qhortonpond - qhortonrootzone -
-            hortonpondingstorage + flextopo.hortonpondingstorage[i][k]
-
-        #update with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.hortonpondingstorage[i][k],
-            k,
-        )
-        flextopo.qhortonpond[i] = setindex(flextopo.qhortonpond[i], qhortonpond, k)
-        flextopo.qhortonrootzone[i] =
-            setindex(flextopo.qhortonrootzone[i], qhortonrootzone, k)
-        flextopo.hortonpondingstorage[i] =
-            setindex(flextopo.hortonpondingstorage[i], hortonpondingstorage, k)
-        flextopo.hortonevap[i] = setindex(flextopo.hortonevap[i], hortonevap, k)
-        flextopo.wb_hortonponding[i] =
-            setindex(flextopo.wb_hortonponding[i], wb_hortonponding, k)
-
-        #average storage over classes
-        flextopo.hortonpondingstorage_m[i] =
-            sum(flextopo.hortonpondingstorage[i] .* flextopo.hrufrac[i])
-        flextopo.hortonevap_m[i] = sum(flextopo.hortonevap[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function hortonponding(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-
-        #calculate reduction of shmax due to frost
-        facc = min(flextopo.facc[i][k] + flextopo.temperature[i] * flextopo.kmf[i][k], 0.0)
-        ft = min(
-            max(
-                facc / (flextopo.facc1[i][k] - flextopo.facc0[i][k]) -
-                flextopo.facc0[i][k] / (flextopo.facc1[i][k] - flextopo.facc0[i][k]),
-                flextopo.shmin[i][k],
-            ),
-            1.0,
-        )
-        shmax_frost = ft * flextopo.shmax[i][k]
-
-        #add effective precipitation from interception to the horton ponding storage.
-        hortonpondingstorage =
-            flextopo.hortonpondingstorage[i][k] + flextopo.precipeffective[i][k]
-        # if soil is filled until max capacity, additional water runs of directly
-        directrunoff_h = max(hortonpondingstorage - shmax_frost, 0.0)
-        #update hortonpondingstorage
-        hortonpondingstorage = hortonpondingstorage - directrunoff_h
-        #net water which infiltrates in horton (careful, shmax_frost can decrease from one timestep to another)
-        netin_hortonpond = max(flextopo.precipeffective[i][k] - directrunoff_h, 0.0)
-
-        #evaporation from the horton ponding storage
-        hortonevap =
-            hortonpondingstorage > shmax_frost * flextopo.lp[i][k] ?
-            min(hortonpondingstorage, flextopo.potsoilevap[i][k]) :
-            min(
-                flextopo.potsoilevap[i][k] *
-                (hortonpondingstorage / (shmax_frost * flextopo.lp[i][k])),
-            )
-        #update storage
-        hortonpondingstorage = hortonpondingstorage - hortonevap
-        #update restevap
-        restevap = max(0.0, flextopo.potsoilevap[i][k] - hortonevap)
-
-        #excess water from beta function of netin_hortonpond (due to rouding hortonpondingstorage could be slightly larger than shmax_frost, make sure it is always less than 1)
-        qhorton_in =
-            netin_hortonpond * (
-                1.0 - pow(
-                    (1.0 - min(hortonpondingstorage / shmax_frost, 1.0)),
-                    flextopo.beta[i][k],
-                )
-            )
-        #update storage
-        hortonpondingstorage = hortonpondingstorage - qhorton_in
-
-        #total water out of horton ponding consists of directrunoff and excess water from beta function
-        qhortonpond = directrunoff_h + qhorton_in
-
-        #infiltration to root-zone
-        qhortonrootzone =
-            min(hortonpondingstorage / shmax_frost, 1.0) > 0.0 ?
-            min(
-                flextopo.fmax[i][k] * exp(
-                    -flextopo.fdec[i][k] *
-                    (1.0 - min(hortonpondingstorage / shmax_frost, 1.0)),
-                ),
-                hortonpondingstorage,
-            ) : 0.0
-        hortonpondingstorage = hortonpondingstorage - qhortonrootzone
-
-        #water balance
-        wb_hortonponding =
-            flextopo.precipeffective[i][k] - hortonevap - qhortonpond - qhortonrootzone -
-            hortonpondingstorage + flextopo.hortonpondingstorage[i][k]
-
-        #update states
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.hortonpondingstorage[i][k],
-            k,
-        )
-        flextopo.potsoilevap[i] = setindex(flextopo.potsoilevap[i], restevap, k)
-        flextopo.qhortonpond[i] = setindex(flextopo.qhortonpond[i], qhortonpond, k)
-        flextopo.qhortonrootzone[i] =
-            setindex(flextopo.qhortonrootzone[i], qhortonrootzone, k)
-        flextopo.hortonpondingstorage[i] =
-            setindex(flextopo.hortonpondingstorage[i], hortonpondingstorage, k)
-        flextopo.hortonevap[i] = setindex(flextopo.hortonevap[i], hortonevap, k)
-        flextopo.wb_hortonponding[i] =
-            setindex(flextopo.wb_hortonponding[i], wb_hortonponding, k)
-
-        #average storage over classes
-        flextopo.hortonpondingstorage_m[i] =
-            sum(flextopo.hortonpondingstorage[i] .* flextopo.hrufrac[i])
-        flextopo.hortonevap_m[i] = sum(flextopo.hortonevap[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function hortonrunoff_no_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        qhortonrun = 0.0
-        hortonrunoffstorage = 0.0
-
-        wb_hortonrunoff =
-            flextopo.qhortonpond[i][k] - qhortonrun - hortonrunoffstorage +
-            flextopo.hortonrunoffstorage[i][k]
-
-        #update with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.hortonrunoffstorage[i][k],
-            k,
-        )
-        flextopo.qhortonrun[i] = setindex(flextopo.qhortonrun[i], qhortonrun, k)
-        flextopo.hortonrunoffstorage[i] =
-            setindex(flextopo.hortonrunoffstorage[i], hortonrunoffstorage, k)
-        flextopo.wb_hortonrunoff[i] =
-            setindex(flextopo.wb_hortonrunoff[i], wb_hortonrunoff, k)
-
-        #average storage over classes
-        flextopo.hortonrunoffstorage_m[i] =
-            sum(flextopo.hortonrunoffstorage[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function hortonrunoff(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-
-        qhortonrun = min(
-            flextopo.hortonrunoffstorage[i][k],
-            flextopo.hortonrunoffstorage[i][k] * flextopo.khf[i][k],
-        )
-        hortonrunoffstorage =
-            flextopo.hortonrunoffstorage[i][k] + flextopo.qhortonpond[i][k] - qhortonrun
-
-        wb_hortonrunoff =
-            flextopo.qhortonpond[i][k] - qhortonrun - hortonrunoffstorage +
-            flextopo.hortonrunoffstorage[i][k]
-
-        #update with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.hortonrunoffstorage[i][k],
-            k,
-        )
-        flextopo.qhortonrun[i] = setindex(flextopo.qhortonrun[i], qhortonrun, k)
-        flextopo.hortonrunoffstorage[i] =
-            setindex(flextopo.hortonrunoffstorage[i], hortonrunoffstorage, k)
-        flextopo.wb_hortonrunoff[i] =
-            setindex(flextopo.wb_hortonrunoff[i], wb_hortonrunoff, k)
-
-        #average storage over classes
-        flextopo.hortonrunoffstorage_m[i] =
-            sum(flextopo.hortonrunoffstorage[i] .* flextopo.hrufrac[i])
-
-    end
-end
-
-function rootzone_no_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        rootevap = 0.0
-        qpercolation = 0.0
-        qcapillary = 0.0
-        qrootzone =
-            max(flextopo.qhortonrootzone[i][k] + flextopo.rootzonestorage[i][k], 0.0) #if store not empty initial conditions
-        rootzonestorage = 0.0
-
-        #compute water balance rootzonestorage storage
-        wb_rootzone =
-            flextopo.qhortonrootzone[i][k] - rootevap - qrootzone - qpercolation +
-            qcapillary - rootzonestorage + flextopo.rootzonestorage[i][k]
-
-        #update states and fluxes with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.rootzonestorage[i][k],
-            k,
-        )
-        flextopo.qrootzone[i] = setindex(flextopo.qrootzone[i], qrootzone, k)
-        flextopo.rootzonestorage[i] =
-            setindex(flextopo.rootzonestorage[i], rootzonestorage, k)
-        flextopo.srootzone_over_srmax[i] = setindex(
-            flextopo.srootzone_over_srmax[i],
-            flextopo.rootzonestorage[i][k] / flextopo.srmax[i][k],
-            k,
-        )
-        flextopo.qcapillary[i] = setindex(flextopo.qcapillary[i], qcapillary, k)
-        flextopo.qpercolation[i] = setindex(flextopo.qpercolation[i], qpercolation, k)
-        flextopo.rootevap[i] = setindex(flextopo.rootevap[i], rootevap, k)
-        flextopo.actevap[i] = setindex(
-            flextopo.actevap[i],
-            flextopo.rootevap[i][k] + flextopo.intevap[i][k] + flextopo.hortonevap[i][k],
-            k,
-        )
-        flextopo.wb_rootzone[i] = setindex(flextopo.wb_rootzone[i], wb_rootzone, k)
-
-        #average storage over classes
-        flextopo.srootzone_m[i] = sum(flextopo.rootzonestorage[i] .* flextopo.hrufrac[i])
-        flextopo.srootzone_over_srmax_m[i] =
-            sum(flextopo.srootzone_over_srmax[i] .* flextopo.hrufrac[i])
-        flextopo.rootevap_m[i] = sum(flextopo.rootevap[i] .* flextopo.hrufrac[i])
-        flextopo.actevap_m[i] = sum(flextopo.actevap[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function rootzone_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-
-        #added water to the root-zone. NB: if no horton storages: qhortonrootzone is in fact Pe!! (effective precip).
-        rootzonestorage = flextopo.rootzonestorage[i][k] + flextopo.qhortonrootzone[i][k]
-        # if soil is filled until max capacity, additional water runs of directly
-        directrunoff = max(rootzonestorage - flextopo.srmax[i][k], 0.0)
-        #update rootzonestorage
-        rootzonestorage = rootzonestorage - directrunoff
-        #net water which infiltrates in root-zone
-        netin_rootzone = flextopo.qhortonrootzone[i][k] - directrunoff
-
-        #evaporation from the rootzone
-        rootevap =
-            rootzonestorage > flextopo.srmax[i][k] * flextopo.lp[i][k] ?
-            min(rootzonestorage, flextopo.potsoilevap[i][k]) :
-            min(
-                flextopo.potsoilevap[i][k] *
-                (rootzonestorage / (flextopo.srmax[i][k] * flextopo.lp[i][k])),
-            )
-        #update storage
-        rootzonestorage = rootzonestorage - rootevap
-
-        #excess water from beta function of netin_rootzone
-        qrootzone_in =
-            netin_rootzone *
-            (1.0 - pow(1.0 - rootzonestorage / flextopo.srmax[i][k], flextopo.beta[i][k]))
-        #update storage
-        rootzonestorage = rootzonestorage - qrootzone_in
-
-        #total water out of root-zone consists of directrunoff and excess water from beta function
-        qrootzone = directrunoff + qrootzone_in
-
-        #percolation
-        qpercolation = flextopo.perc[i][k] * rootzonestorage / flextopo.srmax[i][k]
-        rootzonestorage = rootzonestorage - qpercolation
-
-        #capillary rise
-        qcapillary = min(
-            flextopo.cap[i][k] * (1.0 - rootzonestorage / flextopo.srmax[i][k]),
-            flextopo.slowstorage[i],
-        )
-        rootzonestorage = rootzonestorage + qcapillary
-
-        #compute water balance rootzonestorage storage
-        wb_rootzone =
-            flextopo.qhortonrootzone[i][k] - rootevap - qrootzone - qpercolation +
-            qcapillary - rootzonestorage + flextopo.rootzonestorage[i][k]
-
-        #update states and fluxes with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.rootzonestorage[i][k],
-            k,
-        )
-        flextopo.qrootzone[i] = setindex(flextopo.qrootzone[i], qrootzone, k)
-        flextopo.rootzonestorage[i] =
-            setindex(flextopo.rootzonestorage[i], rootzonestorage, k)
-        flextopo.srootzone_over_srmax[i] = setindex(
-            flextopo.srootzone_over_srmax[i],
-            flextopo.rootzonestorage[i][k] / flextopo.srmax[i][k],
-            k,
-        )
-        flextopo.qcapillary[i] = setindex(flextopo.qcapillary[i], qcapillary, k)
-        flextopo.qpercolation[i] = setindex(flextopo.qpercolation[i], qpercolation, k)
-        flextopo.rootevap[i] = setindex(flextopo.rootevap[i], rootevap, k)
-        flextopo.actevap[i] = setindex(
-            flextopo.actevap[i],
-            flextopo.rootevap[i][k] + flextopo.intevap[i][k] + flextopo.hortonevap[i][k],
-            k,
-        )
-        flextopo.wb_rootzone[i] = setindex(flextopo.wb_rootzone[i], wb_rootzone, k)
-
-        #average storage over classes
-        flextopo.srootzone_m[i] = sum(flextopo.rootzonestorage[i] .* flextopo.hrufrac[i])
-        flextopo.srootzone_over_srmax_m[i] =
-            sum(flextopo.srootzone_over_srmax[i] .* flextopo.hrufrac[i])
-        flextopo.rootevap_m[i] = sum(flextopo.rootevap[i] .* flextopo.hrufrac[i])
-        flextopo.actevap_m[i] = sum(flextopo.actevap[i] .* flextopo.hrufrac[i])
-    end
-end
-
-
-function fast_no_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-        #make sure ds does not exceed 1
-        ds = flextopo.ds[i][k] > 1.0 ? 1.0 : flextopo.ds[i][k]
-        #calc inflow to faststorage
-        qrootzonefast = flextopo.qrootzone[i][k] * (1.0 - ds)
-        qfast = qrootzonefast + flextopo.faststorage[i][k] #if store not empty initial conditions, make sure to empty
-        faststorage = 0.0
-
-        #compute wb faststorage
-        wb_fast = qrootzonefast - qfast - faststorage + flextopo.faststorage[i][k]
-
-        #update with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.faststorage[i][k],
-            k,
-        )
-        flextopo.qrootzonefast[i] = setindex(flextopo.qrootzonefast[i], qrootzonefast, k)
-        flextopo.qfast[i] = setindex(flextopo.qfast[i], qfast, k)
-        flextopo.faststorage[i] = setindex(flextopo.faststorage[i], faststorage, k)
-        flextopo.wb_fast[i] = setindex(flextopo.wb_fast[i], wb_fast, k)
-
-        #average storage over classes
-        flextopo.faststorage_m[i] = sum(flextopo.faststorage[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function fast_storage(flextopo::FLEXTOPO)
-    k = flextopo.kclass[1]
-    for i = 1:flextopo.n
-
-        #make sure ds does not exceed 1
-        ds = flextopo.ds[i][k] > 1.0 ? 1.0 : flextopo.ds[i][k]
-
-        #split part of the outflow from the root-zone to the fast runoff (and part as preferential recharge to the slow reservoir)
-        qrootzonefast = flextopo.qrootzone[i][k] * (1.0 - ds)
-
-        #fast runoff
-        qfast = min(
-            flextopo.faststorage[i][k],
-            pow(flextopo.faststorage[i][k], flextopo.alfa[i][k]) * flextopo.kf[i][k],
-        )
-        #update store
-        faststorage = flextopo.faststorage[i][k] + qrootzonefast - qfast
-
-        #compute wb faststorage
-        wb_fast = qrootzonefast - qfast - faststorage + flextopo.faststorage[i][k]
-
-        #update with setindex
-        flextopo.states_[i] = setindex(
-            flextopo.states_[i],
-            flextopo.states_[i][k] + flextopo.faststorage[i][k],
-            k,
-        )
-        flextopo.qrootzonefast[i] = setindex(flextopo.qrootzonefast[i], qrootzonefast, k)
-        flextopo.qfast[i] = setindex(flextopo.qfast[i], qfast, k)
-        flextopo.faststorage[i] = setindex(flextopo.faststorage[i], faststorage, k)
-        flextopo.wb_fast[i] = setindex(flextopo.wb_fast[i], wb_fast, k)
-
-        #average storage over classes
-        flextopo.faststorage_m[i] = sum(flextopo.faststorage[i] .* flextopo.hrufrac[i])
-    end
-end
-
-function slow_no_storage(flextopo::FLEXTOPO)
-    for i = 1:flextopo.n
-        #make sure ds does not exceed 1
-        pref_recharge = flextopo.qrootzone[i] .* min.(flextopo.ds[i], 1.0)
-        pref_recharge_sum_classes = sum(pref_recharge .* flextopo.hrufrac[i])
-
-        Qcap_sum_classes = sum(flextopo.qcapillary[i] .* flextopo.hrufrac[i])
-        Qperc_sum_classes = sum(flextopo.qpercolation[i] .* flextopo.hrufrac[i])
-
-        Qsin = pref_recharge_sum_classes + Qperc_sum_classes - Qcap_sum_classes
-        qslow = Qsin + flextopo.slowstorage[i] # if at start store is not empty
-        slowstorage = 0.0
-
-        qfast_tot =
-            sum(flextopo.qfast[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.qhortonrun[i] .* flextopo.hrufrac[i])
-        runoff = max(0.0, qslow + qfast_tot)
-
-        wb_slow = Qsin - qslow - slowstorage + flextopo.slowstorage[i]
-
-        #update
-        flextopo.qfast_tot[i] = qfast_tot
-        flextopo.states_m[i] =
-            flextopo.states_m[i] +
-            sum(flextopo.states_[i] .* flextopo.hrufrac[i]) +
-            flextopo.slowstorage[i]
-        flextopo.runoff[i] = runoff
-        flextopo.qrootzoneslow_m[i] = pref_recharge_sum_classes
-        flextopo.qpercolation_m[i] = Qperc_sum_classes
-        flextopo.qcapillary_m[i] = Qcap_sum_classes
-        flextopo.qslow[i] = qslow
-        flextopo.slowstorage[i] = slowstorage
-        flextopo.wb_slow[i] = wb_slow
-    end
-end
-
-function common_slow_storage(flextopo::FLEXTOPO)
-    for i = 1:flextopo.n
-        #make sure ds does not exceed 1
-        pref_recharge = flextopo.qrootzone[i] .* min.(flextopo.ds[i], 1.0)
-        pref_recharge_sum_classes = sum(pref_recharge .* flextopo.hrufrac[i])
-
-        Qcap_sum_classes = sum(flextopo.qcapillary[i] .* flextopo.hrufrac[i])
-        Qperc_sum_classes = sum(flextopo.qpercolation[i] .* flextopo.hrufrac[i])
-
-        Qsin = pref_recharge_sum_classes + Qperc_sum_classes - Qcap_sum_classes
-        slowstorage = flextopo.slowstorage[i] + Qsin
-
-        qslow = min(flextopo.slowstorage[i], flextopo.slowstorage[i] * flextopo.ks[i])
-        slowstorage = slowstorage - qslow
-
-        qfast_tot =
-            sum(flextopo.qfast[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.qhortonrun[i] .* flextopo.hrufrac[i])
-        runoff = max(0.0, qslow + qfast_tot)
-
-        wb_slow = Qsin - qslow - slowstorage + flextopo.slowstorage[i]
-
-        #update
-        flextopo.qfast_tot[i] = qfast_tot
-        flextopo.states_m[i] =
-            flextopo.states_m[i] +
-            sum(flextopo.states_[i] .* flextopo.hrufrac[i]) +
-            flextopo.slowstorage[i]
-        flextopo.runoff[i] = runoff
-        flextopo.qrootzoneslow_m[i] = pref_recharge_sum_classes
-        flextopo.qpercolation_m[i] = Qperc_sum_classes
-        flextopo.qcapillary_m[i] = Qcap_sum_classes
-        flextopo.qslow[i] = qslow
-        flextopo.slowstorage[i] = slowstorage
-        flextopo.wb_slow[i] = wb_slow
-    end
-end
-
-function watbal(flextopo::FLEXTOPO)
-    for i = 1:flextopo.n
-        states =
-            flextopo.snow[i] +
-            flextopo.snowwater[i] +
-            sum(flextopo.interceptionstorage[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.hortonpondingstorage[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.hortonrunoffstorage[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.rootzonestorage[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.faststorage[i] .* flextopo.hrufrac[i]) +
-            sum(flextopo.slowstorage[i] .* flextopo.hrufrac[i])
-        wb_tot =
-            flextopo.precipitation[i] - (
-                sum(flextopo.intevap[i] .* flextopo.hrufrac[i]) +
-                sum(flextopo.hortonevap[i] .* flextopo.hrufrac[i]) +
-                sum(flextopo.rootevap[i] .* flextopo.hrufrac[i])
-            ) - flextopo.runoff[i] - states + flextopo.states_m[i]
-        #update wb
-        flextopo.wb_tot[i] = wb_tot
-    end
-end
diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
deleted file mode 100644
index 5575facba..000000000
--- a/src/flextopo_model.jl
+++ /dev/null
@@ -1,740 +0,0 @@
-"""
-    initialize_flextopo_model(config::Config)
-
-Initial part of the FlexTopo model concept. Reads the input settings and data as defined in the
-Config object. Will return a Model that is ready to run.
-"""
-function initialize_flextopo_model(config::Config)
-    # unpack the paths to the netCDF files
-    static_path = input_path(config, config.input.path_static)
-
-    reader = prepare_reader(config)
-    clock = Clock(config, reader)
-    dt = clock.dt
-
-    do_reservoirs = get(config.model, "reservoirs", false)::Bool
-    do_lakes = get(config.model, "lakes", false)::Bool
-    do_pits = get(config.model, "pits", false)::Bool
-
-    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
-    kw_land_tstep = get(config.model, "kw_land_tstep", 0)
-    kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool
-
-    classes = get(config.model, "classes", "")
-    nclass = length(classes)
-    @info "Classes are set to names `$(join(classes,", "))` and have size `$nclass`."
-    kclass = [1]    #needed to initialize
-
-    # dictionary of available functions for each store
-    dic_function = Dict{String,Function}(
-        "common_snow_hbv" => common_snow_hbv,
-        "common_snow_no_storage" => common_snow_no_storage,
-        "common_glaciers" => common_glaciers,
-        "interception_overflow" => interception_overflow,
-        "interception_no_storage" => interception_no_storage,
-        "hortonponding" => hortonponding,
-        "hortonponding_no_storage" => hortonponding_no_storage,
-        "hortonrunoff" => hortonrunoff,
-        "hortonrunoff_no_storage" => hortonrunoff_no_storage,
-        "rootzone_storage" => rootzone_storage,
-        "rootzone_no_storage" => rootzone_no_storage,
-        "fast_no_storage" => fast_no_storage,
-        "fast_storage" => fast_storage,
-        "slow_no_storage" => slow_no_storage,
-        "common_slow_storage" => common_slow_storage,
-    )
-
-    select_snow = get(config.model, "select_snow", ["common_snow_hbv"])
-    select_interception =
-        get(config.model, "select_interception", ["interception_overflow"])
-    select_hortonponding =
-        get(config.model, "select_hortonponding", ["hortonponding_no_storage"])
-    select_hortonrunoff =
-        get(config.model, "select_hortonrunoff", ["hortonrunoff_no_storage"])
-    select_rootzone = get(config.model, "select_rootzone", ["rootzone_storage"])
-    select_fast = get(config.model, "select_fast", ["fast_storage"])
-    select_slow = get(config.model, "select_slow", ["common_slow_storage"])
-
-    nc = NCDataset(static_path)
-
-    subcatch_2d = ncread(nc, config, "subcatchment"; optional = false, allow_missing = true)
-    # indices based on catchment
-    inds, rev_inds = active_indices(subcatch_2d, missing)
-    n = length(inds)
-
-
-    # glacier parameters
-    g_tt = ncread(
-        nc,
-        config,
-        "vertical.g_tt";
-        sel = inds,
-        defaults = 0.0,
-        type = Float,
-        fill = 0.0,
-    )
-    g_cfmax =
-        ncread(
-            nc,
-            config,
-            "vertical.g_cfmax";
-            sel = inds,
-            defaults = 3.0,
-            type = Float,
-            fill = 0.0,
-        ) .* (dt / basetimestep)
-
-    g_sifrac =
-        ncread(
-            nc,
-            config,
-            "vertical.g_sifrac";
-            sel = inds,
-            defaults = 0.001,
-            type = Float,
-            fill = 0.0,
-        ) .* (dt / basetimestep)
-    glacierfrac = ncread(
-        nc,
-        config,
-        "vertical.glacierfrac";
-        sel = inds,
-        defaults = 0.0,
-        type = Float,
-        fill = 0.0,
-    )
-    glacierstore = ncread(
-        nc,
-        config,
-        "vertical.glacierstore";
-        sel = inds,
-        defaults = 5500.0,
-        type = Float,
-        fill = 0.0,
-    )
-
-    #snow param (single class)
-    cfmax =
-        ncread(
-            nc,
-            config,
-            "vertical.cfmax";
-            sel = inds,
-            defaults = 3.75653,
-            type = Float,
-        ) .* (dt / basetimestep)
-
-    tt = ncread(nc, config, "vertical.tt"; sel = inds, defaults = -1.41934, type = Float)
-
-    tti = ncread(nc, config, "vertical.tti"; sel = inds, defaults = 1.0, type = Float)
-
-    ttm = ncread(nc, config, "vertical.ttm"; sel = inds, defaults = -1.41934, type = Float)
-
-    whc = ncread(nc, config, "vertical.whc"; sel = inds, defaults = 0.1, type = Float)
-
-    cfr = ncread(nc, config, "vertical.cfr"; sel = inds, defaults = 0.05, type = Float)
-
-    #parameters which are not class specific
-    pcorr = ncread(nc, config, "vertical.pcorr"; sel = inds, defaults = 1.0, type = Float)
-    ecorr = ncread(nc, config, "vertical.ecorr"; sel = inds, defaults = 1.0, type = Float)
-    rfcf = ncread(nc, config, "vertical.rfcf"; sel = inds, defaults = 1.0, type = Float)
-    sfcf = ncread(nc, config, "vertical.sfcf"; sel = inds, defaults = 1.0, type = Float)
-    ks =
-        ncread(nc, config, "vertical.ks"; sel = inds, defaults = 0.006, type = Float) .*
-        (dt / basetimestep)
-
-    #initialize parameters that differ per class
-    hrufrac = ncread(
-        nc,
-        config,
-        "vertical.hrufrac";
-        sel = inds,
-        defaults = 1.0 / length(classes),
-        type = Float,
-        dimname = :classes,
-    )
-    imax = ncread(
-        nc,
-        config,
-        "vertical.imax";
-        sel = inds,
-        defaults = 3.0,
-        type = Float,
-        dimname = :classes,
-    )
-    shmax = ncread(
-        nc,
-        config,
-        "vertical.shmax";
-        sel = inds,
-        defaults = 30.0,
-        type = Float,
-        dimname = :classes,
-    )
-    khf =
-        ncread(
-            nc,
-            config,
-            "vertical.khf";
-            sel = inds,
-            defaults = 0.5,
-            type = Float,
-            dimname = :classes,
-        ) .* (dt / basetimestep)
-    facc0 = ncread(
-        nc,
-        config,
-        "vertical.facc0";
-        sel = inds,
-        defaults = -3.0,
-        type = Float,
-        dimname = :classes,
-    )
-    facc1 = ncread(
-        nc,
-        config,
-        "vertical.facc1";
-        sel = inds,
-        defaults = 0.0,
-        type = Float,
-        dimname = :classes,
-    )
-    fdec = ncread(
-        nc,
-        config,
-        "vertical.fdec";
-        sel = inds,
-        defaults = 0.2,
-        type = Float,
-        dimname = :classes,
-    )
-    fmax =
-        ncread(
-            nc,
-            config,
-            "vertical.fmax";
-            sel = inds,
-            defaults = 2.0,
-            type = Float,
-            dimname = :classes,
-        ) .* (dt / basetimestep)
-    shmin = ncread(
-        nc,
-        config,
-        "vertical.shmin";
-        sel = inds,
-        defaults = 0.2,
-        type = Float,
-        dimname = :classes,
-    )
-    kmf = ncread(
-        nc,
-        config,
-        "vertical.kmf";
-        sel = inds,
-        defaults = 1.0,
-        type = Float,
-        dimname = :classes,
-    )
-    srmax = ncread(
-        nc,
-        config,
-        "vertical.srmax";
-        sel = inds,
-        defaults = 260.0,
-        type = Float,
-        dimname = :classes,
-    )
-    beta = ncread(
-        nc,
-        config,
-        "vertical.beta";
-        sel = inds,
-        defaults = 0.3,
-        type = Float,
-        dimname = :classes,
-    )
-    lp = ncread(
-        nc,
-        config,
-        "vertical.lp";
-        sel = inds,
-        defaults = 0.3,
-        type = Float,
-        dimname = :classes,
-    )
-    perc =
-        ncread(
-            nc,
-            config,
-            "vertical.perc";
-            sel = inds,
-            defaults = 0.30,
-            type = Float,
-            dimname = :classes,
-        ) .* (dt / basetimestep)
-    cap =
-        ncread(
-            nc,
-            config,
-            "vertical.cap";
-            sel = inds,
-            defaults = 0.20,
-            type = Float,
-            dimname = :classes,
-        ) .* (dt / basetimestep)
-    kf =
-        ncread(
-            nc,
-            config,
-            "vertical.kf";
-            sel = inds,
-            defaults = 0.1,
-            type = Float,
-            dimname = :classes,
-        ) .* (dt / basetimestep)
-    alfa = ncread(
-        nc,
-        config,
-        "vertical.alfa";
-        sel = inds,
-        defaults = 1.0,
-        type = Float,
-        dimname = :classes,
-    )
-    ds = ncread(
-        nc,
-        config,
-        "vertical.ds";
-        sel = inds,
-        defaults = 0.2,
-        type = Float,
-        dimname = :classes,
-    )
-
-    # read x, y coordinates and calculate cell length [m]
-    y_nc = read_y_axis(nc)
-    x_nc = read_x_axis(nc)
-    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    cellength = abs(mean(diff(x_nc)))
-
-
-    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
-    xl, yl = cell_lengths(y, cellength, sizeinmetres)
-
-
-    interceptionstorage = zeros(Float, nclass, n)
-    hortonpondingstorage = zeros(Float, nclass, n)
-    hortonrunoffstorage = zeros(Float, nclass, n)
-    srootzone_over_srmax = zeros(Float, nclass, n) .+ 1.0
-
-    states_ = fill(mv, nclass, n)
-
-    potsoilevap = fill(mv, nclass, n)
-    intevap = fill(mv, nclass, n)
-    precipeffective = fill(mv, nclass, n)
-    hortonevap = fill(mv, nclass, n)
-    rootevap = fill(mv, nclass, n)
-    qhortonpond = fill(mv, nclass, n)
-    qhortonrootzone = fill(mv, nclass, n)
-    facc = zeros(Float, nclass, n)
-    qhortonrun = fill(mv, nclass, n)
-    qrootzone = fill(mv, nclass, n)
-    qrootzonefast = fill(mv, nclass, n)
-    qfast = fill(mv, nclass, n)
-    actevap = fill(mv, nclass, n)
-    qpercolation = fill(mv, nclass, n)
-    qcapillary = fill(mv, nclass, n)
-
-    wb_interception = fill(mv, nclass, n)
-    wb_hortonponding = fill(mv, nclass, n)
-    wb_hortonrunoff = fill(mv, nclass, n)
-    wb_rootzone = fill(mv, nclass, n)
-    wb_fast = fill(mv, nclass, n)
-
-
-    flextopo = FLEXTOPO{Float,nclass}(
-        dt = Float(tosecond(dt)),
-        nclass = nclass,
-        n = n,
-        dic_function = dic_function,
-        kclass = kclass,
-        classes = classes,
-        select_snow = select_snow,
-        select_interception = select_interception,
-        select_hortonponding = select_hortonponding,
-        select_hortonrunoff = select_hortonrunoff,
-        select_rootzone = select_rootzone,
-        select_fast = select_fast,
-        select_slow = select_slow,
-        hrufrac = svectorscopy(hrufrac, Val{nclass}()),
-        pcorr = pcorr,
-        ecorr = ecorr,
-        #glaciers
-        g_tt = g_tt,
-        g_cfmax = g_cfmax,
-        g_sifrac = g_sifrac,
-        glacierfrac = glacierfrac,
-        glacierstore = glacierstore,
-        # snow
-        tti = tti,
-        tt = tt,
-        ttm = ttm,
-        cfmax = cfmax,
-        whc = whc,
-        cfr = cfr,
-        rfcf = rfcf,
-        sfcf = sfcf,
-        #interception
-        imax = svectorscopy(imax, Val{nclass}()),
-        #horton
-        shmax = svectorscopy(shmax, Val{nclass}()),
-        khf = svectorscopy(khf, Val{nclass}()),
-        facc0 = svectorscopy(facc0, Val{nclass}()),
-        facc1 = svectorscopy(facc1, Val{nclass}()),
-        fdec = svectorscopy(fdec, Val{nclass}()),
-        fmax = svectorscopy(fmax, Val{nclass}()),
-        shmin = svectorscopy(shmin, Val{nclass}()),
-        kmf = svectorscopy(kmf, Val{nclass}()),
-        #root zone
-        srmax = svectorscopy(srmax, Val{nclass}()),
-        lp = svectorscopy(lp, Val{nclass}()),
-        beta = svectorscopy(beta, Val{nclass}()),
-        perc = svectorscopy(perc, Val{nclass}()),
-        cap = svectorscopy(cap, Val{nclass}()),
-        #fast
-        ds = svectorscopy(ds, Val{nclass}()),
-        alfa = svectorscopy(alfa, Val{nclass}()),
-        kf = svectorscopy(kf, Val{nclass}()),
-        #slow
-        ks = ks,
-
-        # # default (cold) states:
-        snow = zeros(Float, n),
-        snowwater = zeros(Float, n),
-        interceptionstorage = svectorscopy(interceptionstorage, Val{nclass}()),
-        hortonpondingstorage = svectorscopy(hortonpondingstorage, Val{nclass}()),
-        hortonrunoffstorage = svectorscopy(hortonrunoffstorage, Val{nclass}()),
-        rootzonestorage = svectorscopy(srmax, Val{nclass}()),
-        faststorage = 0.0 .* svectorscopy(srmax, Val{nclass}()),
-        srootzone_over_srmax = svectorscopy(srootzone_over_srmax, Val{nclass}()),
-        slowstorage = zeros(Float, n) .+ 30.0,
-        #states previous time step
-        states_m = fill(mv, n),
-        states_ = svectorscopy(states_, Val{nclass}()),
-        #states averaged over all classes
-        interceptionstorage_m = zeros(Float, n),
-        hortonpondingstorage_m = zeros(Float, n),
-        hortonrunoffstorage_m = zeros(Float, n),
-        srootzone_m = zeros(Float, n),
-        faststorage_m = zeros(Float, n),
-        srootzone_over_srmax_m = zeros(Float, n),
-
-        # variables:
-        precipitation = fill(mv, n),
-        temperature = fill(mv, n),
-        potential_evaporation = fill(mv, n),
-        epotcorr = fill(mv, n),
-        precipcorr = fill(mv, n),
-        rainfallplusmelt = fill(mv, n),
-        snowfall = fill(mv, n),
-        snowmelt = fill(mv, n),
-        potsoilevap = svectorscopy(potsoilevap, Val{nclass}()),
-        precipeffective = svectorscopy(precipeffective, Val{nclass}()),
-        intevap = svectorscopy(intevap, Val{nclass}()),
-        hortonevap = svectorscopy(hortonevap, Val{nclass}()),
-        rootevap = svectorscopy(rootevap, Val{nclass}()),
-        qhortonpond = svectorscopy(qhortonpond, Val{nclass}()),
-        qhortonrootzone = svectorscopy(qhortonrootzone, Val{nclass}()),
-        facc = svectorscopy(facc, Val{nclass}()),
-        qhortonrun = svectorscopy(qhortonrun, Val{nclass}()),
-        qrootzone = svectorscopy(qrootzone, Val{nclass}()),
-        qrootzonefast = svectorscopy(qrootzonefast, Val{nclass}()),
-        qrootzoneslow_m = fill(mv, n),
-        qfast = svectorscopy(qfast, Val{nclass}()),
-        actevap = svectorscopy(actevap, Val{nclass}()),
-        qpercolation = svectorscopy(qpercolation, Val{nclass}()),
-        qcapillary = svectorscopy(qcapillary, Val{nclass}()),
-        qpercolation_m = fill(mv, n),
-        qcapillary_m = fill(mv, n),
-        # combined for the classes
-        actevap_m = fill(mv, n),
-        intevap_m = fill(mv, n),
-        hortonevap_m = fill(mv, n),
-        rootevap_m = fill(mv, n),
-        qslow = fill(mv, n),
-        qfast_tot = fill(mv, n),
-        runoff = fill(mv, n),
-        wb_snow = fill(mv, n),
-        wb_interception = svectorscopy(wb_interception, Val{nclass}()),
-        wb_hortonponding = svectorscopy(wb_hortonponding, Val{nclass}()),
-        wb_hortonrunoff = svectorscopy(wb_hortonrunoff, Val{nclass}()),
-        wb_rootzone = svectorscopy(wb_rootzone, Val{nclass}()),
-        wb_fast = svectorscopy(wb_fast, Val{nclass}()),
-        wb_slow = fill(mv, n),
-        wb_tot = fill(mv, n),
-    )
-
-    modelsize_2d = size(subcatch_2d)
-    river_2d =
-        ncread(nc, config, "river_location"; optional = false, type = Bool, fill = false)
-    river = river_2d[inds]
-    riverwidth_2d =
-        ncread(nc, config, "lateral.river.width"; optional = false, type = Float, fill = 0)
-    riverwidth = riverwidth_2d[inds]
-    riverlength_2d =
-        ncread(nc, config, "lateral.river.length"; optional = false, type = Float, fill = 0)
-    riverlength = riverlength_2d[inds]
-
-    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
-    nriv = length(inds_riv)
-
-    # reservoirs
-    pits = zeros(Bool, modelsize_2d)
-    if do_reservoirs
-        reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
-    else
-        reservoir = ()
-        reservoirs = nothing
-        resindex = fill(0, nriv)
-    end
-
-    # lakes
-    if do_lakes
-        lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
-    else
-        lake = ()
-        lakes = nothing
-        lakeindex = fill(0, nriv)
-    end
-
-    ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
-    ldd = ldd_2d[inds]
-    if do_pits
-        pits_2d = ncread(nc, config, "pits"; optional = false, type = Bool, fill = false)
-        ldd = set_pit_ldd(pits_2d, ldd, inds)
-    end
-
-    landslope =
-        ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(landslope, 0.00001, Inf)
-
-    dl = map(detdrainlength, ldd, xl, yl)
-    dw = (xl .* yl) ./ dl
-    olf = initialize_surfaceflow_land(
-        nc,
-        config,
-        inds;
-        sl = landslope,
-        dl = dl,
-        width = map(det_surfacewidth, dw, riverwidth, river),
-        iterate = kinwave_it,
-        tstep = kw_land_tstep,
-        dt = dt,
-    )
-
-    graph = flowgraph(ldd, inds, pcr_dir)
-
-    riverlength = riverlength_2d[inds_riv]
-    riverwidth = riverwidth_2d[inds_riv]
-
-    ldd_riv = ldd_2d[inds_riv]
-    if do_pits
-        ldd_riv = set_pit_ldd(pits_2d, ldd_riv, inds_riv)
-    end
-    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)
-
-    # the indices of the river cells in the land(+river) cell vector
-    index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
-
-    rf = initialize_surfaceflow_river(
-        nc,
-        config,
-        inds_riv;
-        dl = riverlength,
-        width = riverwidth,
-        reservoir_index = resindex,
-        reservoir = reservoirs,
-        lake_index = lakeindex,
-        lake = lakes,
-        iterate = kinwave_it,
-        tstep = kw_river_tstep,
-        dt = dt,
-    )
-
-    # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
-    # subdomain is equal to the complete domain
-    toposort = topological_sort_by_dfs(graph)
-    toposort_riv = topological_sort_by_dfs(graph_riv)
-    index_pit_land = findall(x -> x == 5, ldd)
-    index_pit_river = findall(x -> x == 5, ldd_riv)
-    streamorder = stream_order(graph, toposort)
-    min_streamorder_land = get(config.model, "min_streamorder_land", 5)
-    subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
-        graph,
-        toposort,
-        index_pit_land,
-        streamorder,
-        min_streamorder_land,
-    )
-    min_streamorder_river = get(config.model, "min_streamorder_river", 6)
-    subriv_order, indices_subriv, topo_subriv = kinwave_set_subdomains(
-        graph_riv,
-        toposort_riv,
-        index_pit_river,
-        streamorder[index_river],
-        min_streamorder_river,
-    )
-
-    modelmap = (vertical = flextopo, lateral = (land = olf, river = rf))
-    indices_reverse = (
-        land = rev_inds,
-        river = rev_inds_riv,
-        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
-        lake = isempty(lake) ? nothing : lake.reverse_indices,
-    )
-
-    # if fews_run = true, then classes are specified as Float64 index so FEWS can import
-    # NetCDF output (3D data/SVector).
-    fews_run = get(config, "fews_run", false)::Bool
-    classes = fews_run ? Float64.(collect(1:length(flextopo.classes))) : flextopo.classes
-
-    writer = prepare_writer(
-        config,
-        modelmap,
-        indices_reverse,
-        x_nc,
-        y_nc,
-        nc,
-        extra_dim = (name = "classes", value = classes),
-    )
-    close(nc)
-
-    # for each domain save:
-    # - the directed acyclic graph (graph),
-    # - the traversion order (order),
-    # - upstream_nodes,
-    # - subdomains for the kinematic wave domains for parallel execution (execution order of
-    #   subbasins (subdomain_order), traversion order per subbasin (topo_subdomain) and
-    #   Vector indices per subbasin matching the traversion order of the complete domain
-    #   (indices_subdomain))
-    # - the indices that map it back to the two dimensional grid (indices)
-
-    # for the land domain the x and y length [m] of the grid cells are stored
-    # for reservoirs and lakes indices information is available from the initialization
-    # functions
-    land = (
-        graph = graph,
-        upstream_nodes = filter_upsteam_nodes(graph, pits[inds]),
-        subdomain_order = subbas_order,
-        topo_subdomain = topo_subbas,
-        indices_subdomain = indices_subbas,
-        order = toposort,
-        indices = inds,
-        reverse_indices = rev_inds,
-        xl = xl,
-        yl = yl,
-    )
-    river = (
-        graph = graph_riv,
-        upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
-        subdomain_order = subriv_order,
-        topo_subdomain = topo_subriv,
-        indices_subdomain = indices_subriv,
-        order = toposort_riv,
-        indices = inds_riv,
-        reverse_indices = rev_inds_riv,
-    )
-
-    model = Model(
-        config,
-        (; land, river, reservoir, lake, index_river, frac_toriver),
-        (subsurface = nothing, land = olf, river = rf),
-        flextopo,
-        clock,
-        reader,
-        writer,
-        FlextopoModel(),
-    )
-
-    model = set_states(model)
-
-    return model
-end
-
-function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
-    @unpack lateral, vertical, network, clock, config = model
-
-    inds_riv = network.index_river
-
-    #COMMON SNOW
-    vertical.dic_function[vertical.select_snow[1]](vertical)
-
-    #lateral snow transport
-    if get(config.model, "masswasting", false)::Bool
-        lateral_snow_transport!(
-            vertical.snow,
-            vertical.snowwater,
-            lateral.land.sl,
-            network.land,
-        )
-    end
-
-    if get(config.model, "glacier", false)::Bool
-        common_glaciers(vertical, config)
-    end
-
-    for (k, class) in enumerate(vertical.classes)
-        vertical.kclass[1] = k
-
-        #INTERCEPTION
-        vertical.dic_function[vertical.select_interception[k]](vertical)
-
-        #HORTON
-        vertical.dic_function[vertical.select_hortonponding[k]](vertical)
-        vertical.dic_function[vertical.select_hortonrunoff[k]](vertical)
-
-        #ROOT-ZONE
-        vertical.dic_function[vertical.select_rootzone[k]](vertical)
-
-        #FAST
-        vertical.dic_function[vertical.select_fast[k]](vertical)
-
-    end
-
-    #COMMON SLOW
-    vertical.dic_function[vertical.select_slow[1]](vertical)
-
-    # WAT BAL
-    watbal(vertical)
-
-    surface_routing(model)
-
-    return model
-end
-
-function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel}
-    @unpack lateral, config = model
-    reinit = get(config.model, "reinit", true)::Bool
-    # read and set states in model object if reinit=true
-    if reinit == false
-        instate_path = input_path(config, config.state.path_input)
-        set_states(instate_path, model; type = Float, dimname = :classes)
-
-        # update kinematic wave volume for river and land domain
-        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
-        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
-
-        if do_lakes
-            # storage must be re-initialized after loading the state with the current
-            # waterlevel otherwise the storage will be based on the initial water level
-            lakes = lateral.river.lake
-            lakes.storage .=
-                initialize_storage(lakes.storfunc, lakes.area, lakes.waterlevel, lakes.sh)
-        end
-    end
-    return model
-end
diff --git a/src/flow.jl b/src/flow.jl
index 31fa480d7..77b058bf0 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -372,7 +372,8 @@ function stable_timestep(sf::S) where {S<:SurfaceFlow}
             courant = zeros(n)
             for v = 1:n
                 if sf.q[v] > 0.0
-                    sf.cel[v] = 1.0 / (sf.alpha[v] * sf.beta * pow(sf.q[v], (sf.beta - 1.0)))
+                    sf.cel[v] =
+                        1.0 / (sf.alpha[v] * sf.beta * pow(sf.q[v], (sf.beta - 1.0)))
                     courant[v] = (sf.dt / sf.dl[v]) * sf.cel[v]
                 end
             end
@@ -1153,8 +1154,8 @@ function stable_timestep(sw::ShallowWaterLand{T})::T where {T}
     dt_min = T(Inf)
     @batch per = thread reduction = ((min, dt_min),) for i = 1:sw.n
         @fastmath @inbounds dt =
-            sw.rivercells[i] == 0 ? sw.alpha * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]) :
-            T(Inf)
+            sw.rivercells[i] == 0 ?
+            sw.alpha * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]) : T(Inf)
         dt_min = min(dt, dt_min)
     end
     dt_min = isinf(dt_min) ? T(10.0) : dt_min
@@ -1604,12 +1605,12 @@ function initialize_floodplain_1d(
 end
 
 """
-    set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T}
+    set_river_inwater(model, ssf_toriver)
 
-Set `inwater` of the lateral river component for a `Model` with vertical `SBM` concept.
-`ssf_toriver` is the subsurface flow to the river.
+Set `inwater` of the lateral river component for a model `ssf_toriver` is the subsurface
+flow to the river.
 """
-function set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T}
+function set_river_inwater(model, ssf_toriver)
     @unpack lateral, vertical, network = model
     inds = network.index_river
 
@@ -1628,17 +1629,6 @@ function set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<
     )
 end
 
-"""
-    set_river_inwater(model, ssf_toriver)
-
-Set `inwater` of the lateral river component (based on overland flow).
-"""
-function set_river_inwater(model, ssf_toriver)
-    @unpack lateral, network = model
-    inds = network.index_river
-    lateral.river.inwater .= lateral.land.to_river[inds]
-end
-
 """
     set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
 
@@ -1671,17 +1661,6 @@ function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmMode
         lateral.land.dt
 end
 
-"""
-    set_land_inwater(model)
-
-Set `inwater` of the lateral land component, based on `runoff` of the `vertical` concept.
-"""
-function set_land_inwater(model)
-    @unpack lateral, vertical, network = model
-    lateral.land.inwater .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.dt
-end
-
 # Computation of inflow from the lateral components `land` and `subsurface` to water bodies
 # depends on the routing scheme (see different `get_inflow_waterbody` below). For the river
 # kinematic wave, the variables `to_river` can be excluded, because this part is added to
diff --git a/src/hbv.jl b/src/hbv.jl
deleted file mode 100644
index 038514622..000000000
--- a/src/hbv.jl
+++ /dev/null
@@ -1,239 +0,0 @@
-@get_units @exchange @grid_type @grid_location @with_kw struct HBV{T}
-    dt::T | "s" | 0 | "none" | "none"   # Model time step [s]
-    n::Int | "-" | 0 | "none" | "none"  # Number of cells
-    fc::Vector{T} | "mm"                # Field capacity [mm]
-    betaseepage::Vector{T} | "-"        # Exponent in soil runoff generation equation [-]
-    lp::Vector{T} | "-"                 # Fraction of field capacity below which actual evaporation=potential evaporation [-]
-    threshold::Vector{T} | "mm"         # Threshold soilwater storage above which AE=PE [mm]
-    k4::Vector{T} | "dt-1"              # Recession constant baseflow [Δt⁻¹]
-    kquickflow::Vector{T} | "dt-1"      # Recession constant upper reservoir [Δt⁻¹]
-    suz::Vector{T} | "mm"               # Level over which k0 is used [mm]
-    k0::Vector{T} | "dt-1"              # Recession constant upper reservoir [Δt⁻¹]
-    khq::Vector{T} | "dt-1"             # Recession rate at flow hq [Δt⁻¹]
-    hq::Vector{T}                       # High flow rate hq for which recession rate of upper reservoir is known [mm Δt⁻¹]
-    alphanl::Vector{T}                  # Measure of non-linearity of upper reservoir
-    perc::Vector{T}                     # Percolation from upper to lower zone [mm Δt⁻¹]
-    cfr::Vector{T} | "-"                # Refreezing efficiency constant in refreezing of freewater in snow [-]
-    pcorr::Vector{T} | "-"              # Correction factor for precipitation [-]
-    rfcf::Vector{T} | "-"               # Correction factor for rainfall [-]
-    sfcf::Vector{T} | "-"               # Correction factor for snowfall [-]
-    cflux::Vector{T}                    # Maximum capillary rise from runoff response routine to soil moisture routine [mm Δt⁻¹]
-    icf::Vector{T} | "mm"               # Maximum interception storage (in forested and non-forested areas) [mm]
-    cevpf::Vector{T} | "-"              # Correction factor for potential evaporation [-]
-    epf::Vector{T} | "mm-1"             # Exponent of correction factor for evaporation on days with precipitation
-    ecorr::Vector{T} | "-"              # Evap correction [-]
-    tti::Vector{T} | "ᵒC"               # Critical temperature for snowmelt and refreezing [ᵒC]
-    tt::Vector{T} | "ᵒC"                # Defines interval in which precipitation falls as rainfall and snowfall [ᵒC]
-    ttm::Vector{T} | "ᵒC"               # Threshold temperature for snowmelt [ᵒC]
-    cfmax::Vector{T} | "mm ᵒC-1 dt-1"   # Meltconstant in temperature-index [-]
-    whc::Vector{T} | "-"                # Fraction of snow volume that can store water [-]
-    g_tt::Vector{T} | "ᵒC"              # Threshold temperature for snowfall above glacier [ᵒC]
-    g_cfmax::Vector{T} | "mm ᵒC-1 dt-1" # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
-    g_sifrac::Vector{T} | "dt-1"        # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
-    glacierstore::Vector{T} | "mm"      # Water within the glacier [mm]
-    glacierfrac::Vector{T} | "-"        # Fraction covered by a glacier [-]
-    precipitation::Vector{T}            # Precipitation [mm Δt⁻¹]
-    temperature::Vector{T} | "ᵒC"       # Temperature [ᵒC]
-    potential_evaporation::Vector{T}    # Potential evapotranspiration [mm Δt⁻¹]
-    potsoilevap::Vector{T}              # Potential soil evaporation [mm Δt⁻¹]
-    soilevap::Vector{T}                 # Soil evaporation [mm Δt⁻¹]
-    intevap::Vector{T}                  # Evaporation from interception storage [mm Δt⁻¹]
-    actevap::Vector{T}                  # Total actual evapotranspiration (intevap + soilevap) [mm Δt⁻¹]
-    interceptionstorage::Vector{T} | "mm"  # Actual interception storage [mm]
-    snowwater::Vector{T} | "mm"            # Available free water in snow [mm]
-    snow::Vector{T} | "mm"                 # Snow pack [mm]
-    rainfallplusmelt::Vector{T}            # Snow melt + precipitation as rainfall [mm Δt⁻¹]
-    soilmoisture::Vector{T} | "mm"         # Actual soil moisture [mm]
-    directrunoff::Vector{T}             # Direct runoff to upper zone [mm Δt⁻¹]
-    hbv_seepage::Vector{T}              # Recharge to upper zone [mm Δt⁻¹]
-    in_upperzone::Vector{T}             # Water inflow into upper zone [mm Δt⁻¹]
-    upperzonestorage::Vector{T} | "mm"  # Water content of the upper zone [mm]
-    quickflow::Vector{T}                # Specific runoff (quickflow part) [mm Δt⁻¹]
-    real_quickflow::Vector{T}           # Specific runoff (quickflow), if K upper zone is precalculated [mm Δt⁻¹]
-    percolation::Vector{T}              # Actual percolation to the lower zone [mm Δt⁻¹]
-    capflux::Vector{T}                  # Capillary rise [mm Δt⁻¹]
-    lowerzonestorage::Vector{T} | "mm"  # Water content of the lower zone [mm]
-    baseflow::Vector{T}                 # Specific runoff (baseflow part) per cell [mm Δt⁻¹]
-    runoff::Vector{T}                   # Total specific runoff per cell [mm Δt⁻¹]
-end
-
-
-function update_until_snow(hbv::HBV, config)
-
-    for i = 1:hbv.n
-        precipitation = hbv.precipitation[i] * hbv.pcorr[i]
-        # fraction of precipitation which falls as rain
-        rainfrac = if iszero(hbv.tti[i])
-            Float64(hbv.temperature[i] > hbv.tt[i])
-        else
-            frac = (hbv.temperature[i] - (hbv.tt[i] - hbv.tti[i] / 2.0)) / hbv.tti[i]
-            min(frac, 1.0)
-        end
-        rainfrac = max(rainfrac, 0.0)
-
-        # fraction of precipitation which falls as snow
-        snowfrac = 1.0 - rainfrac
-        # different correction for rainfall and snowfall
-        precipitation =
-            hbv.sfcf[i] * snowfrac * precipitation + hbv.rfcf[i] * rainfrac * precipitation
-
-        # interception
-        interception = min(precipitation, hbv.icf[i] - hbv.interceptionstorage[i])
-        # current interception storage
-        interceptionstorage = hbv.interceptionstorage[i] + interception
-        precipitation = precipitation - interception
-
-        # correction for potential evaporation on wet days
-        potevap =
-            exp(-hbv.epf[i] * precipitation) * hbv.ecorr[i] * hbv.potential_evaporation[i]
-        # correct per landuse
-        potevap = hbv.cevpf[i] * potevap
-
-        # evaporation from interception storage
-        intevap = min(interceptionstorage, potevap)
-        interceptionstorage = interceptionstorage - intevap
-        restevap = max(0.0, potevap - intevap)
-
-        snow, snowwater, snowmelt, rainfallplusmelt, snowfall = snowpack_hbv(
-            hbv.snow[i],
-            hbv.snowwater[i],
-            precipitation,
-            hbv.temperature[i],
-            hbv.tti[i],
-            hbv.tt[i],
-            hbv.ttm[i],
-            hbv.cfmax[i],
-            hbv.whc[i],
-        )
-
-        # update the outputs and states
-        hbv.rainfallplusmelt[i] = rainfallplusmelt
-        hbv.snowwater[i] = snowwater
-        hbv.snow[i] = snow
-        hbv.interceptionstorage[i] = interceptionstorage
-        hbv.potsoilevap[i] = restevap
-        hbv.intevap[i] = intevap
-    end
-end
-
-function update_after_snow(hbv::HBV, config)
-
-    modelglacier = get(config.model, "glacier", false)::Bool
-    set_kquickflow = get(config.model, "set_kquickflow", false)::Bool
-    external_qbase = get(config.model, "external_qbase", false)::Bool
-
-    for i = 1:hbv.n
-        if modelglacier
-            # Run Glacier module and add the snowpack on-top of it.
-            # Estimate the fraction of snow turned into ice (HBV-light).
-            # Estimate glacier melt.
-
-            hbv.snow[i], _, hbv.glacierstore[i], glaciermelt = glacier_hbv(
-                hbv.glacierfrac[i],
-                hbv.glacierstore[i],
-                hbv.snow[i],
-                hbv.temperature[i],
-                hbv.g_tt[i],
-                hbv.g_cfmax[i],
-                hbv.g_sifrac[i],
-                Second(hbv.dt),
-            )
-            # Convert to mm per grid cell and add to snowmelt
-            glaciermelt = glaciermelt * hbv.glacierfrac[i]
-            rainfallplusmelt = hbv.rainfallplusmelt[i] + glaciermelt
-        else
-            rainfallplusmelt = hbv.rainfallplusmelt[i]
-        end
-
-        soilmoisture = hbv.soilmoisture[i] + rainfallplusmelt
-        # if soil is filled to capacity: abundant water runs of directly
-        directrunoff = max(soilmoisture - hbv.fc[i], 0.0)
-        soilmoisture = soilmoisture - directrunoff
-        # net water which infiltrates into soil
-        netinsoil = rainfallplusmelt - directrunoff
-
-        # soil evapotranspiration
-        soilevap =
-            soilmoisture > hbv.threshold[i] ? min(soilmoisture, hbv.potsoilevap[i]) :
-            min(hbv.potsoilevap[i] * (soilmoisture / hbv.threshold[i]))
-        # evaporation from soil moisture storage
-        soilmoisture = soilmoisture - soilevap
-        # sum of evaporation components (IntEvap+SoilEvap)
-        actevap = hbv.intevap[i] + soilevap
-        # runoff water from soil
-        hbv_seepage =
-            pow(min(soilmoisture / hbv.fc[i], 1.0), hbv.betaseepage[i]) * netinsoil
-        soilmoisture = soilmoisture - hbv_seepage
-        # correction for extremely wet periods: soil is filled to capacity
-        back_tosoil = min(hbv.fc[i] - soilmoisture, directrunoff)
-        directrunoff = directrunoff - back_tosoil
-        soilmoisture = soilmoisture + back_tosoil
-        # total water available for runoff
-        in_upperzone = directrunoff + hbv_seepage
-        ### calculations for upper zone ###
-        upperzonestorage = hbv.upperzonestorage[i] + in_upperzone
-        percolation = min(hbv.perc[i], upperzonestorage - in_upperzone / 2.0)
-        upperzonestorage = upperzonestorage - percolation
-        # capillary flux flowing back to soil
-        capflux = hbv.cflux[i] * ((hbv.fc[i] - soilmoisture) / hbv.fc[i])
-        capflux = min(hbv.fc[i] - soilmoisture, capflux)
-        upperzonestorage = upperzonestorage - capflux
-        soilmoisture = soilmoisture + capflux
-
-        real_quickflow = 0.0
-        if set_kquickflow == false
-
-            if percolation < hbv.perc[i]
-                quickflow = 0.0
-            else
-                quickflow = min(
-                    hbv.kquickflow[i] * pow(
-                        (upperzonestorage - min(in_upperzone / 2.0, upperzonestorage)),
-                        1.0 + hbv.alphanl[i],
-                    ),
-                    upperzonestorage,
-                )
-
-                upperzonestorage =
-                    percolation < hbv.perc[i] ? upperzonestorage :
-                    max(upperzonestorage - quickflow, 0.0)
-            end
-        else
-            quickflow = hbv.kquickflow[i] * upperzonestorage
-            real_quickflow = max(0.0, hbv.k0[i] * (upperzonestorage - hbv.suz[i]))
-            upperzonestorage = upperzonestorage - quickflow - real_quickflow
-
-        end
-
-        ### calculations for lower zone ###
-        lowerzonestorage = hbv.lowerzonestorage[i] + percolation
-        # baseflow in mm/timestep
-        baseflow = min(lowerzonestorage, hbv.k4[i] * lowerzonestorage)
-        lowerzonestorage = lowerzonestorage - baseflow
-
-        if external_qbase
-            directrunoffstorage = quickflow + hbv_seepage + real_quickflow
-        else
-            directrunoffstorage = quickflow + baseflow + real_quickflow
-        end
-
-        runoff = max(0.0, directrunoffstorage)
-
-        # update the outputs and states
-        hbv.runoff[i] = runoff
-        hbv.rainfallplusmelt[i] = rainfallplusmelt
-        hbv.soilmoisture[i] = soilmoisture
-        hbv.soilevap[i] = soilevap
-        hbv.actevap[i] = actevap
-        hbv.quickflow[i] = quickflow
-        hbv.real_quickflow[i] = real_quickflow
-        hbv.upperzonestorage[i] = upperzonestorage
-        hbv.lowerzonestorage[i] = lowerzonestorage
-        hbv.baseflow[i] = baseflow
-        hbv.hbv_seepage[i] = hbv_seepage
-        hbv.percolation[i] = percolation
-        hbv.capflux[i] = capflux
-        hbv.in_upperzone[i] = in_upperzone
-    end
-
-end
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
deleted file mode 100644
index aab5cbbd7..000000000
--- a/src/hbv_model.jl
+++ /dev/null
@@ -1,449 +0,0 @@
-"""
-    initialize_hbv_model(config::Config)
-
-Initial part of the SBM model concept. Reads the input settings and data as defined in the
-Config object. Will return a Model that is ready to run.
-"""
-function initialize_hbv_model(config::Config)
-    # unpack the paths to the netCDF files
-    static_path = input_path(config, config.input.path_static)
-
-    reader = prepare_reader(config)
-    clock = Clock(config, reader)
-    dt = clock.dt
-
-    do_reservoirs = get(config.model, "reservoirs", false)::Bool
-    do_lakes = get(config.model, "lakes", false)::Bool
-    do_pits = get(config.model, "pits", false)::Bool
-    set_kquickflow = get(config.model, "set_kquickflow", false)::Bool
-
-    kw_river_tstep = get(config.model, "kw_river_tstep", 0)
-    kw_land_tstep = get(config.model, "kw_land_tstep", 0)
-    kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool
-
-    nc = NCDataset(static_path)
-
-    subcatch_2d = ncread(nc, config, "subcatchment"; optional = false, allow_missing = true)
-    # indices based on catchment
-    inds, rev_inds = active_indices(subcatch_2d, missing)
-    n = length(inds)
-
-    cfmax =
-        ncread(
-            nc,
-            config,
-            "vertical.cfmax";
-            sel = inds,
-            defaults = 3.75653,
-            type = Float,
-        ) .* (dt / basetimestep)
-    tt = ncread(nc, config, "vertical.tt"; sel = inds, defaults = -1.41934, type = Float)
-    tti = ncread(nc, config, "vertical.tti"; sel = inds, defaults = 1.0, type = Float)
-    ttm = ncread(nc, config, "vertical.ttm"; sel = inds, defaults = -1.41934, type = Float)
-    whc = ncread(nc, config, "vertical.whc"; sel = inds, defaults = 0.1, type = Float)
-    # glacier parameters
-    g_tt = ncread(
-        nc,
-        config,
-        "vertical.g_tt";
-        sel = inds,
-        defaults = 0.0,
-        type = Float,
-        fill = 0.0,
-    )
-    g_cfmax =
-        ncread(
-            nc,
-            config,
-            "vertical.g_cfmax";
-            sel = inds,
-            defaults = 3.0,
-            type = Float,
-            fill = 0.0,
-        ) .* (dt / basetimestep)
-    g_sifrac =
-        ncread(
-            nc,
-            config,
-            "vertical.g_sifrac";
-            sel = inds,
-            defaults = 0.001,
-            type = Float,
-            fill = 0.0,
-        ) .* (dt / basetimestep)
-    glacierfrac = ncread(
-        nc,
-        config,
-        "vertical.glacierfrac";
-        sel = inds,
-        defaults = 0.0,
-        type = Float,
-        fill = 0.0,
-    )
-    glacierstore = ncread(
-        nc,
-        config,
-        "vertical.glacierstore";
-        sel = inds,
-        defaults = 5500.0,
-        type = Float,
-        fill = 0.0,
-    )
-    fc = ncread(nc, config, "vertical.fc"; sel = inds, defaults = 260.0, type = Float)
-    betaseepage =
-        ncread(nc, config, "vertical.betaseepage"; sel = inds, defaults = 1.8, type = Float)
-    lp = ncread(nc, config, "vertical.lp"; sel = inds, defaults = 0.53, type = Float)
-    k4 =
-        ncread(nc, config, "vertical.k4"; sel = inds, defaults = 0.02307, type = Float) .*
-        (dt / basetimestep)
-    kquickflow =
-        ncread(
-            nc,
-            config,
-            "vertical.kquickflow";
-            sel = inds,
-            defaults = 0.09880,
-            type = Float,
-        ) .* (dt / basetimestep)
-    suz = ncread(nc, config, "vertical.suz"; sel = inds, defaults = 100.0, type = Float)
-    k0 =
-        ncread(nc, config, "vertical.k0"; sel = inds, defaults = 0.30, type = Float) .*
-        (dt / basetimestep)
-    khq =
-        ncread(nc, config, "vertical.khq"; sel = inds, defaults = 0.09880, type = Float) .*
-        (dt / basetimestep)
-    hq =
-        ncread(nc, config, "vertical.hq"; sel = inds, defaults = 3.27, type = Float) .*
-        (dt / basetimestep)
-    alphanl =
-        ncread(nc, config, "vertical.alphanl"; sel = inds, defaults = 1.1, type = Float)
-    perc =
-        ncread(nc, config, "vertical.perc"; sel = inds, defaults = 0.4, type = Float) .*
-        (dt / basetimestep)
-    cfr = ncread(nc, config, "vertical.cfr"; sel = inds, defaults = 0.05, type = Float)
-    pcorr = ncread(nc, config, "vertical.pcorr"; sel = inds, defaults = 1.0, type = Float)
-    rfcf = ncread(nc, config, "vertical.rfcf"; sel = inds, defaults = 1.0, type = Float)
-    sfcf = ncread(nc, config, "vertical.sfcf"; sel = inds, defaults = 1.0, type = Float)
-    cflux =
-        ncread(nc, config, "vertical.cflux"; sel = inds, defaults = 2.0, type = Float) .*
-        (dt / basetimestep)
-    icf = ncread(nc, config, "vertical.icf"; sel = inds, defaults = 2.0, type = Float)
-    cevpf = ncread(nc, config, "vertical.cevpf"; sel = inds, defaults = 1.0, type = Float)
-    epf = ncread(nc, config, "vertical.epf"; sel = inds, defaults = 1.0, type = Float)
-    ecorr = ncread(nc, config, "vertical.ecorr"; sel = inds, defaults = 1.0, type = Float)
-
-    # read x, y coordinates and calculate cell length [m]
-    y_nc = read_y_axis(nc)
-    x_nc = read_x_axis(nc)
-    y = permutedims(repeat(y_nc, outer = (1, length(x_nc))))[inds]
-    cellength = abs(mean(diff(x_nc)))
-
-
-    sizeinmetres = get(config.model, "sizeinmetres", false)::Bool
-    xl, yl = cell_lengths(y, cellength, sizeinmetres)
-
-    threshold = fc .* lp
-
-    hbv = HBV{Float}(
-        dt = Float(tosecond(dt)),
-        n = n,
-        fc = fc,
-        betaseepage = betaseepage,
-        lp = lp,
-        threshold = threshold,
-        k4 = k4,
-        kquickflow = set_kquickflow ? kquickflow :
-                     pow.(khq, 1.0 .+ alphanl) .* pow.(hq, -alphanl),
-        suz = suz,
-        k0 = k0,
-        khq = khq,
-        hq = hq,
-        alphanl = alphanl,
-        perc = perc,
-        cfr = cfr,
-        pcorr = pcorr,
-        rfcf = rfcf,
-        sfcf = sfcf,
-        cflux = cflux,
-        icf = icf,
-        cevpf = cevpf,
-        epf = epf,
-        ecorr = ecorr,
-        tti = tti,
-        tt = tt,
-        ttm = ttm,
-        cfmax = cfmax,
-        whc = whc,
-        # glacier parameters
-        g_tt = g_tt,
-        g_sifrac = g_sifrac,
-        g_cfmax = g_cfmax,
-        glacierstore = glacierstore,
-        glacierfrac = glacierfrac,
-        # default (cold) states:
-        interceptionstorage = zeros(Float, n),
-        snow = zeros(Float, n),
-        snowwater = zeros(Float, n),
-        soilmoisture = copy(fc),
-        upperzonestorage = 0.2 .* fc,
-        lowerzonestorage = 1.0 ./ (3.0 .* k4),
-        # variables:
-        precipitation = fill(mv, n),
-        temperature = fill(mv, n),
-        potential_evaporation = fill(mv, n),
-        potsoilevap = fill(mv, n),
-        soilevap = fill(mv, n),
-        intevap = fill(mv, n),
-        actevap = fill(mv, n),
-        rainfallplusmelt = fill(mv, n),
-        directrunoff = fill(mv, n),
-        hbv_seepage = fill(mv, n),
-        in_upperzone = fill(mv, n),
-        quickflow = fill(mv, n),
-        real_quickflow = fill(mv, n),
-        percolation = fill(mv, n),
-        capflux = fill(mv, n),
-        baseflow = fill(mv, n),
-        runoff = fill(mv, n),
-    )
-
-    modelsize_2d = size(subcatch_2d)
-    river_2d =
-        ncread(nc, config, "river_location"; optional = false, type = Bool, fill = false)
-    river = river_2d[inds]
-    riverwidth_2d =
-        ncread(nc, config, "lateral.river.width"; optional = false, type = Float, fill = 0)
-    riverwidth = riverwidth_2d[inds]
-    riverlength_2d =
-        ncread(nc, config, "lateral.river.length"; optional = false, type = Float, fill = 0)
-    riverlength = riverlength_2d[inds]
-
-    inds_riv, rev_inds_riv = active_indices(river_2d, 0)
-    nriv = length(inds_riv)
-
-    # reservoirs
-    pits = zeros(Bool, modelsize_2d)
-    if do_reservoirs
-        reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
-    else
-        reservoir = ()
-        reservoirs = nothing
-        resindex = fill(0, nriv)
-    end
-
-    # lakes
-    if do_lakes
-        lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
-    else
-        lake = ()
-        lakes = nothing
-        lakeindex = fill(0, nriv)
-    end
-
-    ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
-    ldd = ldd_2d[inds]
-    if do_pits
-        pits_2d = ncread(nc, config, "pits"; optional = false, type = Bool, fill = false)
-        ldd = set_pit_ldd(pits_2d, ldd, inds)
-    end
-
-    landslope =
-        ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(landslope, 0.00001, Inf)
-
-    dl = map(detdrainlength, ldd, xl, yl)
-    dw = (xl .* yl) ./ dl
-    olf = initialize_surfaceflow_land(
-        nc,
-        config,
-        inds;
-        sl = landslope,
-        dl = dl,
-        width = map(det_surfacewidth, dw, riverwidth, river),
-        iterate = kinwave_it,
-        tstep = kw_land_tstep,
-        dt = dt,
-    )
-
-    graph = flowgraph(ldd, inds, pcr_dir)
-
-    riverlength = riverlength_2d[inds_riv]
-    riverwidth = riverwidth_2d[inds_riv]
-    minimum(riverlength) > 0 || error("river length must be positive on river cells")
-    minimum(riverwidth) > 0 || error("river width must be positive on river cells")
-
-    ldd_riv = ldd_2d[inds_riv]
-    if do_pits
-        ldd_riv = set_pit_ldd(pits_2d, ldd_riv, inds_riv)
-    end
-    graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)
-
-    # the indices of the river cells in the land(+river) cell vector
-    index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
-
-    rf = initialize_surfaceflow_river(
-        nc,
-        config,
-        inds_riv;
-        dl = riverlength,
-        width = riverwidth,
-        reservoir_index = resindex,
-        reservoir = reservoirs,
-        lake_index = lakeindex,
-        lake = lakes,
-        iterate = kinwave_it,
-        tstep = kw_river_tstep,
-        dt = dt,
-    )
-
-    # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
-    # subdomain is equal to the complete domain
-    toposort = topological_sort_by_dfs(graph)
-    toposort_riv = topological_sort_by_dfs(graph_riv)
-    index_pit_land = findall(x -> x == 5, ldd)
-    index_pit_river = findall(x -> x == 5, ldd_riv)
-    streamorder = stream_order(graph, toposort)
-    min_streamorder_land = get(config.model, "min_streamorder_land", 5)
-    subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
-        graph,
-        toposort,
-        index_pit_land,
-        streamorder,
-        min_streamorder_land,
-    )
-    min_streamorder_river = get(config.model, "min_streamorder_river", 6)
-    subriv_order, indices_subriv, topo_subriv = kinwave_set_subdomains(
-        graph_riv,
-        toposort_riv,
-        index_pit_river,
-        streamorder[index_river],
-        min_streamorder_river,
-    )
-
-    modelmap = (vertical = hbv, lateral = (land = olf, river = rf))
-    indices_reverse = (
-        land = rev_inds,
-        river = rev_inds_riv,
-        reservoir = isempty(reservoir) ? nothing : reservoir.reverse_indices,
-        lake = isempty(lake) ? nothing : lake.reverse_indices,
-    )
-    writer = prepare_writer(config, modelmap, indices_reverse, x_nc, y_nc, nc)
-    close(nc)
-
-    # for each domain save:
-    # - the directed acyclic graph (graph),
-    # - the traversion order (order),
-    # - upstream_nodes,
-    # - subdomains for the kinematic wave domains for parallel execution (execution order of
-    #   subbasins (subdomain_order), traversion order per subbasin (topo_subdomain) and
-    #   Vector indices per subbasin matching the traversion order of the complete domain
-    #   (indices_subdomain))
-    # - the indices that map it back to the two dimensional grid (indices)
-
-    # for the land domain the x and y length [m] of the grid cells are stored
-    # for reservoirs and lakes indices information is available from the initialization
-    # functions
-    land = (
-        graph = graph,
-        upstream_nodes = filter_upsteam_nodes(graph, pits[inds]),
-        subdomain_order = subbas_order,
-        topo_subdomain = topo_subbas,
-        indices_subdomain = indices_subbas,
-        order = toposort,
-        indices = inds,
-        reverse_indices = rev_inds,
-        xl = xl,
-        yl = yl,
-    )
-    river = (
-        graph = graph_riv,
-        upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
-        subdomain_order = subriv_order,
-        topo_subdomain = topo_subriv,
-        indices_subdomain = indices_subriv,
-        order = toposort_riv,
-        indices = inds_riv,
-        reverse_indices = rev_inds_riv,
-    )
-
-    model = Model(
-        config,
-        (; land, river, reservoir, lake, index_river, frac_toriver),
-        (subsurface = nothing, land = olf, river = rf),
-        hbv,
-        clock,
-        reader,
-        writer,
-        HbvModel(),
-    )
-
-    model = set_states(model)
-
-    return model
-end
-
-function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
-    @unpack lateral, vertical, network, clock, config = model
-
-    inds_riv = network.index_river
-
-    # vertical hbv concept is updated until snow state, after that (optional)
-    # snow transport is possible
-    update_until_snow(vertical, config)
-
-    # lateral snow transport
-    if get(config.model, "masswasting", false)::Bool
-        lateral_snow_transport!(
-            vertical.snow,
-            vertical.snowwater,
-            lateral.land.sl,
-            network.land,
-        )
-    end
-
-    # update vertical hbv concept
-    update_after_snow(vertical, config)
-
-    surface_routing(model)
-
-    return model
-end
-
-function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
-    @unpack lateral, config = model
-
-    reinit = get(config.model, "reinit", true)::Bool
-    do_lakes = get(config.model, "lakes", false)::Bool
-    # read and set states in model object if reinit=true
-    if reinit == false
-        instate_path = input_path(config, config.state.path_input)
-        @info "Set initial conditions from state file `$instate_path`."
-        set_states(instate_path, model; type = Float)
-        # update kinematic wave volume for river and land domain
-        @unpack lateral = model
-        # makes sure land cells with zero flow width are set to zero q and h
-        for i in eachindex(lateral.land.width)
-            if lateral.land.width[i] <= 0.0
-                lateral.land.q[i] = 0.0
-                lateral.land.h[i] = 0.0
-            end
-        end
-        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
-        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
-
-        if do_lakes
-            # storage must be re-initialized after loading the state with the current
-            # waterlevel otherwise the storage will be based on the initial water level
-            lakes = lateral.river.lake
-            lakes.storage .=
-                initialize_storage(lakes.storfunc, lakes.area, lakes.waterlevel, lakes.sh)
-        end
-    else
-        @info "Set initial conditions from default values."
-    end
-
-    return model
-end
diff --git a/src/io.jl b/src/io.jl
index 16d89baed..60c8c0dae 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -464,12 +464,6 @@ function set_extradim_netcdf(
     if extra_dim.name == "layer"
         attributes =
             ["long_name" => "layer_index", "standard_name" => "layer_index", "axis" => "Z"]
-    elseif extra_dim.name == "classes"
-        attributes = [
-            "long_name" => extra_dim.name,
-            "standard_name" => extra_dim.name,
-            "axis" => "Z",
-        ]
     end
     defVar(ds, extra_dim.name, extra_dim.value, (extra_dim.name,), attrib = attributes)
     return nothing
@@ -1437,7 +1431,7 @@ function internal_dim_name(name::Symbol)
         return :x
     elseif name in (:y, :lat, :latitude)
         return :y
-    elseif name in (:time, :layer, :flood_depth, :classes)
+    elseif name in (:time, :layer, :flood_depth)
         return name
     elseif startswith(string(name), "time")
         return :time
@@ -1522,8 +1516,6 @@ function permute_data(data, dim_names)
     @assert ndims(data) == length(dim_names)
     if :layer in dim_names
         desired_order = (:x, :y, :layer)
-    elseif :classes in dim_names
-        desired_order = (:x, :y, :classes)
     elseif :flood_depth in dim_names
         desired_order = (:x, :y, :flood_depth)
     else
@@ -1567,12 +1559,6 @@ function reverse_data!(data, dims_increasing)
             y = dims_increasing.y,
             flood_depth = dims_increasing.flood_depth,
         )
-    elseif length(dims_increasing) == 3 && haskey(dims_increasing, :classes)
-        dims_increasing_ordered = (
-            x = dims_increasing.x,
-            y = dims_increasing.y,
-            classes = dims_increasing.classes,
-        )
     else
         error("Unsupported number of dimensions")
     end
@@ -1628,29 +1614,24 @@ function read_y_axis(ds::CFDataset)::Vector{Float64}
     error("no y axis found in $(path(ds))")
 end
 
-"Get `index` for dimension name `layer` or `classes` based on `model`"
+"Get `index` for dimension name `layer` based on `model`"
 function get_index_dimension(var, model)::Int
     @unpack vertical = model
     if haskey(var, "layer")
         inds = collect(1:vertical.maxlayers)
         index = inds[var["layer"]]
-    elseif haskey(var, "class")
-        index = findfirst(x -> x == var["class"], vertical.classes)
     else
         error("Unrecognized or missing dimension name to index $(var)")
     end
     return index
 end
 
-"Get `index` for dimension name `layer` or `classes` based on `config` (TOML file)"
+"Get `index` for dimension name `layer` based on `config` (TOML file)"
 function get_index_dimension(var, config::Config, dim_value)::Int
     if haskey(var, "layer")
         v = get(config.model, "thicknesslayers", Float[])
         inds = collect(1:length(v)+1)
         index = inds[dim_value]
-    elseif haskey(var, "class")
-        classes = get(config.model, "classes", "")
-        index = findfirst(x -> x == dim_value, classes)
     else
         error("Unrecognized or missing dimension name to index $(var)")
     end
diff --git a/src/states.jl b/src/states.jl
index 9af01cc0e..2c5eb4d70 100644
--- a/src/states.jl
+++ b/src/states.jl
@@ -29,26 +29,6 @@ function get_vertical_states(model_type::AbstractString; snow = false, glacier =
         else
             vertical_states = (:satwaterdepth, :ustorelayerdepth, :canopystorage)
         end
-    elseif model_type == "hbv"
-        vertical_states = (
-            :soilmoisture,
-            :snow,
-            :snowwater,
-            :upperzonestorage,
-            :lowerzonestorage,
-            :interceptionstorage,
-        )
-    elseif model_type == "flextopo"
-        vertical_states = (
-            :snow,
-            :snowwater,
-            :interceptionstorage,
-            :hortonpondingstorage,
-            :hortonrunoffstorage,
-            :rootzonestorage,
-            :faststorage,
-            :slowstorage,
-        )
     elseif model_type == "sediment"
         vertical_states = ()
     else
diff --git a/src/utils.jl b/src/utils.jl
index 57563b7b6..c91442e97 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -154,8 +154,6 @@ function set_states(instate_path, model; type = nothing, dimname = nothing)
             if dims == 4
                 if dimname == :layer
                     dimensions = (x = :, y = :, layer = :, time = 1)
-                elseif dimname == :classes
-                    dimensions = (x = :, y = :, classes = :, time = 1)
                 else
                     error("Unrecognized dimension name $dimname")
                 end
@@ -249,8 +247,6 @@ function ncread(
     # first timestep), that is later updated with the `update_cyclic!` function.
     if isnothing(dimname)
         dim_sel = (x = :, y = :, time = 1)
-    elseif dimname == :classes
-        dim_sel = (x = :, y = :, classes = :, time = 1)
     elseif dimname == :layer
         dim_sel = (x = :, y = :, layer = :, time = 1)
     elseif dimname == :flood_depth
@@ -810,7 +806,8 @@ end
 function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
     for i in eachindex(ssf.ssf)
         ssf.kh[i] = kh_layered_profile(sbm, ssf.khfrac[i], i, ksat_profile)
-        ssf.ssf[i] = ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.slope[i] * ssf.dw[i]
+        ssf.ssf[i] =
+            ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.slope[i] * ssf.dw[i]
         kh_max = 0.0
         for j = 1:sbm.nlayers[i]
             if j <= sbm.nlayers_kv[i]


Test patch:
diff --git a/build/create_binaries/download_test_data.jl b/build/create_binaries/download_test_data.jl
index d61310f3b..8ad77e30f 100644
--- a/build/create_binaries/download_test_data.jl
+++ b/build/create_binaries/download_test_data.jl
@@ -20,12 +20,8 @@ end
 
 staticmaps_rhine_path = testdata(v"0.1", "staticmaps.nc", "staticmaps-rhine.nc")
 staticmaps_moselle_path =
-    testdata(v"0.2.8", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
-staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn.nc")
-staticmaps_meuse_path =
-    testdata(v"0.2.8", "staticmaps_flex_meuse.nc", "staticmaps_flex_meuse.nc")
+    testdata(v"0.2.9", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
 forcing_moselle_path = testdata(v"0.2.6", "forcing-moselle.nc", "forcing-moselle.nc")
-forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
 forcing_moselle_sed_path =
     testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
 staticmaps_moselle_sed_path =
@@ -43,7 +39,6 @@ forcing_sbm_gw_path = testdata(
     "forcing-sbm-groundwater-part2.nc",
     "forcing-sbm-groundwater-part2.nc",
 )
-forcing_meuse_path = testdata(v"0.2.8", "forcing_meuse.nc", "forcing_meuse.nc")
 staticmaps_sbm_gw_path =
     testdata(v"0.2.3", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
 instates_sbm_gw_path =
diff --git a/test/flextopo_config.toml b/test/flextopo_config.toml
deleted file mode 100644
index aba6659af..000000000
--- a/test/flextopo_config.toml
+++ /dev/null
@@ -1,523 +0,0 @@
-casename = "wflow_meuse"
-calendar = "proleptic_gregorian"
-starttime = "2009-12-31T00:00:00"
-endtime = "2010-07-01T00:00:00"
-time_units = "days since 1900-01-01 00:00:00"
-timestepsecs = 86400
-dir_input = "data/input"
-dir_output = "data/output"
-loglevel = "info"
-
-[state]
-path_input = "instates.nc"
-path_output = "outstates-meuse.nc"
-
-[input]
-path_forcing = "forcing_meuse.nc"
-path_static = "staticmaps_flex_meuse.nc"
-gauges = "wflow_gauges"
-ldd = "wflow_ldd"
-river_location = "wflow_river"
-subcatchment = "wflow_subcatch"
-forcing = [ "vertical.precipitation", "vertical.temperature", "vertical.potential_evaporation",]
-#cyclic = [ "vertical.leaf_area_index",]
-gauges_grdc = "wflow_gauges_grdc"
-gauges_S01 = "wflow_gauges_S01"
-gauges_S02 = "wflow_gauges_S02"
-gauges_S03 = "wflow_gauges_S03"
-gauges_S04 = "wflow_gauges_S04"
-gauges_S05 = "wflow_gauges_S05"
-gauges_S06 = "wflow_gauges_S06"
-gauges_Sall = "wflow_gauges_Sall"
-
-sub_S01 = "wflow_subcatch_S01"
-sub_S02 = "wflow_subcatch_S02"
-sub_S03 = "wflow_subcatch_S03"
-sub_S04 = "wflow_subcatch_S04"
-sub_S05 = "wflow_subcatch_S05"
-sub_S06 = "wflow_subcatch_S06"
-sub_Sall = "wflow_subcatch_Sall"
-
-[model]
-type = "flextopo"
-masswasting = true
-snow = true
-reinit = true
-reservoirs = false
-lakes = false
-glacier = false
-kin_wave_iteration = true
-kw_river_tstep = 900
-kw_land_tstep = 3600
-classes = ["h", "p", "w"]
-
-
-select_snow = ["common_snow_hbv"]
-select_interception = ["interception_overflow", "interception_overflow", "interception_overflow"]
-select_hortonponding = ["hortonponding_no_storage", "hortonponding_no_storage", "hortonponding_no_storage"]
-select_hortonrunoff = ["hortonrunoff_no_storage", "hortonrunoff_no_storage", "hortonrunoff_no_storage"]
-select_rootzone = ["rootzone_storage", "rootzone_storage", "rootzone_storage"]
-select_fast = ["fast_storage", "fast_storage", "fast_storage"]
-select_slow = ["common_slow_storage"]
-
-[input.vertical]
-altitude = "wflow_dem"
-potential_evaporation = "PET"
-precipitation = "P"
-temperature = "TEMP"
-tt = "tth"
-ttm = "tmh"
-cfmax = "fmh"
-whc = "WHC"
-#alfa = "alfa"
-beta = "beta"
-cap = "cap"
-#ds = "d"
-fdec = "decf"
-fmax = "fmax"
-imax = "imax"
-#kf = "kf"
-#ks = "ksh"
-#lp = "lp"
-perc = "perc"
-srmax = "sumax"
-hrufrac = "hrufrac_lu"
-
-[input.vertical.lp]
-netcdf.variable.name = "lp"
-scale = [0.4, 0.4, 0.4]
-offset = [0.0, 0.0, 0.0]
-class = ["h", "p", "w"]
-
-[input.vertical.ds]
-netcdf.variable.name = "d"
-scale = [1.2, 1.2, 1.0]
-offset = [0.0, 0.0, 0.0]
-class = ["h", "p", "w"]
-
-[input.vertical.kf]
-netcdf.variable.name = "kf"
-scale = [1.0, 3.0, 3.0]
-offset = [0.0, 0.0, 0.0]
-class = ["h", "p", "w"]
-
-[input.vertical.alfa]
-netcdf.variable.name = "alfa"
-scale = 1.3
-offset = 0
-class = "p"
-
-[input.vertical.ks]
-netcdf.variable.name = "ksh"
-scale = 0.5
-offset = 0.0
-
-
-[input.lateral.river]
-length = "wflow_riverlength"
-#n = "N_River"
-n = "n_river_uniform"
-slope = "RiverSlope"
-width = "wflow_riverwidth"
-bankfull_depth = "RiverDepth"
-
-[input.lateral.land]
-#n = "N"
-n = "n_uniform"
-slope = "Slope"
-
-[input.lateral.river.reservoir]
-area = "ResSimpleArea"
-areas = "wflow_reservoirareas"
-demand = "ResDemand"
-locs = "wflow_reservoirlocs"
-maxrelease = "ResMaxRelease"
-maxvolume = "ResMaxVolume"
-targetfullfrac = "ResTargetFullFrac"
-targetminfrac = "ResTargetMinFrac"
-
-[state.lateral.river.reservoir]
-#volume = "volume_reservoir"
-
-[state.vertical]
-snow = "snow"
-snowwater = "snowwater"
-interceptionstorage = "interceptionstorage"
-hortonpondingstorage = "hortonpondingstorage"
-hortonrunoffstorage = "hortonrunoffstorage"
-rootzonestorage = "rootzonestorage"
-faststorage = "faststorage"
-slowstorage = "slowstorage"
-
-
-[state.lateral.river]
-q = "q_river"
-h = "h_river"
-h_av = "h_av_river"
-
-
-[state.lateral.land]
-q = "q_land"
-h = "h_land"
-h_av = "h_av_land"
-
-
-[output]
-path = "output-flex-meuse.nc"
-
-[output.vertical]
-#precipitation = "prec"
-#temperature = "temp"
-#potential_evaporation = "pet"
-faststorage = "faststorage"
-
-
-[output.lateral.river]
-q_av = "q_river"
-#h = "h_river"
-
-
-[output.lateral.land]
-#q = "q_land"
-#h = "h_land"
-
-
-[csv]
-path = "output-flex-meuse.csv"
-
-[[csv.column]]
-header = "Q"
-map = "gauges"
-parameter = "lateral.river.q_av"
-
-[[csv.column]]
-header = "Q"
-map = "gauges_grdc"
-parameter = "lateral.river.q_av"
-
-
-[[csv.column]]
-header = "Q"
-map = "gauges_Sall"
-parameter = "lateral.river.q_av"
-
-[[csv.column]]
-header = "H"
-map = "gauges_Sall"
-parameter = "lateral.river.h_av"
-
-[[csv.column]]
-header = "P"
-map = "sub_S06"
-parameter = "vertical.precipitation"
-reducer = "mean"
-
-[[csv.column]]
-header = "Ep"
-map = "sub_S06"
-parameter = "vertical.potential_evaporation"
-reducer = "mean"
-
-[[csv.column]]
-header = "T"
-map = "sub_S06"
-parameter = "vertical.temperature"
-reducer = "mean"
-
-[[csv.column]]
-header = "Ea"
-map = "sub_S06"
-parameter = "vertical.actevap_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Ei"
-map = "sub_S06"
-parameter = "vertical.intevap_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Er"
-map = "sub_S06"
-parameter = "vertical.rootevap_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Eh"
-map = "sub_S06"
-parameter = "vertical.hortonevap_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sw"
-map = "sub_S06"
-parameter = "vertical.snow"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sww"
-map = "sub_S06"
-parameter = "vertical.snowwater"
-reducer = "mean"
-
-[[csv.column]]
-header = "Si"
-map = "sub_S06"
-parameter = "vertical.interceptionstorage_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sr"
-map = "sub_S06"
-parameter = "vertical.srootzone_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sh"
-map = "sub_S06"
-parameter = "vertical.hortonpondingstorage_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Shf"
-map = "sub_S06"
-parameter = "vertical.hortonrunoffstorage_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sf"
-map = "sub_S06"
-parameter = "vertical.faststorage_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "Ss"
-map = "sub_S06"
-parameter = "vertical.slowstorage"
-reducer = "mean"
-
-[[csv.column]]
-header = "Sr_over_srmax_p"
-map = "sub_S06"
-parameter = "vertical.srootzone_over_srmax"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "Sr_over_srmax_m"
-map = "sub_S06"
-parameter = "vertical.srootzone_over_srmax_m"
-reducer = "mean"
-
-
-
-
-[[csv.column]]
-header = "Qftotal"
-map = "sub_S06"
-parameter = "vertical.qfast_tot"
-reducer = "mean"
-
-[[csv.column]]
-header = "Qh"
-map = "sub_S06"
-parameter = "vertical.qhortonpond"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "Qhf"
-map = "sub_S06"
-parameter = "vertical.qhortonrun"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "Qs"
-map = "sub_S01"
-parameter = "vertical.qslow"
-reducer = "mean"
-
-[[csv.column]]
-header = "Qs"
-map = "sub_S06"
-parameter = "vertical.qslow"
-reducer = "mean"
-
-[[csv.column]]
-header = "Qcap"
-map = "sub_S06"
-parameter = "vertical.qcapillary_m"
-reducer = "mean"
-
-[[csv.column]]
-header = "QfP"
-map = "sub_S06"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "QfW"
-map = "sub_S06"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "w"
-
-[[csv.column]]
-header = "QfH"
-map = "sub_S06"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "h"
-
-
-[[csv.column]]
-header = "percentageP"
-map = "sub_S01"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "percentageW"
-map = "sub_S01"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "w"
-
-[[csv.column]]
-header = "percentageH"
-map = "sub_S01"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "h"
-
-[[csv.column]]
-header = "percentageP"
-map = "sub_S06"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "percentageW"
-map = "sub_S06"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "w"
-
-[[csv.column]]
-header = "percentageH"
-map = "sub_S06"
-parameter = "vertical.hrufrac"
-reducer = "mean"
-class = "h"
-
-
-
-
-[[csv.column]]
-header = "QfP"
-map = "sub_S01"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "QfW"
-map = "sub_S01"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "w"
-
-[[csv.column]]
-header = "QfH"
-map = "sub_S01"
-parameter = "vertical.qfast"
-reducer = "mean"
-class = "h"
-
-
-[[csv.column]]
-header = "kfP"
-map = "sub_Sall"
-parameter = "vertical.kf"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "kfW"
-map = "sub_Sall"
-parameter = "vertical.kf"
-reducer = "mean"
-class = "w"
-
-[[csv.column]]
-header = "kfH"
-map = "sub_Sall"
-parameter = "vertical.kf"
-reducer = "mean"
-class = "h"
-
-
-
-[[csv.column]]
-header = "wbtot"
-map = "sub_S06"
-parameter = "vertical.wb_tot"
-reducer = "mean"
-
-[[csv.column]]
-header = "wbSi_p"
-map = "sub_S06"
-parameter = "vertical.wb_interception"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "wbShf_p"
-map = "sub_S06"
-parameter = "vertical.wb_hortonrunoff"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "wbSh_p"
-map = "sub_S06"
-parameter = "vertical.wb_hortonponding"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "wbSr_p"
-map = "sub_S06"
-parameter = "vertical.wb_rootzone"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "wbSf_p"
-map = "sub_S06"
-parameter = "vertical.wb_fast"
-reducer = "mean"
-class = "p"
-
-[[csv.column]]
-header = "wbSs"
-map = "sub_S06"
-parameter = "vertical.wb_slow"
-reducer = "mean"
-
-[[csv.column]]
-header = "to_river_land"
-map = "gauges_S06"
-parameter = "lateral.land.to_river"
-
-[[csv.column]]
-header = "land_inwater"
-map = "sub_S06"
-parameter = "lateral.land.inwater"
-reducer = "mean"
-
-
diff --git a/test/hbv_config.toml b/test/hbv_config.toml
deleted file mode 100644
index 1a6a584b5..000000000
--- a/test/hbv_config.toml
+++ /dev/null
@@ -1,137 +0,0 @@
-# This is a TOML configuration file for Wflow.
-# Relative file paths are interpreted as being relative to this TOML file.
-# Wflow documentation https://deltares.github.io/Wflow.jl/dev/
-# TOML documentation: https://github.com/toml-lang/toml
-
-calendar = "proleptic_gregorian"
-endtime = 2000-02-01T00:00:00
-starttime = 1999-12-31T00:00:00
-time_units = "days since 1900-01-01 00:00:00"
-timestepsecs = 86400
-
-[state]
-path_input = "data/input/instates-lahn.nc"
-path_output = "data/output/outstates-lahn.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]
-interceptionstorage = "interceptionstorage"
-lowerzonestorage = "lowerzonestorage"
-snow = "snow"
-snowwater = "snowwater"
-soilmoisture = "soilmoisture"
-upperzonestorage = "upperzonestorage"
-
-[state.lateral.river]
-h = "h_river"
-q = "q_river"
-
-[state.lateral.land]
-h = "h_land"
-q = "q_land"
-
-[input]
-path_forcing = "data/input/forcing-lahn.nc"
-path_static = "data/input/staticmaps-lahn.nc"
-
-# these are not directly part of the model
-gauges = "wflow_gauges"
-ldd = "wflow_ldd"
-river_location = "wflow_river"
-subcatchment = "wflow_subcatch"
-
-# specify the internal IDs of the parameters which vary over time
-# the external name mapping needs to be below together with the other mappings
-forcing = [
-  "vertical.precipitation",
-  "vertical.temperature",
-  "vertical.potential_evaporation",
-]
-
-[input.vertical]
-alphanl = "AlphaNL"
-altitude = "wflow_dem"
-betaseepage = "BetaSeepage"
-cevpf = "CEVPF"
-cflux = "Cflux"
-cfmax = "Cfmax"
-cfr = "CFR"
-epf = "EPF"
-fc = "FC"
-hq = "HQ"
-icf = "ICF"
-k4 = "K4"
-khq = "KHQ"
-lp = "LP"
-perc = "PERC"
-potential_evaporation = "PET"
-precipitation = "P"
-sfcf = "SFCF"
-temperature = "TEMP"
-tt = "TT"
-tti = "TTI"
-ttm = "TTM"
-whc = "WHC"
-
-[input.lateral.river]
-length = "wflow_riverlength"
-n = "N_River"
-slope = "RiverSlope"
-width = "wflow_riverwidth"
-
-[input.lateral.land]
-n = "N"
-slope = "Slope"
-
-[model]
-kin_wave_iteration = true
-masswasting = true
-reinit = true
-snow = true
-type = "hbv"
-
-[output]
-path = "data/output/output_lahn.nc"
-
-[output.vertical]
-lowerzonestorage = "lowerzonestorage"
-snow = "snow"
-snowwater = "snowwater"
-soilmoisture = "soilmoisture"
-upperzonestorage = "upperzonestorage"
-
-[output.lateral.river]
-q = "q"
-
-[csv]
-path = "data/output/output_lahn.csv"
-
-[[csv.column]]
-header = "Q"
-parameter = "lateral.river.q"
-reducer = "maximum"
-
-[[csv.column]]
-coordinate.x = 8.279
-coordinate.y = 50.534
-header = "temp_bycoord"
-parameter = "vertical.temperature"
-
-[[csv.column]]
-header = "temp_byindex"
-index.x = 88
-index.y = 95
-parameter = "vertical.temperature"
-
-[[csv.column]]
-header = "Q"
-map = "gauges"
-parameter = "lateral.river.q"
-
-[[csv.column]]
-header = "perc"
-map = "subcatchment"
-parameter = "vertical.perc"
-reducer = "mean"
diff --git a/test/io.jl b/test/io.jl
index f692efed3..73d681ee6 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -47,14 +47,6 @@ config = Wflow.Config(tomlpath)
           joinpath(@__DIR__, "data", "input", "instates-moselle.nc")
     @test Wflow.output_path(config, config.state.path_output) ==
           joinpath(@__DIR__, "data", "output", "outstates-moselle.nc")
-    # hbv_config doesn't use dir_input and dir_output
-    hbv_config = Wflow.Config(joinpath(@__DIR__, "hbv_config.toml"))
-    @test !haskey(hbv_config, "dir_input")
-    @test !haskey(hbv_config, "dir_output")
-    @test Wflow.input_path(hbv_config, hbv_config.state.path_input) ==
-          joinpath(@__DIR__, "data", "input", "instates-lahn.nc")
-    @test Wflow.output_path(hbv_config, hbv_config.state.path_output) ==
-          joinpath(@__DIR__, "data", "output", "outstates-lahn.nc")
 end
 
 @testset "Clock constructor" begin
diff --git a/test/run_flextopo.jl b/test/run_flextopo.jl
deleted file mode 100644
index cfb4bbffb..000000000
--- a/test/run_flextopo.jl
+++ /dev/null
@@ -1,116 +0,0 @@
-
-tomlpath = joinpath(@__DIR__, "flextopo_config.toml")
-config = Wflow.Config(tomlpath)
-
-model = Wflow.initialize_flextopo_model(config)
-@unpack network = model
-
-model = Wflow.run_timestep(model)
-
-# test if the first timestep was written to the CSV file
-flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
-@testset "CSV output" begin
-    row = csv_first_row(model.writer.csv_path)
-
-    @test row.time == DateTime("2010-01-01T00:00:00")
-    @test row.Q_16 ≈ 0.0000370817920883859f0
-    @test row.percentageH_16 ≈ 0.340213f0
-    @test row.QfP_1011 ≈ 0.0f0
-    @test row.kfW_10 ≈ 0.136334478855133f0
-    @test row.kfH_10 ≈ 0.0454448275268077f0
-    @test row.Qs_503 ≈ 0.1575031550601124f0
-    @test row.Si_16 ≈ 0.0f0
-    @test row.Ss_16 ≈ 29.89794354323524f0
-    @test row.Ea_16 ≈ 0.0110627892408438f0
-end
-
-@testset "first timestep" begin
-    flextopo = model.vertical
-    @test flextopo.tt[3500] ≈ 1.3f0
-    @test model.clock.iteration == 1
-    @test flextopo.rootzonestorage[3500] ≈
-          [147.11238663084805f0, 79.08369375691255f0, 79.23637697443984f0]
-    @test flextopo.runoff[3500] ≈ 0.19008129369467497f0
-    @test flextopo.rootevap[3500] ≈
-          [0.009225917980074883f0, 0.009225917980074883f0, 0.009225917980074883f0]
-    @test flextopo.snow[3500] == 0.0
-end
-
-# run the second timestep
-model = Wflow.run_timestep(model)
-
-@testset "second timestep" begin
-    flextopo = model.vertical
-    @test flextopo.rootzonestorage[3500] ≈
-          [146.73222385533154f0, 78.55190222246516f0, 78.85739340140256f0]
-    @test flextopo.runoff[3500] ≈ 0.18887692333975817f0
-    @test flextopo.rootevap[3500] ≈
-          [0.38016277551651f0, 0.38016277551651f0, 0.38016277551651f0]
-    @test flextopo.snow[3500] == 0.0f0
-end
-
-@testset "overland domain" begin
-    q = model.lateral.land.q_av
-    land = model.lateral.land
-    @test sum(q) ≈ 202.1162497378859f0
-    @test q[10354] ≈ 0.001624108000162103f0
-    @test land.volume[10354] ≈ 66.41398729880544f0
-    @test land.inwater[10354] ≈ 0.001658182220814705f0
-    @test q[network.land.order[end]] ≈ 0.0033337667005815565f0
-end
-
-@testset "river flow" begin
-    q = model.lateral.river.q_av
-    @test sum(q) ≈ 108.45732176546284f0
-    @test q[651] ≈ 0.014755293483684526f0
-    @test q[1056] ≈ 0.0030933658056867823f0
-    @test q[network.river.order[end]] ≈ 0.0005541137960384742f0
-end
-
-#change the config to use other functions for the storages for several classes:
-config = Wflow.Config(tomlpath)
-config["model"]["select_snow"] = ["common_snow_no_storage"]
-config["model"]["select_interception"] =
-    ["interception_no_storage", "interception_overflow", "interception_overflow"]
-config["model"]["select_hortonponding"] =
-    ["hortonponding_no_storage", "hortonponding", "hortonponding_no_storage"]
-config["model"]["select_hortonrunoff"] =
-    ["hortonrunoff_no_storage", "hortonrunoff", "hortonrunoff_no_storage"]
-config["model"]["select_rootzone"] =
-    ["rootzone_storage", "rootzone_storage", "rootzone_no_storage"]
-config["model"]["select_fast"] = ["fast_no_storage", "fast_storage", "fast_storage"]
-config["model"]["select_slow"] = ["common_slow_storage"]
-
-model = Wflow.initialize_flextopo_model(config)
-@unpack network = model
-
-model = Wflow.run_timestep(model)
-
-@testset "first timestep" begin
-    flextopo = model.vertical
-    @test flextopo.rootzonestorage[3500] ≈
-          [147.11238663084805f0, 79.08369375691255f0, 0.0f0]
-    @test flextopo.runoff[3500] ≈ 0.19008129369467497f0
-    @test flextopo.rootevap[3500] ≈ [0.009225917980074883f0, 0.009225917980074883f0, 0.0f0]
-    @test flextopo.snow[3500] == 0.0f0
-end
-
-@testset "overland domain" begin
-    q = model.lateral.land.q_av
-    land = model.lateral.land
-    @test sum(q) ≈ 69.51434665542678f0
-    @test q[10354] ≈ 0.0009376627801121855f0
-    @test land.volume[10354] ≈ 63.46997463442048f0
-    @test land.inwater[10354] ≈ 0.0016722689680105701f0
-    @test q[network.land.order[end]] ≈ 0.0006175366185364747f0
-end
-
-@testset "river flow" begin
-    q = model.lateral.river.q_av
-    @test sum(q) ≈ 22.236556345676156f0
-    @test q[651] ≈ 0.00035906990499781237f0
-    @test q[1056] ≈ 0.0002647264689318855f0
-    @test q[network.river.order[end]] ≈ 0.000015647838277869928f0
-end
-
-Wflow.close_files(model, delete_output = false)
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
deleted file mode 100644
index 12d6eb169..000000000
--- a/test/run_hbv.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-
-tomlpath = joinpath(@__DIR__, "hbv_config.toml")
-config = Wflow.Config(tomlpath)
-
-model = Wflow.initialize_hbv_model(config)
-@unpack network = model
-
-model = Wflow.run_timestep(model)
-
-# test if the first timestep was written to the CSV file
-flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
-@testset "CSV output" begin
-    row = csv_first_row(model.writer.csv_path)
-
-    @test row.time == DateTime("2000-01-01T00:00:00")
-    @test row.Q ≈ 521.8433822888003f0
-    @test row.temp_bycoord ≈ 2.965437173843384f0
-    @test row.temp_byindex ≈ 1.1716821193695068f0
-    @test row.Q_1 ≈ 505.1935875677504f0
-    @test row.perc_33 ≈ 2.308000087738037f0
-    @test row.perc_34 ≈ 1.8980000019073486f0
-    @test row.perc_35 ≈ 2.7100000381469727f0
-    @test row.perc_36 ≈ 3.818000078201294f0
-    @test row.perc_37 ≈ 2.1440000534057617f0
-end
-
-@testset "first timestep" begin
-    hbv = model.vertical
-    @test hbv.tt[4377] ≈ 0.0
-    @test model.clock.iteration == 1
-    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
-    @test hbv.runoff[4377] ≈ 7.406898120121746f0
-    @test hbv.soilevap[4377] == 0.0
-    @test hbv.snow[4377] == 0.0
-end
-
-# run the second timestep
-model = Wflow.run_timestep(model)
-
-@testset "second timestep" begin
-    hbv = model.vertical
-    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
-    @test hbv.runoff[4377] ≈ 4.3533463f0
-    @test hbv.soilevap[4377] == 0.0
-    @test hbv.snow[4377] == 0.0
-end
-
-@testset "overland domain" begin
-    q = model.lateral.land.q_av
-    land = model.lateral.land
-    @test sum(q) ≈ 3264.157286198723f0
-    @test q[10354] ≈ 0.2350958683414288f0
-    @test land.volume[10354] ≈ 2057.7314802432425f0
-    @test land.inwater[10354] ≈ 0.027351853491789667f0
-    @test q[network.land.order[end]] ≈ 0.27158713754634217f0
-end
-
-@testset "river flow" begin
-    q = model.lateral.river.q_av
-    @test sum(q) ≈ 52686.97949266187f0
-    @test q[651] ≈ 5.7624898709967125f0
-    @test q[1056] ≈ 8.892291220438524f0
-    @test q[network.river.order[end]] ≈ 360.40425076323913f0
-end
-
-Wflow.close_files(model, delete_output = false)
diff --git a/test/runtests.jl b/test/runtests.jl
index 18771d1de..440c640a3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -38,11 +38,7 @@ end
 staticmaps_rhine_path = testdata(v"0.1", "staticmaps.nc", "staticmaps-rhine.nc")
 staticmaps_moselle_path =
     testdata(v"0.2.9", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
-staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn.nc")
-staticmaps_meuse_path =
-    testdata(v"0.2.8", "staticmaps_flex_meuse.nc", "staticmaps_flex_meuse.nc")
 forcing_moselle_path = testdata(v"0.2.6", "forcing-moselle.nc", "forcing-moselle.nc")
-forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
 forcing_moselle_sed_path =
     testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
 staticmaps_moselle_sed_path =
@@ -60,7 +56,6 @@ forcing_sbm_gw_path = testdata(
     "forcing-sbm-groundwater-part2.nc",
     "forcing-sbm-groundwater-part2.nc",
 )
-forcing_meuse_path = testdata(v"0.2.8", "forcing_meuse.nc", "forcing_meuse.nc")
 staticmaps_sbm_gw_path =
     testdata(v"0.2.3", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
 instates_sbm_gw_path =
@@ -85,7 +80,6 @@ with_logger(NullLogger()) do
         include("vertical_process.jl")
         include("reservoir_lake.jl")
         include("run_sbm.jl")
-        include("run_hbv.jl")
         include("run_sbm_gwf.jl")
         include("run.jl")
         include("groundwater.jl")
@@ -93,7 +87,6 @@ with_logger(NullLogger()) do
         include("bmi.jl")
         include("run_sediment.jl")
         include("subdomains.jl")
-        include("run_flextopo.jl")
 
         Aqua.test_all(Wflow; ambiguities = false)
     end