diff --git a/docs/developments/julia_structs.qmd b/docs/developments/julia_structs.qmd
index 0ed155b9d..c4827cffb 100644
--- a/docs/developments/julia_structs.qmd
+++ b/docs/developments/julia_structs.qmd
@@ -35,11 +35,11 @@ The `subsurface` part is mapped to the lateral subsurface flow kinematic wave co
 `land` part is mapped the overland flow kinematic wave concept and the `river` part is mapped
 to the river flow kinematic wave concept. Knowledge of this specific mapping is required to
 understand and correctly set input, output and state variables in the TOML configuration file,
-see also [Config and TOML](@ref config_toml). This mapping is described in more detail for each
+see also [Config and TOML](../user_guide/toml_file.qmd). This mapping is described in more detail for each
 model in the section Models. Also the `struct` of each mapped concept is provided, so one can
 check the internal variables in the code. These structs are defined as a parametric composite
 type, with type parameters between curly braces after the `struct` name. See also the next
-paragraph [Vertical and lateral models](@ref) for a more detailed description.
+paragraph [Vertical and lateral models](#vertical-and-lateral-models) for a more detailed description.
 
 ## Vertical and lateral models
 The different model concepts used in wflow are defined as parametric [composite
diff --git a/docs/getting_started/building_a_model.qmd b/docs/getting_started/building_a_model.qmd
index 5f7b1c8b0..17c8afa83 100644
--- a/docs/getting_started/building_a_model.qmd
+++ b/docs/getting_started/building_a_model.qmd
@@ -11,9 +11,9 @@ model schematization, (dynamic) forcing data, results and states.
 The wflow plugin [hydroMT-wflow](https://github.com/Deltares/hydromt_wflow) of hydroMT can be
 used to build and analyze the following model configurations:
 
-- [wflow\_sbm + kinematic wave routing](../model_docs/model_configurations.html#sbm-kinematic-wave)
-- [wflow\_sbm + local inertial river and floodplain](../model_docs/model_configurations.html#sbm-local-inertial-river)
-- [wflow\_sbm + local inertial river (1D) and land (2D)](../model_docs/model_configurations.html#sbm-local-inertial-river-1d-and-land-2d)
+- [wflow\_sbm + kinematic wave routing](../model_docs/model_configurations.qmd#sbm-kinematic-wave)
+- [wflow\_sbm + local inertial river and floodplain](../model_docs/model_configurations.qmd#sbm-local-inertial-river)
+- [wflow\_sbm + local inertial river (1D) and land (2D)](../model_docs/model_configurations.qmd#sbm-local-inertial-river-1d-and-land-2d)
 - [wflow\_sediment](../model_docs/model_configurations.html#wflow_sediment)
 
 To learn more about the wflow plugin of this Python package, we refer to the [hydroMT-wflow
@@ -21,8 +21,8 @@ documentation](https://deltares.github.io/hydromt_wflow/latest/index.html).
 
 To inspect or modify (for example in QGIS) the netCDF static data of these wflow models it is
 convenient to export the maps to a raster format. This can be done as part of the hydroMT-wflow
-plugin, see also the following [example]
-(https://deltares.github.io/hydromt_wflow/latest/_examples/convert_staticmaps_to_mapstack.html).
+plugin, see also the following
+[example](https://deltares.github.io/hydromt_wflow/latest/_examples/convert_staticmaps_to_mapstack.html).
 It is also possible to create again the netCDF static data file based on the modified raster
 map stack.
 
@@ -31,7 +31,8 @@ map stack.
 The actual data requirements depend on the application of the Model and the Model type. Both
 forcing and static data should be provided in netCDF format, with the same grid definition for
 forcing and static data. The only exception is storage and rating curves for lakes, that should
-be provided in CSV format, see also [Additional settings](@ref).
+be provided in CSV format, see also [Additional settings for
+waterbodies](../model_docs/lateral/waterbodies.qmd#additional-settings).
 
 * Forcing data:
   - Precipitation
@@ -57,65 +58,80 @@ An approach to generate `ldd` data is to make use of the Python package
   data](https://deltares.github.io/pyflwdir/latest/_examples/from_dem.html),
 
 see also Eilander et al. (2021) for more information. Pyflwdir is also used by the
-[hydroMT](@ref) Python package described in the next paragraph. Another approach to generate
-`ldd` data is to make use of PCRaster functionality, see for example
+[hydroMT](https://github.com/Deltares/hydromt) Python package described in the next paragraph.
+Another approach to generate `ldd` data is to make use of PCRaster functionality, see for
+example
 [lddcreate](https://pcraster.geo.uu.nl/pcraster/4.3.1/documentation/pcraster_manual/sphinx/op_lddcreate.html).
 
 Optionally, but also not directly part of a model component are `gauge` locations, that are
 used to extract gridded data from certain locations.
 
 The different supported model configurations are described in the section [Model
-configurations](@ref). Wflow\_sbm models have the vertical concept [SBM](@ref vert_sbm) in
-common and input parameters for this component are described in the [SBM](@ref params_sbm)
-section of Model parameters. For wflow\_sbm models there are two ways to include subsurface
-flow:
-
-1. The kinematic wave approach (see section [Subsurface flow routing](@ref)) as part of the
-   `sbm` model type. Parameters that are part of this component are described in the [Lateral
-   subsurface flow](@ref params_ssf) section of Model parameters. Input parameters for this
-   component are derived from the SBM vertical concept and the land slope. One external
-   parameter [`ksathorfrac`](@ref params_ssf) is used to calculate the horizontal hydraulic
-   conductivity at the soil surface `kh_0`.
-2. Groundwater flow (see section [Groundwater flow component](@ref lateral_gwf)) as part of the
-   `sbm_gwf` model type. For the unconfined aquifer the input parameters are described in the
-   section [Unconfined aquifer](@ref) of Model parameters. The bottom (`bottom`) of the
-   groundwater layer is derived from from the `soilthickness` [mm] parameter of `SBM` and the
-   provided surface elevation `altitude` [m] as part of the static input. The `area` parameter
-   is derived from the model grid. Parameters that are part of the boundary conditions of the
-   unconfined aquifer are listed under [Constant Head](@ref) and [Boundary conditions](@ref) of
-   the Model parameters section.
+configurations](../model_docs/model_configurations.qmd). Wflow\_sbm models have the vertical
+concept [SBM](../model_docs/vertical/sbm.qmd) in common and input parameters for this component
+are described in the [SBM model parameters table](../model_docs/parameters_vertical.qmd#sbm).
+For wflow\_sbm models there are two ways to include subsurface flow:
+
+1. The kinematic wave approach (see section [Subsurface flow
+   routing](../model_docs/lateral/kinwave.qmd#subsurface-flow-routing)) as part of the `sbm`
+   model type. Parameters that are part of this component are described in the [Lateral
+   subsurface flow](../model_docs/parameters_lateral.qmd#lateral-subsurface-flow) section of
+   Model parameters. Input parameters for this component are derived from the SBM vertical
+   concept and the land slope. One external parameter (`ksathorfrac`) is used to calculate the
+   horizontal hydraulic conductivity at the soil surface `kh_0`.
+2. Groundwater flow (see section [Groundwater flow component](../model_docs/lateral/gwf.qmd))
+   as part of the `sbm_gwf` model type. For the unconfined aquifer the input parameters are
+   described in the section [Unconfined
+   aquifer](../model_docs/parameters_lateral.qmd#unconfined-aquifer) of Model parameters. The
+   bottom (`bottom`) of the groundwater layer is derived from from the `soilthickness` [mm]
+   parameter of `SBM` and the provided surface elevation `altitude` [m] as part of the static
+   input. The `area` parameter is derived from the model grid. Parameters that are part of the
+   boundary conditions of the unconfined aquifer are listed under [Constant
+   Head](../model_docs/lateral/gwf.qmd#head-boundary) and [Boundary
+   conditions](../model_docs/parameters_lateral.qmd#boundary-conditions) of the Model
+   parameters section.
 
 Most hydrological model configurations make use of the kinematic wave surface routing (river
 flow, overland flow or both) and input data required for the river and overland flow components
-is described in [Surface flow](@ref).  There is also the option to use the local inertial model
-as part of the wflow\_sbm models (model types `sbm` and `sbm_gwf`):
-+ for river flow, see also the [Local inertial river and floodplain](@ref
-  config_sbm_gwf_lie_river) model.
-+ for 1D river flow and 2D overland flow combined, see also the [Local inertial river (1D)
-  and land (2D)](@ref config_sbm_gwf_lie_river_land) model.
-
-Input parameters for this approach are described in [River flow (local inertial)](@ref
-local-inertial_river_params), including the optional 1D [floodplain schematization](@ref
-local-inertial_floodplain_params), and [Overland flow (local inertial)](@ref
-local-inertial_land_params) of the Model parameters section.
+is described in [Surface flow](../model_docs/parameters_lateral.qmd#surface-flow). There is
+also the option to use the local inertial model as part of the wflow\_sbm models (model types
+`sbm` and `sbm_gwf`):
+
++ for river flow, see also the [Local inertial river and
+  floodplain](../model_docs/model_configurations.qmd#sbm-local-inertial-river) model.
++ for 1D river flow and 2D overland flow combined, see also the [Local inertial river (1D) and
+  land (2D)](../model_docs/model_configurations.qmd#sbm-local-inertial-river-1d-and-land-2d)
+  model.
+
+Input parameters for this approach are described in [River flow (local
+inertial)](../model_docs/parameters_lateral.qmd#local-inertial), including the optional 1D
+[floodplain schematization](../model_docs/parameters_lateral.qmd#d-floodplain), and [Overland
+flow (local inertial)](../model_docs/parameters_lateral.qmd#overland-flow) of the Model
+parameters section.
 
 Reservoirs or lakes can be part of the kinematic wave or local inertial model for river flow
-and input parameters are described in [Reservoirs](@ref reservoir_params) and [Lakes](@ref
-lake_params).
-
-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 lateral
-[Sediment Flux in overland flow](@ref) concept are described in the [Overland flow](@ref)
-section of the Model parameters. Parameters of this component are not directly set by data from
-static input. The input parameters of the lateral concept [River Sediment Model](@ref) are
-listed in [River flow](@ref) of the Model parameters section.
+and input parameters are described in
+[Reservoirs](../model_docs/parameters_lateral.qmd#reservoirs) and
+[Lakes](../model_docs/parameters_lateral.qmd#lakes).
+
+The [wflow\_sediment](../model_docs/model_configurations.qmd#wflow_sediment) model
+configuration consists of the vertical [Soil Erosion](../model_docs/vertical/sediment.qmd)
+concept and the input parameters for this concept are described in the
+[Sediment](../model_docs/parameters_vertical.qmd#sediment) section of the Model parameters. The
+parameters of the lateral [Sediment Flux in overland
+flow](../model_docs/lateral/sediment_flux.qmd#sediment-flux-in-overland-flow) concept are
+described in the [Overland flow](../model_docs/parameters_lateral.qmd#overland-flow-1) section
+of the Model parameters. Parameters of this component are not directly set by data from static
+input. The input parameters of the lateral concept [River Sediment
+Model](../model_docs/lateral/sediment_flux.qmd#river-sediment-model) are listed in [River
+flow](../model_docs/parameters_lateral.qmd#river-flow-1) of the Model parameters section.
 
 The Model parameters section lists all the parameters per Model component and these Tables can
 also be used to check which parameters can be part of the output, see also [Output netCDF
-section](@ref) and [Output CSV section](@ref).
+section](../user_guide/toml_file.qmd#output-netcdf-section) and [Output CSV
+section](../user_guide/toml_file.qmd#output-csv-section).
 
-Example models can be found in the [Example models section](@ref sample_data).
+Example models can be found in the [Example models section](./download_example_models.qmd).
 
 ## References
 + Yamazaki, D., Ikeshima, D., Sosa, J., Bates, P. D., Allen, G. H. and Pavelsky, T. M.: MERIT
diff --git a/docs/getting_started/download_example_models.qmd b/docs/getting_started/download_example_models.qmd
index 0c48853a5..67917dd1c 100644
--- a/docs/getting_started/download_example_models.qmd
+++ b/docs/getting_started/download_example_models.qmd
@@ -14,8 +14,8 @@ listed in the Table below:
 
 The associated Model files (input static, forcing and state files) can easily be downloaded and
 for this we share the following Julia code (per Model) that downloads the required files to
-your current working directory. For running these test model see also [Usage](@ref run_wflow)
-and [Command Line Interface](@ref cli).
+your current working directory. For running these test model see also [Usage](./running_wflow.qmd#running-a-simulation)
+and [Command Line Interface](./running_wflow.qmd#using-the-command-line-interface).
 
 ## wflow\_sbm + kinematic wave
 ```julia
diff --git a/docs/getting_started/running_wflow.qmd b/docs/getting_started/running_wflow.qmd
index b59194a4f..50f47fb1d 100644
--- a/docs/getting_started/running_wflow.qmd
+++ b/docs/getting_started/running_wflow.qmd
@@ -12,7 +12,7 @@ julia -e 'using Wflow; Wflow.run()' path/to/config.toml
 ```
 
 Furthermore, it is possible to write a Julia script to run a simulation. Example data to
-explore how this works can be found [here](@ref sample_data).
+explore how this works can be found [here](./download_example_models.qmd).
 
 ```julia
 using Wflow
diff --git a/docs/model_docs/index.qmd b/docs/model_docs/index.qmd
index 54ddc17ab..969f89cb1 100644
--- a/docs/model_docs/index.qmd
+++ b/docs/model_docs/index.qmd
@@ -8,8 +8,8 @@ modelling framework of wflow. Descriptions are given regarding the model concept
 the original scientific papers which explain the concepts in more detail. The model parameters
 which influence the processes are also shown, using inline code blocks. An overview of all
 model parameters is also provided for easy reference, including their short names, long
-descriptions and their units (see [parameters vertical concepts](@ref params_vert) and
-[parameters lateral concepts](@ref params_lat)).
+descriptions and their units (see [parameters vertical concepts](./parameters_vertical.qmd) and
+[parameters lateral concepts](./parameters_lateral.qmd)).
 
 ## Division between vertical and lateral
 
diff --git a/docs/model_docs/lateral/gwf.qmd b/docs/model_docs/lateral/gwf.qmd
index 7e17e8f41..d24c8a5bd 100644
--- a/docs/model_docs/lateral/gwf.qmd
+++ b/docs/model_docs/lateral/gwf.qmd
@@ -38,7 +38,7 @@ where $\SIb{S}{m m^{-1}}$ is storativity (or specific yield), $\SIb{\phi}{m}$ is
 head, $t$ is time, $\SIb{k}{m t^{-1}}$ is horizontal hydraulic conductivity, $\SIb{H}{m}$ is
 the (saturated) aquifer height: groundwater level - aquifer bottom elevation and $\SIb{Q}{m
 t^{-1}}$ represents fluxes from boundary conditions (e.g. recharge or abstraction), see also
-[Aquifer boundary conditions](@ref).
+[Aquifer boundary conditions](#aquifer-boundary-conditions).
 
 The simplest finite difference formulation is forward in time, central in space, and can be
 written as:
@@ -50,7 +50,7 @@ $$
 where $i$ is the cell index, $t$ is time, $\Delta t$ is the step size, $C_{i-1}$ is the the
 intercell conductance between cell $i-1$ and $i$ and $C_i$ is the intercell conductance between
 cell $i$ and $i+1$. The connection data between cells is stored as part of the `Connectivity`
-struct, see also [Connectivity](@ref) for more information.
+struct, see also [Connectivity](#connectivity) for more information.
 
 Conductance $C$ is defined as:
 
@@ -76,7 +76,7 @@ saturated thickness of a cell-to-cell is approximated using the cell with the hi
 This results in a consistent overestimation of the saturated thickness, but it avoids
 complexities related with cell drying and rewetting, such as having to define a "wetting
 threshold" or a "wetting factor". See also the documentation for MODFLOW-NWT (Niswonger et al.,
-2011) or MODFLOW6 (Langevin et al., 2017) for more background information. For more background
+1)    or MODFLOW6 (Langevin et al., 2017) for more background information. For more background
 on drying and rewetting, see for example McDonald et al. (1991).
 
 For the finite difference formulation, there is only one unknown, $\phi_i^{t+1}$. Reshuffling
@@ -164,9 +164,9 @@ the river bed exfiltration conductance, $\SIb{\subtext{B}{riv}}{L}$ the bottom o
 bed, $\SIb{\subtext{h}{riv}}{L}$ is the river stage and $\SIb{\phi}{L}$ is the hydraulic head
 in the river cell.
 
-The Table in the Groundwater flow [river boundary condition](@ref gwf_river_params) section of
-the Model parameters provides the parameters of the struct `River`. Parameters that can be set
-directly from the static input data (netCDF) are marked in this Table.
+The Table in the Groundwater flow [river boundary condition](../parameters_lateral.qmd#river)
+section of the Model parameters provides the parameters of the struct `River`. Parameters that
+can be set directly from the static input data (netCDF) are marked in this Table.
 
 The exchange flux (river to aquifer) $\subtext{Q}{riv}$ is an output variable (field `flux` of
 the `River` struct), and is used to update the total flux in a river cell. For the model `SBM +
@@ -184,9 +184,10 @@ where $\SIb{\subtext{Q}{drain}}{L^3 T^{-1}}$ is the exchange flux from drains to
 $\SIb{\subtext{C}{drain}}{L^2 T^{-1}}$ is the drain conductance, $\SIb{\subtext{h}{drain}}{L}$
 is the drain elevation and $\SIb{\phi}{L}$ is the hydraulic head in the cell with drainage.
 
-The table in the Groundwater flow [drainage boundary condition](@ref gwf_drainage_params)
-section of the Model parameters provides the parameters of the struct `Drainage`. Parameters
-that can be set directly from the static input data (netCDF) are marked in this Table.
+The table in the Groundwater flow [drainage boundary
+condition](../parameters_lateral.qmd#drainage) section of the Model parameters provides the
+parameters of the struct `Drainage`. Parameters that can be set directly from the static input
+data (netCDF) are marked in this Table.
 
 The exchange flux (drains to aquifer) $\subtext{Q}{drain}$ is an output variable (field `flux`
 of struct `Drainage`), and is used to update the total flux in a cell with drains. For the
@@ -208,10 +209,10 @@ $$
 
 with $\SIb{}{L T^{-1}}$ the recharge rate and $\SIb{A}{L^2}$ the area of the aquifer cell.
 
-The table in the Groundwater flow [recharge boundary condition](@ref gwf_recharge_params)
-section of the Model parameters section provides the parameters of the struct `Recharge`.
-Parameters that can be set directly from the static input data (netCDF) are marked in this
-Table.
+The table in the Groundwater flow [recharge boundary
+condition](../parameters_lateral.qmd#recharge) section of the Model parameters section provides
+the parameters of the struct `Recharge`. Parameters that can be set directly from the static
+input data (netCDF) are marked in this Table.
 
 The recharge flux $Q_r$ is an output variable (field `flux` of struct `Recharge`), and is used
 to update the total flux in a cell where recharge occurs. For the model `SBM + Groundwater
@@ -232,8 +233,9 @@ $$
 with $\SIb{C_{hb}}{L^2 T^{-1}}$ the conductance of the head boundary, $\SIb{\phi_{hb}}{L}$ the
 head of the head boundary and $\phi$ the head of the aquifer cell.
 
-The table in the Groundwater flow [head boundary condition](@ref gwf_headboundary_params)
-section of the Model parameters provides the parameters of the struct `HeadBoundary`.
+The table in the Groundwater flow [head boundary
+condition](../parameters_lateral.qmd#head-boundary) section of the Model parameters provides
+the parameters of the struct `HeadBoundary`.
 
 The head boundary flux $Q_{hb}$ is an output variable (field `flux` of struct `HeadBoundary`),
 and is used to update the total flux in a cell where this type of boundary occurs. The
@@ -246,8 +248,8 @@ This boundary is not (yet) part of the model `SBM + Groundwater flow`.
 ### Well boundary
 A volumetric well rate $\SIb{}{L^3 T^{-1}}$ can be specified as a boundary condition.
 
-The Table in the [well boundary condition](@ref well_boundary_params) section of the Model
-parameters provides the parameters of the struct `Well`.
+The Table in the [well boundary condition](../parameters_lateral.qmd#well-boundary) section of
+the Model parameters provides the parameters of the struct `Well`.
 
 The volumetric well rate $\subtext{Q}{well}$ can be can be specified as a fixed or time
 dependent value. If a cell is dry, the actual well flux `flux` is set to zero (see also the
diff --git a/docs/model_docs/lateral/kinwave.qmd b/docs/model_docs/lateral/kinwave.qmd
index 0c3e4ccb5..b026d158a 100644
--- a/docs/model_docs/lateral/kinwave.qmd
+++ b/docs/model_docs/lateral/kinwave.qmd
@@ -53,13 +53,13 @@ bankfull_depth = "wflow_riverdepth"
 The wetted perimeter of the river is based on half bankfull river depth. For the land part the
 wetted perimeter is based on the flow width.
 
-Simplified [reservoir and lake](@ref reservoir_lake) models can be included as part of the
-river kinematic wave network.
+Simplified [reservoir and lake](./waterbodies.qmd) models can be included as part of the river
+kinematic wave network.
 
 ## Inflow
 External water (supply/abstraction) `inflow` $\SIb{}{m^3 s^{-1}}$  can be added to the
 kinematic wave for surface water routing, as a cyclic parameter or as part of forcing (see also
-[Input section](@ref)).
+[Input section](../../user_guide/required_files.qmd)).
 
 ## Abstractions
 Abstractions from the river through the variable `abstraction` $\SIb{}{m^3 s^{-1}}$ are
@@ -71,9 +71,9 @@ routing scheme for river flow.
 ## Subsurface flow routing
 In the SBM model the kinematic wave approach is used to route subsurface flow laterally.
 Different vertical hydraulic conductivity depth profiles are possible as part of the vertical
-[SBM](@ref soil) concept, and these profiles (after unit conversion) are also used to compute
-lateral subsurface flow. The following profiles (see [SBM](@ref soil) for a detailed
-description) are available:
+[SBM](../vertical/sbm.qmd#soil-processes) concept, and these profiles (after unit conversion)
+are also used to compute lateral subsurface flow. The following profiles (see
+[SBM](../vertical/sbm.qmd#soil-processes) for a detailed description) are available:
 - `exponential` (default)
 - `exponential_constant`
 - `layered`
@@ -190,10 +190,11 @@ scheme. In summary the following limitations apply:
 ## External inflows
 External inflows, for example water supply or abstractions, can be added to the kinematic wave
 via the `inflow` variable. For this, the user can supply a 2D map of the inflow, as a cyclic
-parameter or as part of forcing (see also [Input section](@ref)). These inflows are added or
-abstracted from the upstream inflow `qin` before running the kinematic wave to solve the impact
-on resulting `q`. In case of a negative inflow (abstractions), a minimum of zero is applied to
-the upstream flow `qin`.
+parameter or as part of forcing (see also [Input
+section](../../user_guide/required_files.qmd)). These inflows are added or abstracted from the
+upstream inflow `qin` before running the kinematic wave to solve the impact on resulting `q`.
+In case of a negative inflow (abstractions), a minimum of zero is applied to the upstream flow
+`qin`.
 
 ## References
 + Chow, V., Maidment, D. and Mays, L., 1988, Applied Hydrology. McGraw-Hill Book Company,
diff --git a/docs/model_docs/lateral/local-inertial.qmd b/docs/model_docs/lateral/local-inertial.qmd
index 8f9bc6466..00987d27e 100644
--- a/docs/model_docs/lateral/local-inertial.qmd
+++ b/docs/model_docs/lateral/local-inertial.qmd
@@ -93,13 +93,13 @@ river cell. The river length $\SIb{}{m}$ of the boundary cell can be set through
 with `riverlength_bc`, and has a default value of $\SI{10}{km}$. The water depth at the
 boundary cell is fixed at $\SI{0.0}{m}$.
 
-Simplified [reservoir and lake](@ref reservoir_lake) models can be included as part of the
-local inertial model for river flow (1D) and river and overland flow combined (see next
-section). Reservoir and lake models are included as a boundary point with zero water depth for
-both river and overland flow. For river flow the reservoir or lake model replaces the local
-inertial model at the reservoir or lake location, and $Q$ is set by the outflow from the
-reservoir or lake. Overland flow at a reservoir or lake location is not allowed to or from the
-downstream river grid cell.
+Simplified [reservoir and lake](./waterbodies.qmd) models can be included as part of the local
+inertial model for river flow (1D) and river and overland flow combined (see next section).
+Reservoir and lake models are included as a boundary point with zero water depth for both river
+and overland flow. For river flow the reservoir or lake model replaces the local inertial model
+at the reservoir or lake location, and $Q$ is set by the outflow from the reservoir or lake.
+Overland flow at a reservoir or lake location is not allowed to or from the downstream river
+grid cell.
 
 ## Overland flow (2D)
 For the simulation of 2D overland flow on a staggered grid the numerical scheme proposed by de
@@ -148,12 +148,13 @@ h_thresh = 0.1                   # water depth [m] threshold for calculating flo
 The properties `inertial_flow_alpha`, `froude_limit` and `h_thresh` apply to 1D river routing
 as well as 2D overland flow. The properties `inertial_flow_alpha` and `froude_limit`, and the
 adaptive model time step $\Delta t$ are explained in more detail in the [River and floodplain
-routing](@ref) section of the local inertial model.
+routing](#river-and-floodplain-routing) section of the local inertial model.
 
 ## Inflow
 External water (supply/abstraction) `inflow` $\SIb{}{m^3 s^{-1}}$ can be added to the local
 inertial model for river flow (1D) and river and overland flow combined (1D-2D), as a cyclic
-parameter or as part of forcing (see also [Input section](@ref)).
+parameter or as part of forcing (see also [Input
+section](../../user_guide/required_files.qmd)).
 
 ## Abstractions
 Abstractions from the river through the variable `abstraction` $\SIb{}{m^3 s^{-1}}$ are
diff --git a/docs/model_docs/model_configurations.qmd b/docs/model_docs/model_configurations.qmd
index 148194413..5fa436500 100644
--- a/docs/model_docs/model_configurations.qmd
+++ b/docs/model_docs/model_configurations.qmd
@@ -4,7 +4,7 @@ title: Model configurations
 
 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` [SBM](@ref vert_sbm) concept in combination
+wflow model configuration consists of a `vertical` [SBM](./vertical/sbm.qmd) 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:
@@ -18,17 +18,18 @@ model configurations are supported in wflow:
  - wflow\_sediment as post processing of wflow\_sbm output
 
 Below, some explanation 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 data](@ref sample_data).
+for other model configurations is provided in the section with [sample
+data](../getting_started/download_example_models.qmd).
 
 ## wflow\_sbm
 
 Wflow\_sbm represents hydrological models derived from the CQflow model (Köhler et al., 2006)
-that have the [SBM](@ref vert_sbm) vertical concept in common, but can have different lateral
-concepts that control how water is routed for example over the land or river domain. The soil
-part of SBM is largely based on the Topog\_SBM model but has had considerable changes over
-time. Topog\_SBM is specifically designed to simulate fast runoff processes in small catchments
-while the wflow\_sbm model can be applied more widely. The main differences are for the
-vertical concept SBM of wflow\_sbm:
+that have the [SBM](./vertical/sbm.qmd) vertical concept in common, but can have different
+lateral concepts that control how water is routed for example over the land or river domain.
+The soil part of SBM is largely based on the Topog\_SBM model but has had considerable changes
+over time. Topog\_SBM is specifically designed to simulate fast runoff processes in small
+catchments while the wflow\_sbm model can be applied more widely. The main differences are for
+the vertical concept SBM of wflow\_sbm:
 
 - The unsaturated zone can be split-up in different layers
 - The addition of evapotranspiration losses
@@ -38,13 +39,13 @@ vertical concept SBM of wflow\_sbm:
 The water demand and allocation computations are supported by the wflow\_sbm model
 configurations:
 
-- [SBM + Kinematic wave](@ref config_sbm)
-- [SBM + Groundwater flow](@ref config_sbm_gwf)
-- [SBM + Local inertial river](@ref config_sbm_gwf_lie_river)
-- [SBM + Local inertial river (1D) and land (2D)](@ref config_sbm_gwf_lie_river_land)
+- SBM + Kinematic wave
+- SBM + Groundwater flow
+- SBM + Local inertial river
+- SBM + Local inertial river (1D) and land (2D)
 
 The vertical SBM concept is explained in more detail in the following section [SBM vertical
-concept](@ref vert_sbm).
+concept](./vertical/sbm.qmd).
 
 Topog\_SBM uses an element network based on contour lines and trajectories for water routing.
 Wflow\_sbm models differ in how the lateral components river, land, and subsurface are solved.
@@ -53,7 +54,7 @@ Below the different wflow\_sbm model configurations are described.
 ### SBM + Kinematic wave
 For the lateral components of this wflow\_sbm model water is routed over a D8 network, and the
 kinematic wave approach is used for river, overland and lateral subsurface flow. This is
-described in more detail in the section [Kinematic wave](@ref kin_wave).
+described in more detail in the section [Kinematic wave](./lateral/kinwave.qmd).
 
 An overview of the different processes and fluxes in the wflow_sbm model with the kinematic
 wave approach for river, overland and lateral subsurface flow:
@@ -78,7 +79,7 @@ lateral.river.reservoir => struct SimpleReservoir{T} # optional
 For river and overland flow the kinematic wave approach over a D8 network is used for this
 wflow\_sbm model. For the subsurface domain, an unconfined aquifer with groundwater flow in
 four directions (adjacent cells) is used. This is described in more detail in the section
-[Groundwater flow](@ref lateral_gwf).
+[Groundwater flow](./lateral/gwf.qmd).
 
 ```toml
 [model]
@@ -151,7 +152,7 @@ lateral.land => struct ShallowWaterLand{T}
 ```
 
 The local inertial approach is described in more detail in the section [Local inertial
-model](@ref local_inertial).
+model](./lateral/local-inertial.qmd).
 
 ## wflow\_sediment
 The processes and fate of many particles and pollutants impacting water quality at the
@@ -171,8 +172,10 @@ In order to model the exports of terrestrial sediment to the coast through the L
 Aquatic Continuum or LOAC (inland waters network such as streams, lakes...), two different
 modelling parts were considered. The first part, called the inland sediment model, is the
 modelling and estimation of soil loss and sediment yield to the river system by land erosion,
-separated into vertical [Soil Erosion](@ref) processes and lateral [Sediment Flux in overland
-flow](@ref). The second part, called the [River Sediment Model](@ref) is the transport and
+separated into vertical [Soil Erosion](./vertical/sediment.qmd#soil-erosion) processes and
+lateral [Sediment Flux in overland
+flow](./lateral/sediment_flux.qmd#sediment-flux-in-overland-flow). The second part, called the
+[River Sediment Model](./lateral/sediment_flux.qmd#river-sediment-model) is the transport and
 processes of the sediment in the river system. The two parts together constitute the
 wflow\_sediment model.
 
diff --git a/docs/model_docs/parameters_lateral.qmd b/docs/model_docs/parameters_lateral.qmd
index f97607965..e5913806c 100644
--- a/docs/model_docs/parameters_lateral.qmd
+++ b/docs/model_docs/parameters_lateral.qmd
@@ -108,10 +108,10 @@ The `khfrac` parameter compensates for anisotropy, small scale `kv_0` measuremen
 that do not represent larger scale hydraulic conductivity, and smaller flow length scales
 (hillslope) in reality, not represented by the model resolution.
 
-For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity depth
-profiles are possible, and these also determine which `LateralSSF` parameters are used
-including the input requirements for the computation of lateral subsurface flow. For the
-`exponential` profile the model parameters `kh_0` and `f` are used. For the
+For the vertical [SBM](./parameters_vertical.qmd#sbm) concept different vertical hydraulic
+conductivity depth profiles are possible, and these also determine which `LateralSSF`
+parameters are used including the input requirements for the computation of lateral subsurface
+flow. For the `exponential` profile the model parameters `kh_0` and `f` are used. For the
 `exponential_constant` profile `kh_0` and `f` are used, and `z_exp` is required as part of
 `[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for the
 `layered_exponential` profile `kv` is used and `z_exp` is required as part of
@@ -311,6 +311,88 @@ internal model parameter `z`, and is listed in the Table below between parenthes
 | `rivercells`  |  river cells| - | - |
 | `h_av` | average water depth| m | - |
 
+## Water bodies
+
+### Reservoirs
+The Table below shows the parameters (fields) of struct `SimpleReservoir`, 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 input data (netCDF), and can be
+listed in the TOML configuration file under `[input.lateral.river.reservoir]`, to map the
+internal model parameter to the external netCDF variable.
+
+Two parameters reservoir coverage `areas` and the outlet of reservoirs (unique id) `locs` that
+are not part of the `SimpleReservoir` struct are also required, and can be set as follows
+through the TOML file:
+
+```toml
+[input.lateral.river.reservoir]
+areas = "wflow_reservoirareas"
+locs = "wflow_reservoirlocs"
+```
+
+|  parameter | description  	        | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| **`area`**        | area | m``^2`` | - |
+| **`demand`**  | minimum (environmental) flow requirement downstream of the reservoir | m``^3`` s``^{-1}``| - |
+| **`maxrelease`** | maximum amount that can be released if below spillway | m``^3`` s``^{-1}`` | - |
+| **`maxvolume`** | maximum storage (above which water is spilled) | m``^3`` | - |
+| **`targetfullfrac`** | target fraction full (of max storage)| -        | - |
+| **`targetminfrac`** | target minimum full fraction (of max storage) | -  | - |
+| `demandrelease`| minimum (environmental) flow released from reservoir  | m``^3`` s``^{-1}``| - |
+| `dt`             | model time step     | s | - |
+| `volume`             | volume     | m``^3`` | - |
+| `inflow`             | total inflow into reservoir | m``^3`` | - |
+| `outflow`             | outflow of reservoir | m``^3`` s``^{-1}`` | - |
+| `totaloutflow`        | total outflow of reservoir | m``^3`` | - |
+| `percfull`             | fraction full (of max storage) | - | - |
+| `precipitation`             | average precipitation for reservoir area | mm Δt⁻¹ | - |
+| `evaporation`             | average potential evaporation for reservoir area | mm Δt⁻¹ | - |
+| `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
+
+### Lakes
+The Table below shows the parameters (fields) of struct `Lake`, 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 input data (netCDF), and can be
+listed in the TOML configuration file under `[input.lateral.river.lake]`, to map the
+internal model parameter to the external netCDF variable.
+
+Two parameters lake coverage `areas` and the outlet of lakes (unique id) `locs` that are not
+part of the `Lake` struct are also required, and can be set as follows through the
+TOML file:
+
+```toml
+[input.lateral.river.lake]
+areas = "wflow_lakeareas"
+locs = "wflow_lakelocs"
+```
+
+The input parameter `linkedlakelocs` (listed under `[input.lateral.river.lake]`) is not
+equal to the internal model parameter `lowerlake_ind`, and is listed in the Table below
+between parentheses.
+
+|  parameter | description  	        | unit | default |
+|:---------------| --------------- | ---------------------- | ----- |
+| **`area`**         | area| m``^2`` | - |
+| **`b`**       | Rating curve coefficient  | - | - |
+| **`e`**      | Rating curve exponent | -  | - |
+| **`outflowfunc`**     | type of lake rating curve | - | - |
+| **`storfunc`** | type of lake storage curve| - | - |
+| **`threshold`**  | water level threshold ``H_0`` below that level outflow is zero | m  | - |
+| **`waterlevel`**  | waterlevel ``H`` of lake | m  | - |
+| **`lowerlake_ind`** (`linkedlakelocs`) | Index of lower lake (linked lakes) | - | 0 |
+| **`sh`** | data for storage curve | - | - |
+| **`hq`** | data rating curve | - | - |
+| `dt`             | model time step     | s |  - |
+| `inflow` | total inflow to the lake | m``^3``  | - |
+| `storage` | storage lake | m``^3``  | - |
+| `maxstorage`| maximum storage lake with rating curve type 1 | m``^3``  | - |
+| `outflow` | outflow lake | m``^3`` s``^{-1}``  | - |
+| `totaloutflow` | total outflow lake | m``^3``  | - |
+| `precipitation` | average precipitation for lake area | mm Δt⁻¹  | - |
+| `evaporation` | average potential evaporation for lake area | mm Δt⁻¹  | - |
+| `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
+
+
 ## Water allocation river
 The Table below shows the parameters (fields) of struct `AllocationRiver`, used when water
 demand and allocation is computed (optional), including a description of these parameters, the
diff --git a/docs/model_docs/parameters_vertical.qmd b/docs/model_docs/parameters_vertical.qmd
index 4823dae03..1c320d9ce 100644
--- a/docs/model_docs/parameters_vertical.qmd
+++ b/docs/model_docs/parameters_vertical.qmd
@@ -17,9 +17,9 @@ parameter `sl` is mapped as follows in the TOML file to the external netCDF vari
 specific_leaf = "Sl"
 ```
 
-Different [vertical hydraulic conductivity depth profiles](@ref soil): `exponential` (default),
-`exponential_constant`, `layered` and `layered_exponential` can be provided through the TOML
-file. Below an example for the `exponential_constant` profile:
+Different [vertical hydraulic conductivity depth profiles](./vertical/sbm.qmd#soil-processes):
+`exponential` (default), `exponential_constant`, `layered` and `layered_exponential` can be
+provided through the TOML file. Below an example for the `exponential_constant` profile:
 
 ```toml
 [input.vertical]
diff --git a/docs/model_docs/vertical/sbm.qmd b/docs/model_docs/vertical/sbm.qmd
index e553c9bd1..2d0953319 100644
--- a/docs/model_docs/vertical/sbm.qmd
+++ b/docs/model_docs/vertical/sbm.qmd
@@ -19,9 +19,9 @@ performed based on the air temperature. If the temperature is below a threshold
 (`tt`), precipitation will fall as snow. An interval parameter (`tti`) defines the range over
 which precipitation is partly falling as snow, and partly as rain. Snowfall is added to the
 snowpack, where it is subject to melting and refreezing (see the section on [snow and
-glaciers](@ref snow)). The amount of rainfall is subject to [interception](@ref interception),
-and ultimately becomes available for [evaporation](@ref evap) and/or [soil processes](@ref
-soil).
+glaciers](./shared_processes.qmd#snow-and-glaciers)). The amount of rainfall is subject to
+[interception](#rainfall-interception), and ultimately becomes available for
+[evaporation](#evaporation) and/or [soil processes](#soil-processes).
 
 ![Division between snow and precipitation based on the threshold temperature](../../images/snowfall.png)
 
@@ -308,8 +308,9 @@ satwaterdepth = satwaterdepth - actevapsat
 
 ## Snow and glaciers
 
-The snow and glacier model is described in [Snow and glaciers](@ref snow_and_glac). Both
-options can be enabled by specifying the following in the TOML file:
+The snow and glacier model is described in [Snow and
+glaciers](./shared_processes.qmd#snow-and-glaciers). Both options can be enabled by specifying
+the following in the TOML file:
 
 ```toml
 [model]
diff --git a/docs/model_docs/vertical/sediment.qmd b/docs/model_docs/vertical/sediment.qmd
index f1c9a1220..690a80d1d 100644
--- a/docs/model_docs/vertical/sediment.qmd
+++ b/docs/model_docs/vertical/sediment.qmd
@@ -172,7 +172,8 @@ slope angle.
 
 ## Delivery to the river system
 Once soil is detached, it can be transported by overland flow and reach the river system. This
-process is described in [Sediment Flux in overland flow](@ref).
+process is described in [Sediment Flux in overland
+flow](../lateral/sediment_flux.qmd#sediment-flux-in-overland-flow).
 
 ## References
 + D.B Beasley and L.F Huggins. ANSWERS - Users Manual. Technical report, EPA, 1991.
diff --git a/docs/user_guide/multi_threading.qmd b/docs/user_guide/multi_threading.qmd
index 6437ce6fa..99933b6fc 100644
--- a/docs/user_guide/multi_threading.qmd
+++ b/docs/user_guide/multi_threading.qmd
@@ -7,11 +7,14 @@ title: Multithreading
 Wflow supports multi-threading execution of the wflow\_sbm model that uses the kinematic wave
 approach for river, overland and lateral subsurface flow. Both the vertical SBM concept and the
 kinematic wave components of this model can run on multiple threads. The optional [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 the wflow\_sbm model [SBM + Groundwater flow](@ref config_sbm_gwf).
-The multi-threading functionality in wflow is considered experimental, see also the following
+inertial model for river
+flow](../model_docs/model_configurations.qmd#sbm--local-inertial-river-1d-and-land-2d) and the
+optional [local inertial model for river (1D) and land
+(2D)](../model_docs/model_configurations.qmd#sbm--local-inertial-river-1d-and-land-2d), both
+part of wflow\_sbm, can also run on multiple threads. The threading functionality for the
+kinematic wave may also be useful for the wflow\_sbm model [SBM + Groundwater
+flow](../model_docs/model_configurations.qmd#sbm--groundwater-flow). 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).
diff --git a/docs/user_guide/required_files.qmd b/docs/user_guide/required_files.qmd
index 05cf3986e..7ca047598 100644
--- a/docs/user_guide/required_files.qmd
+++ b/docs/user_guide/required_files.qmd
@@ -36,8 +36,9 @@ file is presented below.
 
 The list below contains a brief overview of several essential static maps required to run
 wflow. These NC variables names refer to the example data of the wflow\_sbm + kinematic wave
-model (see [here](@ref wflow_sbm_data)). Example data for the other model configurations can be
-found [here](@ref sample_data).
+model (see [here](../getting_started/download_example_models.qmd#wflow_sbm--kinematic-wave)).
+Example data for the other model configurations can be found
+[here](../getting_started/download_example_models.qmd).
 
 Description | NC variable name | unit
 --- | --- | ---
@@ -51,8 +52,8 @@ River slope | `RiverSlope` | m m$^{-1}$
 
 As mentioned before, the model parameters can also be defined as spatial maps. They can be
 included in the same netCDF file, as long as their variable names are correctly mapped in the
-TOML settings file. See the section on [example models](@ref sample_data) on how to use this
-functionality.
+TOML settings file. See the section on [example
+models](../getting_started/download_example_models.qmd) on how to use this functionality.
 
 ::: callout-important
 When using cyclic data, three different options are supported:
@@ -71,10 +72,12 @@ provided) or 166 (when 366 values are provided).
 The forcing data typically contains the meteorological boundary conditions. This data is
 provided as a single netCDF file, with several variables containing the forcing data for
 precipitation, temperature and potential evaporation. The code snippet below shows the contents
-of the example file (downloaded [here](@ref wflow_sbm_data)), and displaying the content with
-`NCDatasets` in Julia. As can be seen, each forcing variable (`precip`, `pet` and `temp`)
-consists of a three-dimensional dataset (`x`, `y`, and `time`), and each timestep consists of a
-two-dimensional map with values at each grid cell. Only values within the basin are required.
+of the example file (downloaded
+[here](../getting_started/download_example_models.qmd#wflow_sbm--kinematic-wave)), and
+displaying the content with `NCDatasets` in Julia. As can be seen, each forcing variable
+(`precip`, `pet` and `temp`) consists of a three-dimensional dataset (`x`, `y`, and `time`),
+and each timestep consists of a two-dimensional map with values at each grid cell. Only values
+within the basin are required.
 
 ::: callout-important
 Wflow expects the forcing to be "right-labelled". This means that e.g. daily precipitation at
diff --git a/docs/user_guide/toml_file.qmd b/docs/user_guide/toml_file.qmd
index 1c7211d6d..307ca16fa 100644
--- a/docs/user_guide/toml_file.qmd
+++ b/docs/user_guide/toml_file.qmd
@@ -53,7 +53,7 @@ debug-level messages. Note that for finer control, you can also pass an integer
 details, see Julia's
 [Logging](https://docs.julialang.org/en/v1/stdlib/Logging/#Log-event-structure) documentation.
 `path_log` sets the desired output path for the log file. For information on `fews_run`, see
-[Run from Delft-FEWS](@ref run_fews).
+[Run from Delft-FEWS](./fews.qmd).
 
 ## Model section
 Model-specific settings can be included in the model section of the TOML file.
@@ -77,7 +77,7 @@ output states of the model. This section is mostly relevant if the model needs t
 "warm" state (i.e. based on the results of a previous simulation). The example below shows how
 to save the output states of the current simulation, so it can be used to initialize another
 model in the future. Details on the settings required to start a model with a warm state can be
-found in the [additional model options](@ref reinit). If it is not required to store the
+found in the [additional model options](./warm_states.qmd). If it is not required to store the
 outstates of the current simulation, the entire `state` section can be removed.
 
 ```toml
@@ -238,21 +238,21 @@ In addition to gridded data, scalar data can also be written to a netCDF file. B
 example that shows how to write scalar data to the file "output\_scalar\_moselle.nc". For each
 netCDF variable, a `name` (external variable name) and a `parameter` (internal model parameter)
 are required. A `reducer` can be specified to apply to the model output. See more details in
-the [Output CSV section](@ref) section. If a `map` is provided to extract data for specific
-locations (e.g. `gauges`) or areas (e.g. `subcatchment`), the 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 `layer` dimension and are part of the
-vertical `sbm` concept, an internal layer index can be specified (an example is provided
-below). Similarly, for model parameters and variables that have an extra dimension `classes`
-and are part of the vertical `FLEXTopo` concept, the class name can be specified. If multiple
-layers or classes are desired, this can be specified in separate `[[netcdf.variable]]` entries.
-Note that the additional dimension should be specified when wflow is integrated with
+the [Output CSV section](#output-csv-section) section. If a `map` is provided to extract data
+for specific locations (e.g. `gauges`) or areas (e.g. `subcatchment`), the 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 `layer` dimension and
+are part of the vertical `sbm` concept, an internal layer index can be specified (an example is
+provided below). Similarly, for model parameters and variables that have an extra dimension
+`classes` and are part of the vertical `FLEXTopo` concept, the class name can be specified. If
+multiple layers or classes are desired, this can be specified in separate `[[netcdf.variable]]`
+entries. Note that the additional dimension should be specified 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, it
-is recommended to write scalar data to netCDF format, as the General Adapter of Delft-FEWS can
-directly ingest data from netCDF files. For more information, see [Run from Delft-FEWS](@ref
-run_fews).
+section](#output-csv-section), similar functionality is available for CSV. For integration with
+Delft-FEWS, it is recommended to write scalar data to netCDF format, as the General Adapter of
+Delft-FEWS can directly ingest data from netCDF files. For more information, see [Run from
+Delft-FEWS](./fews.qmd).
 
 ```toml
 [netcdf]