diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 0b1f22c54..cd78cf337 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -17,6 +17,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   lake storage. The variable `actevap` has also been added to the reservoir module.
 - The `set_states` function for model type `sbm` with local inertial routing for river and
   land component.
+- Inflow to reservoir and lake locations for local inertial routing: 1) with floodplain
+  routing, the floodplain discharge was not added to the inflow of these locations, 2) the
+  `to_river` variable from overland flow and lateral subsurface flow was not added to the
+  inflow of these locations.
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was
@@ -32,6 +36,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Checks to see if all states are covered in the .toml file. If not all states are covered,
   an error is thrown. If there are more states specified than required, these states are
   ignored (with a warning in the logging), and the simulation will continue.
+- Support different vertical hydraulic conductivity profiles for the `SBM` concept. See also
+  the following sections: [The SBM soil water accounting scheme](@ref) and [Subsurface flow
+  routing](@ref) for a short description.
 
 ## v0.7.3 - 2024-01-12
 
diff --git a/docs/src/images/sbm_ksat_profiles.png b/docs/src/images/sbm_ksat_profiles.png
new file mode 100644
index 000000000..6e6bbeda2
Binary files /dev/null and b/docs/src/images/sbm_ksat_profiles.png differ
diff --git a/docs/src/model_docs/lateral/kinwave.md b/docs/src/model_docs/lateral/kinwave.md
index 3e7b92727..1fdde084a 100644
--- a/docs/src/model_docs/lateral/kinwave.md
+++ b/docs/src/model_docs/lateral/kinwave.md
@@ -55,30 +55,64 @@ kinematic wave for surface water routing, as a cyclic parameter or as part of fo
 also [Input section](@ref)).
 
 ## Subsurface flow routing
-In the SBM model the kinematic wave approach is used to route subsurface flow laterally. The
-saturated store ``S`` can be drained laterally by saturated downslope subsurface flow per
-unit width of slope ``w`` [m] according to:
+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:
+- `exponential` (default)
+- `exponential_constant`
+- `layered`
+- `layered_exponential` 
+For the profiles `exponential` and `exponential_constant`, the saturated store ``S`` is
+drained laterally by saturated downslope subsurface flow for a slope with width ``w`` [m]
+according to:
 ```math
-    q=\frac{K_{0}\mathit{tan(\beta)}}{f}(e^{(-fz_{i})}-e^{(-fz_{t})})
+    Q = \begin{cases}
+    \frac{K_0\tan(\beta)}{f}\left(e^{(-fz_{i})}-e^{(-fz_\mathrm{exp})}\right) w + 
+    K_0e^{(-fz_\mathrm{exp})}(z_t-z_\mathrm{exp})\tan(\beta) w & \text{if $z_i < z_\mathrm{exp}$}\\
+    \\
+    K_0e^{(-fz_\mathrm{exp})}(z_t - z_i)\tan(\beta) w & \text{if $z_i \ge z_\mathrm{exp}$},
+    \end{cases}
 ```
-where ``\beta`` is element slope angle [deg.], ``q`` is subsurface flow [m``^{2}``/t],
-``K_{0}`` is the saturated hydraulic conductivity at the soil surface [m/t], ``z_{i}`` is
-the water table depth [m], ``z_{t}`` is total soil depth [m], and ``f`` is a scaling
-parameter [m``^{-1}``], that controls the decrease of vertical saturated conductivity with
-depth.
+where ``\beta`` is element slope angle, ``Q`` is subsurface flow [m``^{3}`` d``^{-1}``],
+``K_0`` is the saturated hydraulic conductivity at the soil surface [m d``^{-1}``], ``z_i``
+is the water table depth [m], ``z_{t}`` is the total soil depth [m], ``f`` is a scaling
+parameter [m``^{-1}``] that controls the decrease of ``K_0`` with depth and
+``z_\mathrm{exp}`` [m] is the depth from soil surface for which the exponential decline of
+``K_0`` is valid. For the `exponential` profile, ``z_\mathrm{exp}`` is equal to ``z_t``.
 
 Combining with the following continuity equation:
 ```math
-    (\theta_s-\theta_r)\frac{\partial h}{\partial t} = -w\frac{\partial q}{\partial x} + wr
+    (\theta_s-\theta_r)w\frac{\partial h}{\partial t} = -\frac{\partial Q}{\partial x} + wr
 ```
-where ``h`` is the water table height [m], ``x`` is the distance downslope [m], and ``r``
-is the net input rate [m/t] to the saturated store. Substituting for ``h (\frac{\partial
-q}{\partial h})``, gives:
+where ``h`` is the water table height [m], ``x`` is the distance downslope [m], and ``r`` is
+the net input rate [m d``^{-1}``] to the saturated store. Substituting for ``h
+(\frac{\partial Q}{\partial h})``, gives:
 ```math
-  w \frac{\partial q}{\partial t} = -cw\frac{\partial q}{\partial x} + cwr
+  \frac{\partial Q}{\partial t} = -c\frac{\partial Q}{\partial x} + cwr
 ```
 
-where celerity ``c = \frac{K_{0}\mathit{tan(\beta)}}{(\theta_s-\theta_r)} e^{(-fz_{i})}``
+where celerity ``c`` is calculated as follows:
+```math
+    c = \begin{cases}
+    \frac{K_0e^{(-fz_{i})}\tan(\beta)}{(\theta_s-\theta_r)} 
+    + \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)}  & \text{if $z_i < z_\mathrm{exp}$}\\
+    \\
+    \frac{K_0e^{(-fz_\mathrm{exp})}\tan(\beta)}{(\theta_s-\theta_r)} & \text{if $z_i \ge z_\mathrm{exp}$}.
+    \end{cases}
+```
+
+For the `layered` and `layered_exponential` profiles the equivalent horizontal hydraulic
+conductivity ``K_h`` [m d``^{-1}``] is calculated for water table height ``h = z_t-z_i``
+[m], and lateral subsurface flow is calculated as follows:
+```math
+  Q = K_h h \tan(\beta) w,
+```
+and celerity ``c`` is given by:
+```math
+    c = \frac{K_h \tan(\beta)}{(\theta_s-\theta_r)}.
+```
 
 The kinematic wave equation for lateral subsurface flow is solved iteratively using Newton's
 method.
diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 97234d044..7ad19afaa 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -155,11 +155,11 @@ between parentheses.
 The Table below shows the parameters (fields) of struct `LateralSSF`, 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). The
-soil related parameters `f`, `soilthickness`, `θₛ` and `θᵣ` are derived from the vertical
-`SBM` concept (including unit conversion for `f` and `soilthickness`), and can be listed in
-the TOML configuration file under `[input.vertical]`, to map the internal model parameter to
-the external netCDF variable. The internal slope model parameter `βₗ` is set through the
-TOML file as follows:
+soil related parameters `f`, `soilthickness`, `z_exp`, `θₛ` and `θᵣ` are derived from the
+vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and
+can be listed in the TOML configuration file under `[input.vertical]`, to map the internal
+model parameter to the external netCDF variable. The internal slope model parameter `βₗ` is
+set through the TOML file as follows:
 
 ```toml
 [input.lateral.land]
@@ -168,22 +168,33 @@ slope = "Slope"
 
 The parameter `kh₀` is computed by multiplying the vertical hydraulic conductivity at the
 soil surface `kv₀` (including unit conversion) of the vertical `SBM` concept with the
-external parameter `ksathorfrac` \[-\] (default value of 1.0). The external parameter
-`ksathorfrac` can be set as follows through the TOML file:
+internal parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter
+`khfrac` is set through the TOML file as follows:
 
 ```toml
 [input.lateral.subsurface]
 ksathorfrac = "KsatHorFrac"
 ```
 
-The `ksathorfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil
+The `khfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil
 core) 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₀` and `f` are used. For the
+`exponential_constant` profile `kh₀` 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
+`[input.vertical]`.
+
 | parameter | description  	        | unit | default |
 |:---------------| --------------- | ---------------------- | ----- |
 | `kh₀`          | horizontal hydraulic conductivity at soil surface  | m d``^{-1}`` |  3.0 |
-| **`f`** | a scaling parameter (controls exponential decline of kh₀) | m``^{-1}``  | 1.0 |
+| **`f`** | a scaling parameter (controls exponential decline of `kh₀`) | m``^{-1}``  | 1.0 |
+| `kh` | horizontal hydraulic conductivity | m d``^{-1}``  | - |
+| **`khfrac`** (`ksathorfrac`) | a muliplication factor applied to vertical hydraulic conductivity `kv` | -  | 100.0 |
 | **`soilthickness`** | soil thickness | m  | 2.0 |
 | **`θₛ`** | saturated water content (porosity) | -  | 0.6 |
 | **`θᵣ`** | residual water content  | -  | 0.01 |
@@ -192,13 +203,13 @@ scales (hillslope) in reality, not represented by the model resolution.
 | `dl` | drain length | m | - |
 | `dw` | drain width | m | - |
 | `zi` | pseudo-water table depth (top of the saturated zone) | m | - |
+| **`z_exp`** | depth from soil surface for which exponential decline of `kh₀` is valid | m | - |
 | `exfiltwater` | exfiltration (groundwater above surface level, saturated excess conditions) | m Δt⁻¹ | - |
 | `recharge` | net recharge to saturated store  | m``^2`` Δt⁻¹ | - |
 | `ssf` | subsurface flow  | m``^3`` d``{-1}``  | - |
 | `ssfin` | inflow from upstream cells | m``^3`` d``{-1}``  | - |
 | `ssfmax` | maximum subsurface flow | m``^2`` d``{-1}``  | - |
 | `to_river` | part of subsurface flow that flows to the river | m``^3`` d``{-1}``  | - |
-| `wb_pit` | boolean location (0 or 1) of a waterbody (wb, reservoir or lake) | -  | - |
 
 ## Local inertial
 
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 03a7817c0..fde92c15b 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -16,6 +16,20 @@ external netCDF variable `Sl`:
 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:
+
+```toml
+[input.vertical]
+ksat_profile = "exponential_constant"
+```
+
+For the `exponential` profile the input parameters `kv₀` and `f` are used. For the
+`exponential_constant` profile `kv₀` and `f` are used, and `z_exp` is required as input. For
+the `layered` profile, input parameter `kv` is used, and for the `layered_exponential`
+profile `kv` is used and `z_layered` is required as input.
+
 |  parameter | description    | unit | default |
 |:---------------| --------------- | ---------------------- | ----- |
 | **`cfmax`**  | degree-day factor | mm ᵒC``^{-1}`` Δt``^{-1}`` | 3.75653 mm ᵒC``^{-1}`` day``^{-1}``  |
@@ -33,7 +47,10 @@ specific_leaf = "Sl"
 | **`θₛ`** (`theta_s`) | saturated water content (porosity) | - | 0.6  |
 | **`θᵣ`** (`theta_r`) | residual water content | - | 0.01  |
 | **`kv₀`** (`kv_0`) | Vertical hydraulic conductivity at soil surface | mm Δt``^{-1}`` | 3000.0 mm day``^{-1}``|
-| **`f`** | scaling parameter (controls exponential decline of kv₀) | mm``^{-1}`` | 0.001  |
+| **`kv`** |  Vertical hydraulic conductivity per soil layer | mm Δt``^{-1}`` | 1000.0 mm day``^{-1}``|
+| **`f`** | scaling parameter (controls exponential decline of `kv₀`) | mm``^{-1}`` | 0.001  |
+| **`z_exp`** | Depth from soil surface for which exponential decline of `kv₀` is valid | mm | -  |
+| **`z_layered`** | Depth from soil surface for which layered profile (of `layered_exponential`) is valid | mm | -  |
 | **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | 10.0  |
 | **`soilthickness`** | soil thickness | mm | 2000.0  |
 | **`infiltcappath`** | infiltration capacity of the compacted areas  | mm Δt``^{-1}`` | 10.0 mm day``^{-1}`` |
@@ -59,6 +76,7 @@ specific_leaf = "Sl"
 | `n`      |  number of grid cells    | - | - |
 | `nlayers`      |  number of soil layers    | - | - |
 | `n_unsatlayers`      |  number of unsaturated soil layers    | - | - |
+| `nlayers_kv`      |  number of soil layers with vertical hydraulic conductivity value `kv` | - | - |
 | `riverfrac`      |  fraction of river    | - | - |
 | `act_thickl`      |  thickness of soil layers    | mm | - |
 | `sumlayers`      |  cumulative sum of soil layers thickness, starting at soil surface | mm | - |
@@ -69,7 +87,6 @@ specific_leaf = "Sl"
 | `zi`|  pseudo-water table depth (top of the saturated zone) | mm | - |
 | `soilwatercapacity`|  soilwater capacity | mm | - |
 | `canopystorage`|  canopy storage | mm | - |
-| `canopygapfraction` | canopygapfraction | - | - |
 |**`precipitation`** | precipitation | mm Δt``^{-1}``|  - |
 | **`temperature`** | temperature | ᵒC | - |
 | **`potential_evaporation`** | potential evaporation | mm Δt``^{-1}`` | - |
@@ -111,7 +128,6 @@ specific_leaf = "Sl"
 | `snow` | snow storage  | mm | - |
 | `snowwater` | liquid water content in the snow pack  | mm | - |
 | `rainfallplusmelt` | snowmelt + precipitation as rainfall | mm Δt``^{-1}`` | - |
-| `glacierstore` | water within the glacier | mm | - |
 | `tsoil` | top soil temperature | ᵒC | - |
 | **`leaf_area_index`** | leaf area index | m``^2`` m``{-2}`` | - |
 | `waterlevel_land` | water level land | mm | - |
diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md
index 14b11338d..7f174a441 100644
--- a/docs/src/model_docs/vertical/sbm.md
+++ b/docs/src/model_docs/vertical/sbm.md
@@ -350,14 +350,28 @@ t``^{-1}``]) is in that case controlled by the saturated hydraulic conductivity
 ```math
     st=K_{\mathit{sat}}\frac{U_{s}}{S_{d}}
 ```
-Saturated conductivity (``K_{sat}`` [mm t``^{-1}``]) declines with soil depth (``z`` [mm])
-in the model according to:
+
+Four different saturated hydraulic conductivity depth profiles (`ksat_profile`) are
+available and a `ksat_profile` can be specified in the TOML file as follows:
+
+```toml
+[input.vertical]
+ksat_profile = "exponential_constant" # optional, one of ("exponential", "exponential_constant", "layered", "layered_exponential"), default is "exponential"
+```
+
+Soil measurements are often available for about the upper 1.5-2 m of the soil column to
+estimate the saturated hydraulic conductivity, while these measurements are often lacking
+for soil depths beyond 1.5-2 m. These different profiles allow to extent the saturated
+hydraulic conductivity profile based on measurements (either an exponential fit or hydraulic
+conductivity value per soil layer) with an exponential or constant profile. By default, with
+`ksat_profile` "exponential", the saturated hydraulic conductivity (``K_{sat}`` [mm
+t``^{-1}``]) declines with soil depth (``z`` [mm]) in the model according to:
 
 ```math
-    K_{sat}=K_{0}e^{(-fz)}
+    K_{sat}=K_{0}e^{(-fz)},
 ```
-where ``K_{0}`` [mm t``^{-1}``] is the saturated conductivity at the soil surface and ``f``
-is a scaling parameter [mm``^{-1}``].
+where ``K_{0}`` [mm t``^{-1}``] is the saturated hydraulic conductivity at the soil surface
+and ``f`` is a scaling parameter [mm``^{-1}``].
 
 The plot below shows the relation between soil depth ``z`` and saturated hydraulic
 conductivity ``K_{sat}`` for different values of ``f``.
@@ -385,6 +399,30 @@ conductivity ``K_{sat}`` for different values of ``f``.
     end                                                                                     # hide
 ```
 
+With `ksat_profile` "exponential\_constant", ``K_{sat}`` declines exponentially with soil
+depth ``z`` until ``z_\mathrm{exp}`` [mm] below the soil surface, and stays constant at and
+beyond soil depth ``z_\mathrm{exp}``:
+
+```math
+    K_{sat} = \begin{cases}
+    K_{0}e^{(-fz)} & \text{if $z < z_\mathrm{exp}$}\\
+    K_{0}e^{(-fz_\mathrm{exp})} & \text{if $z \ge z_\mathrm{exp}$}.
+    \end{cases}
+```
+
+It is also possible to provide a ``K_{sat}`` value per soil layer by specifying
+`ksat_profile` "layered", these ``K_{sat}`` values are used directly to compute the vertical
+transfer of water between soil layers and to the saturated store ``S``. Finally, with the
+`ksat_profile` "layered\_exponential" a ``K_{sat}`` value per soil layer is used until depth
+``z_\mathrm{layered}`` below the soil surface, and beyond ``z_\mathrm{layered}`` an
+exponential decline of ``K_{sat}`` (of the soil layer with bottom ``z_\mathrm{layered}``)
+controlled by ``f`` occurs. The different available `ksat_profle` options are schematized in
+the figure below where the blue line represents the ``K_{sat}`` value.
+
+![ksat_profiles](../../images/sbm_ksat_profiles.png)
+
+*Overview of available `ksat_profile` options, for a soil column with five layers*
+
 ### Infiltration
 
 The water available for infiltration is taken as the rainfall including meltwater.
diff --git a/server/test/client.jl b/server/test/client.jl
index 41fc1182d..b6e000ac8 100644
--- a/server/test/client.jl
+++ b/server/test/client.jl
@@ -33,15 +33,15 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 181)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 181)
-    @test request((fn = "get_input_var_names",))["input_var_names"][[1, 5, 151, 175]] == [
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 192)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 192)
+    @test request((fn = "get_input_var_names",))["input_var_names"][[1, 6, 161, 186]] == [
         "vertical.nlayers",
         "vertical.θᵣ",
         "lateral.river.q",
         "lateral.river.reservoir.outflow",
     ]
-    @test request((fn = "get_output_var_names",))["output_var_names"][[1, 5, 151, 175]] == [
+    @test request((fn = "get_output_var_names",))["output_var_names"][[1, 6, 161, 186]] == [
         "vertical.nlayers",
         "vertical.θᵣ",
         "lateral.river.q",
diff --git a/src/flow.jl b/src/flow.jl
index c9c8d25e8..a07f15991 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -13,6 +13,7 @@ abstract type SurfaceFlow end
     qlat::Vector{T} | "m2 s-1"                   # Lateral inflow per unit length [m² s⁻¹]
     inwater::Vector{T} | "m3 s-1"                # Lateral inflow [m³ s⁻¹]
     inflow::Vector{T} | "m3 s-1"                 # External inflow (abstraction/supply/demand) [m³ s⁻¹]
+    inflow_wb::Vector{T} | "m3 s-1"              # inflow waterbody (lake or reservoir model) from land part [m³ s⁻¹]
     volume::Vector{T} | "m3"                     # Kinematic wave volume [m³] (based on water level h)
     h::Vector{T} | "m"                           # Water level [m]
     h_av::Vector{T} | "m"                        # Average water level [m]
@@ -158,6 +159,7 @@ function initialize_surfaceflow_river(
         qlat = zeros(Float, n),
         inwater = zeros(Float, n),
         inflow = zeros(Float, n),
+        inflow_wb = zeros(Float, n),
         volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
@@ -246,7 +248,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
     sf.volume .= sf.dl .* sf.width .* sf.h
 end
 
-function update(sf::SurfaceFlowRiver, network, inflow_wb, doy)
+function update(sf::SurfaceFlowRiver, network, doy)
 
     @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
         network
@@ -306,7 +308,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy)
                         # run reservoir model and copy reservoir outflow to inflow (qin) of
                         # downstream river cell
                         i = sf.reservoir_index[v]
-                        update(sf.reservoir, i, sf.q[v] + inflow_wb[v], Δt)
+                        update(sf.reservoir, i, sf.q[v] + sf.inflow_wb[v], Δt)
 
                         downstream_nodes = outneighbors(graph, v)
                         n_downstream = length(downstream_nodes)
@@ -327,7 +329,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy)
                         # run lake model and copy lake outflow to inflow (qin) of downstream river
                         # cell
                         i = sf.lake_index[v]
-                        update(sf.lake, i, sf.q[v] + inflow_wb[v], doy, Δt)
+                        update(sf.lake, i, sf.q[v] + sf.inflow_wb[v], doy, Δt)
 
                         downstream_nodes = outneighbors(graph, v)
                         n_downstream = length(downstream_nodes)
@@ -389,6 +391,8 @@ end
 @get_units @exchange @grid_type @grid_location @with_kw struct LateralSSF{T}
     kh₀::Vector{T} | "m d-1"               # Horizontal hydraulic conductivity at soil surface [m d⁻¹]
     f::Vector{T} | "m-1"                   # A scaling parameter [m⁻¹] (controls exponential decline of kh₀)
+    kh::Vector{T} | "m d-1"                # Horizontal hydraulic conductivity [m d⁻¹]
+    khfrac::Vector{T} | "-"                # A muliplication factor applied to vertical hydraulic conductivity `kv` [-]
     soilthickness::Vector{T} | "m"         # Soil thickness [m]
     θₛ::Vector{T} | "-"                     # Saturated water content (porosity) [-]
     θᵣ::Vector{T} | "-"                    # Residual water content [-]
@@ -397,6 +401,7 @@ end
     dl::Vector{T} | "m"                    # Drain length [m]
     dw::Vector{T} | "m"                    # Flow width [m]
     zi::Vector{T} | "m"                    # Pseudo-water table depth [m] (top of the saturated zone)
+    z_exp::Vector{T} | "m"                 # Depth [m] from soil surface for which exponential decline of kv₀ is valid
     exfiltwater::Vector{T} | "m Δt-1"      # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions)
     recharge::Vector{T} | "m2 Δt-1"        # Net recharge to saturated store [m² Δt⁻¹]
     ssf::Vector{T} | "m3 d-1"              # Subsurface flow [m³ d⁻¹]
@@ -411,7 +416,7 @@ end
 end
 
 
-function update(ssf::LateralSSF, network, frac_toriver)
+function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
     @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network
 
     ns = length(subdomain_order)
@@ -433,22 +438,42 @@ function update(ssf::LateralSSF, network, frac_toriver)
                     upstream_nodes[n],
                     eltype(ssf.to_river),
                 )
-
-                ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
-                    ssf.ssfin[v],
-                    ssf.ssf[v],
-                    ssf.zi[v],
-                    ssf.recharge[v],
-                    ssf.kh₀[v],
-                    ssf.βₗ[v],
-                    ssf.θₛ[v] - ssf.θᵣ[v],
-                    ssf.f[v],
-                    ssf.soilthickness[v],
-                    ssf.Δt,
-                    ssf.dl[v],
-                    ssf.dw[v],
-                    ssf.ssfmax[v],
-                )
+                if (ksat_profile == "exponential") ||
+                   (ksat_profile == "exponential_constant")
+                    ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
+                        ssf.ssfin[v],
+                        ssf.ssf[v],
+                        ssf.zi[v],
+                        ssf.recharge[v],
+                        ssf.kh₀[v],
+                        ssf.βₗ[v],
+                        ssf.θₛ[v] - ssf.θᵣ[v],
+                        ssf.f[v],
+                        ssf.soilthickness[v],
+                        ssf.Δt,
+                        ssf.dl[v],
+                        ssf.dw[v],
+                        ssf.ssfmax[v],
+                        ssf.z_exp[v],
+                        ksat_profile,
+                    )
+                elseif (ksat_profile == "layered") ||
+                       (ksat_profile == "layered_exponential")
+                    ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
+                        ssf.ssfin[v],
+                        ssf.ssf[v],
+                        ssf.zi[v],
+                        ssf.recharge[v],
+                        ssf.kh[v],
+                        ssf.βₗ[v],
+                        ssf.θₛ[v] - ssf.θᵣ[v],
+                        ssf.soilthickness[v],
+                        ssf.Δt,
+                        ssf.dl[v],
+                        ssf.dw[v],
+                        ssf.ssfmax[v],
+                    )
+                end
             end
         end
     end
@@ -668,14 +693,16 @@ function initialize_shallowwater_river(
     return sw_river, nodes_at_link
 end
 
-function shallowwater_river_update(
-    sw::ShallowWaterRiver,
-    network,
-    Δt,
-    inflow_wb,
-    doy,
-    update_h,
-)
+"Return the upstream inflow for a waterbody in `ShallowWaterRiver`"
+function get_inflow_waterbody(sw::ShallowWaterRiver, src_edge)
+    q_in = sum_at(sw.q, src_edge)
+    if !isnothing(sw.floodplain)
+        q_in = q_in + sum_at(sw.floodplain.q, src_edge)
+    end
+    return q_in
+end
+
+function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, update_h)
 
     @unpack nodes_at_link, links_at_node = network
 
@@ -683,32 +710,6 @@ function shallowwater_river_update(
     if !isnothing(sw.floodplain)
         sw.floodplain.q0 .= sw.floodplain.q
     end
-    # For reservoir and lake locations the local inertial solution is replaced by the
-    # reservoir or lake model. These locations are handled as boundary conditions in the
-    # local inertial model (fixed h).
-    for v in eachindex(sw.reservoir_index)
-        i = sw.reservoir_index[v]
-        update(
-            sw.reservoir,
-            v,
-            sum_at(sw.q0, links_at_node.src[i]) + inflow_wb[i] + sw.inflow_wb[i],
-            Δt,
-        )
-        sw.q[i] = sw.reservoir.outflow[v]
-        sw.q_av[i] += sw.q[i] * Δt
-    end
-    for v in eachindex(sw.lake_index)
-        i = sw.lake_index[v]
-        update(
-            sw.lake,
-            v,
-            sum_at(sw.q0, links_at_node.src[i]) + inflow_wb[i] + sw.inflow_wb[i],
-            doy,
-            Δt,
-        )
-        sw.q[i] = sw.lake.outflow[v]
-        sw.q_av[i] += sw.q[i] * Δt
-    end
     @tturbo for j in eachindex(sw.active_e)
         i = sw.active_e[j]
         i_src = nodes_at_link.src[i]
@@ -838,6 +839,25 @@ function shallowwater_river_update(
             sw.floodplain.q_av[i] += sw.floodplain.q[i] * Δt
         end
     end
+    # For reservoir and lake locations the local inertial solution is replaced by the
+    # reservoir or lake model. These locations are handled as boundary conditions in the
+    # local inertial model (fixed h).
+    for v in eachindex(sw.reservoir_index)
+        i = sw.reservoir_index[v]
+
+        q_in = get_inflow_waterbody(sw, links_at_node.src[i])
+        update(sw.reservoir, v, q_in + sw.inflow_wb[i], Δt)
+        sw.q[i] = sw.reservoir.outflow[v]
+        sw.q_av[i] += sw.q[i] * Δt
+    end
+    for v in eachindex(sw.lake_index)
+        i = sw.lake_index[v]
+
+        q_in = get_inflow_waterbody(sw, links_at_node.src[i])
+        update(sw.lake, v, q_in + sw.inflow_wb[i], doy, Δt)
+        sw.q[i] = sw.lake.outflow[v]
+        sw.q_av[i] += sw.q[i] * Δt
+    end
     if update_h
         @batch per = thread minbatch = 2000 for i in sw.active_n
 
@@ -885,13 +905,7 @@ function shallowwater_river_update(
     end
 end
 
-function update(
-    sw::ShallowWaterRiver{T},
-    network,
-    inflow_wb,
-    doy;
-    update_h = true,
-) where {T}
+function update(sw::ShallowWaterRiver{T}, network, doy; update_h = true) where {T}
     @unpack nodes_at_link, links_at_node = network
 
     if !isnothing(sw.reservoir)
@@ -917,7 +931,7 @@ function update(
         if t + Δt > sw.Δt
             Δt = sw.Δt - t
         end
-        shallowwater_river_update(sw, network, Δt, inflow_wb, doy, update_h)
+        shallowwater_river_update(sw, network, Δt, doy, update_h)
         t = t + Δt
     end
     sw.q_av ./= sw.Δt
@@ -967,6 +981,7 @@ const dirs = (:yd, :xd, :xu, :yu)
     volume::Vector{T} | "m3"                                # total volume of cell (including river volume for river cells)
     error::Vector{T} | "m3"                                 # error volume
     runoff::Vector{T} | "m3 s-1"                            # runoff from hydrological model
+    inflow_wb::Vector{T} | "m3 s-1"                         # inflow to water body from hydrological model
     h::Vector{T} | "m"                                      # water depth of cell (for river cells the reference is the river bed elevation `zb`)
     z::Vector{T} | "m"                                      # elevation of cell
     froude_limit::Bool | "-" | 0 | "none" | "none"          # if true a check is performed if froude number > 1.0 (algorithm is modified)
@@ -1083,6 +1098,7 @@ function initialize_shallowwater_land(
         volume = zeros(n),
         error = zeros(n),
         runoff = zeros(n),
+        inflow_wb = zeros(n),
         h = zeros(n),
         h_av = zeros(n),
         z = elevation,
@@ -1129,7 +1145,6 @@ function update(
     sw::ShallowWaterLand{T},
     swr::ShallowWaterRiver{T},
     network,
-    inflow_wb,
     doy;
     update_h = false,
 ) where {T}
@@ -1158,8 +1173,8 @@ function update(
         if t + Δt > swr.Δt
             Δt = swr.Δt - t
         end
-        shallowwater_river_update(swr, network.river, Δt, inflow_wb, doy, update_h)
-        update(sw, swr, network, Δt)
+        shallowwater_river_update(swr, network.river, Δt, doy, update_h)
+        shallowwater_update(sw, swr, network, Δt)
         t = t + Δt
     end
     swr.q_av ./= swr.Δt
@@ -1167,7 +1182,12 @@ function update(
     sw.h_av ./= sw.Δt
 end
 
-function update(sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, Δt) where {T}
+function shallowwater_update(
+    sw::ShallowWaterLand{T},
+    swr::ShallowWaterRiver{T},
+    network,
+    Δt,
+) where {T}
 
     indices = network.land.staggered_indices
     inds_riv = network.land.index_river
@@ -1273,7 +1293,9 @@ function update(sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, Δt
                 # for reservoir or lake set inflow from land part, these are boundary points
                 # and update of volume and h is not required
                 swr.inflow_wb[inds_riv[i]] =
-                    sw.runoff[i] + (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i])
+                    sw.inflow_wb[i] +
+                    sw.runoff[i] +
+                    (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i])
             else
                 sw.volume[i] +=
                     (
@@ -1564,7 +1586,7 @@ end
 """
     set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<:SBM,R,W,T}
 
-Set `inwater` of the river component for a `Model` with vertical `SBM` concept.
+Set `inwater` of the lateral river component for a `Model` with vertical `SBM` concept.
 `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}
@@ -1589,7 +1611,7 @@ end
 """
     set_river_inwater(model, ssf_toriver)
 
-Set `inwater` of the river component (based on overland flow).
+Set `inwater` of the lateral river component (based on overland flow).
 """
 function set_river_inwater(model, ssf_toriver)
     @unpack lateral, network = model
@@ -1600,7 +1622,7 @@ end
 """
     set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
 
-Set `inwater` of the land component for the `SbmGwgModel` type.
+Set `inwater` of the lateral land component for the `SbmGwgModel` type.
 """
 function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     @unpack lateral, vertical, network, config = model
@@ -1620,7 +1642,7 @@ end
 """
     set_land_inwater(model)
 
-Set `inwater` of the land component, based on `runoff` of the `vertical` concept.
+Set `inwater` of the lateral land component, based on `runoff` of the `vertical` concept.
 """
 function set_land_inwater(model)
     @unpack lateral, vertical, network = model
@@ -1628,66 +1650,90 @@ function set_land_inwater(model)
         (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
 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
+# the river kinematic wave (kinematic wave is also solved for the water body cell). For
+# local inertial river routing, `to_river` is included, because for the local inertial
+# solution the water body cells are excluded (boundary condition). For `GroundwaterFlow`
+# (Darcian flow in 4 directions), the lateral subsurface flow is excluded (for now) and
+# inflow consists of overland flow.
 """
-    get_inflow_waterbody(model)
+    set_inflow_waterbody(
+        model::Model{N,L,V,R,W,T},
+    ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,SurfaceFlow}},V,R,W,T}
 
-Get inflow to a water body (reservoir or lake) `inflow_wb` based on overland flow.
+Set inflow from the subsurface and land components to a water body (reservoir or lake)
+`inflow_wb` from a model type that contains the lateral components `SurfaceFlow`.
 """
-function get_inflow_waterbody(model)
+function set_inflow_waterbody(
+    model::Model{N,L,V,R,W,T},
+) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,SurfaceFlow}},V,R,W,T}
     @unpack lateral, network = model
+    @unpack subsurface, land, river = lateral
     inds = network.index_river
+
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
-        inflow_wb = lateral.land.q_av[inds]
-    else
-        inflow_wb = nothing
+        if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange
+            @. river.inflow_wb =
+                subsurface.ssf[inds] / tosecond(basetimestep) + land.q_av[inds]
+        elseif typof(subsurface.flow) <: GroundwaterFlow || isnothing(subsurface)
+            river.inflow_wb .= land.q_av[inds]
+        end
     end
-    return inflow_wb
 end
 
 """
-    get_inflow_waterbody(
+    set_inflow_waterbody(
         model::Model{N,L,V,R,W,T},
-    ) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,SurfaceFlow,Any}},V,R,W,T}
+    ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,ShallowWaterRiver}},V,R,W,T}
 
-Get inflow to a water body (reservoir or lake) `inflow_wb` from a model type that contains
-the lateral components `LateralSSF` and `SurfaceFlow`.
+Set inflow from the subsurface and land components to a water body (reservoir or lake)
+`inflow_wb` from a model type that contains the lateral components `SurfaceFlow` and
+`ShallowWaterRiver`.
 """
-function get_inflow_waterbody(
+function set_inflow_waterbody(
     model::Model{N,L,V,R,W,T},
-) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,SurfaceFlow,Any}},V,R,W,T}
+) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,SurfaceFlow,ShallowWaterRiver}},V,R,W,T}
     @unpack lateral, network = model
-
+    @unpack subsurface, land, river = lateral
     inds = network.index_river
+
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
-        inflow_wb =
-            lateral.subsurface.ssf[inds] ./ tosecond(basetimestep) .+
-            lateral.land.q_av[inds]
-    else
-        inflow_wb = nothing
+        if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange
+            @. river.inflow_wb =
+                (subsurface.ssf[inds] + subsurface.to_river[inds]) /
+                tosecond(basetimestep) +
+                land.q_av[inds] +
+                land.to_river[inds]
+        elseif typeof(subsurface.flow) <: GroundwaterFlow || isnothing(subsurface)
+            @. river.inflow_wb = lateral.land.q_av[inds] + lateral.land.to_river[inds]
+        end
     end
-    return inflow_wb
 end
 
 """
-    get_inflow_waterbody(
+    set_inflow_waterbody(
         model::Model{N,L,V,R,W,T},
-    ) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,ShallowWaterLand,Any}},V,R,W,T}
+    ) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T}
 
-Get inflow to a water body (reservoir or lake) `inflow_wb` from a model type that contains
-the lateral components `LateralSSF` and `ShallowWaterLand`.
+Set inflow from the subsurface and land components to a water body (reservoir or lake)
+`inflow_wb` from a model type that contains the lateral components `ShallowWaterLand` and
+`ShallowWaterRiver`.
 """
-function get_inflow_waterbody(
+function set_inflow_waterbody(
     model::Model{N,L,V,R,W,T},
-) where {N,L<:NamedTuple{<:Any,<:Tuple{LateralSSF,ShallowWaterLand,Any}},V,R,W,T}
+) where {N,L<:NamedTuple{<:Any,<:Tuple{Any,ShallowWaterLand,ShallowWaterRiver}},V,R,W,T}
     @unpack lateral, network = model
-
+    @unpack subsurface, land, river = lateral
     inds = network.index_river
+
     if !isnothing(lateral.river.reservoir) || !isnothing(lateral.river.lake)
-        inflow_wb = lateral.subsurface.ssf[inds] ./ tosecond(basetimestep)
-    else
-        inflow_wb = nothing
+        if typeof(subsurface) <: LateralSSF || typeof(subsurface) <: GroundwaterExchange
+            @. land.inflow_wb[inds] =
+                (subsurface.ssf[inds] + subsurface.to_river[inds]) / tosecond(basetimestep)
+        end
     end
-    return inflow_wb
 end
 
 """
@@ -1705,12 +1751,8 @@ function surface_routing(model; ssf_toriver = 0.0)
 
     # run river flow
     set_river_inwater(model, ssf_toriver)
-    update(
-        lateral.river,
-        network.river,
-        get_inflow_waterbody(model),
-        julian_day(clock.time - clock.Δt),
-    )
+    set_inflow_waterbody(model)
+    update(lateral.river, network.river, julian_day(clock.time - clock.Δt))
 end
 
 """
@@ -1738,12 +1780,6 @@ function surface_routing(
             vertical.Δt
         )
     )
-
-    update(
-        lateral.land,
-        lateral.river,
-        network,
-        get_inflow_waterbody(model),
-        julian_day(clock.time - clock.Δt),
-    )
+    set_inflow_waterbody(model)
+    update(lateral.land, lateral.river, network, julian_day(clock.time - clock.Δt))
 end
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index b61585608..380edd39e 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -68,8 +68,65 @@ function kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt)
     return Q
 end
 
-"Kinematic wave for lateral subsurface flow for a single cell and timestep"
-function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, dw, ssfmax)
+"Returns water table depth `zi` based on lateral subsurface flow `ssf` and hydraulic conductivity profile `ksat_profile`"
+function ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile)
+    if ksat_profile == "exponential"
+        zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f
+    elseif ksat_profile == "exponential_constant"
+        ssf_constant = kh₀ * β * exp(-f * z_exp) * (d - z_exp) * dw
+        if ssf > ssf_constant
+            zi = log((f * (ssf - ssf_constant)) / (dw * kh₀ * β) + exp(-f * z_exp)) / -f
+        else
+            zi = d - ssf / (dw * kh₀ * β * exp(-f * z_exp))
+        end
+    end
+    return zi
+end
+
+"Returns kinematic wave celecity `Cn` of lateral subsurface flow based on hydraulic conductivity profile `ksat_profile`"
+function ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile)
+    if ksat_profile == "exponential"
+        Cn = (kh₀ * exp(-f * zi) * β) / θₑ
+    elseif ksat_profile == "exponential_constant"
+        Cn_const = (kh₀ * exp(-f * z_exp) * β) / θₑ
+        if zi < z_exp
+            Cn = (kh₀ * exp(-f * zi) * β) / θₑ + Cn_const
+        else
+            Cn = Cn_const
+        end
+    end
+    return Cn
+end
+
+"""
+    kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, dw, ssfmax, z_exp, ksat_profile)
+
+Kinematic wave for lateral subsurface flow for a single cell and timestep. An exponential
+decline of hydraulic conductivity at the soil surface `kh₀`, controllled by parameter `f`,
+is assumed. The hydraulic conductivity profile `ksat_profile` is either `exponential` or
+`exponential_constant`, with `z_exp` the depth from the soil surface for which the
+exponential decline of `kh₀` is valid.
+
+Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate
+`exfilt`.
+"""
+function kinematic_wave_ssf(
+    ssfin,
+    ssfₜ₋₁,
+    ziₜ₋₁,
+    r,
+    kh₀,
+    β,
+    θₑ,
+    f,
+    d,
+    Δt,
+    Δx,
+    dw,
+    ssfmax,
+    z_exp,
+    ksat_profile,
+)
 
     ϵ = 1.0e-12
     max_iters = 3000
@@ -81,10 +138,10 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ
         ssf = (ssfₜ₋₁ + ssfin) / 2.0
         count = 1
 
-        # Estimate zi on the basis of the relation between subsurfacel flow and zi
-        zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f
+        # Estimate zi on the basis of the relation between subsurface flow and zi
+        zi = ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile)
         # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation)
-        Cn = (kh₀ * exp(-f * zi) * β) / θₑ
+        Cn = ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile)
         # Term of the continuity equation for Newton-Raphson iteration for iteration 1
         # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
         # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
@@ -104,10 +161,9 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ
         # Start while loop of Newton-Raphson iteration m until continuity equation approaches zero
         while true
             # Estimate zi on the basis of the relation between lateral flow rate and groundwater level
-            zi = log((f * ssf) / (dw * kh₀ * β) + exp(-f * d)) / -f
+            zi = ssf_water_table_depth(ssf, kh₀, β, f, d, dw, z_exp, ksat_profile)
             # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation
-            Cn = (kh₀ * exp(-f * zi) * β) / θₑ
-
+            Cn = ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile)
             # Term of the continuity equation for given Newton-Raphson iteration m
             # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
             # then (1./Cn)*ssfₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
@@ -142,6 +198,71 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θ
     end
 end
 
+"""
+    kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, dw, ssfmax)
+
+Kinematic wave for lateral subsurface flow for a single cell and timestep, based on
+(average) hydraulic conductivity `kh`.
+
+Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate
+`exfilt`.
+"""
+function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, dw, ssfmax)
+
+    ϵ = 1.0e-12
+    max_iters = 3000
+
+    if ssfin + ssfₜ₋₁ ≈ 0.0 && r <= 0.0
+        return 0.0, d, 0.0
+    else
+        # initial estimate
+        ssf = (ssfₜ₋₁ + ssfin) / 2.0
+        count = 1
+        # celerity (Cn)
+        Cn = (β * kh) / θₑ
+        # constant term of the continuity equation for Newton-Raphson
+        c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssfₜ₋₁ + r
+        # continuity equation of which solution should be zero
+        fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
+        # Derivative of the continuity equation
+        dfQ = (Δt / Δx) + 1.0 / Cn
+        # Update lateral subsurface flow estimate ssf
+        ssf = ssf - (fQ / dfQ)
+        if isnan(ssf)
+            ssf = 0.0
+        end
+        ssf = max(ssf, 1.0e-30)
+
+        # Start while loop of Newton-Raphson iteration
+        while true
+            fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
+            dfQ = (Δt / Δx) + 1.0 / Cn
+            ssf = ssf - (fQ / dfQ)
+            if isnan(ssf)
+                ssf = 0.0
+            end
+            ssf = max(ssf, 1.0e-30)
+            if (abs(fQ) <= ϵ) || (count >= max_iters)
+                break
+            end
+            count += 1
+        end
+
+        # Constrain the lateral subsurface flow rate ssf
+        ssf = min(ssf, (ssfmax * dw))
+        # On the basis of the lateral flow rate, estimate the amount of groundwater level above surface (saturation excess conditions), then rest = negative
+        zi = ziₜ₋₁ - (ssfin * Δt + r * Δx - ssf * Δt) / (dw * Δx) / θₑ
+        if zi > d
+            ssf = max(ssf - (dw * Δx) * θₑ * (zi - d), 1.0e-30)
+        end
+        # Exfiltration rate and set zi to zero.
+        exfilt = min(zi, 0.0) * -θₑ
+        zi = clamp(zi, 0.0, d)
+
+        return ssf, zi, exfilt
+    end
+end
+
 """
     accucapacitystate!(material, network, capacity) -> material
 
diff --git a/src/sbm.jl b/src/sbm.jl
index 77c44347c..adc103df4 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -9,6 +9,8 @@
     nlayers::Vector{Int} | "-"
     # Number of unsaturated soil layers
     n_unsatlayers::Vector{Int} | "-"
+    # Number of soil layers with vertical hydraulic conductivity value `kv`
+    nlayers_kv::Vector{Int} | "-"
     # Fraction of river [-]
     riverfrac::Vector{T} | "-"
     # Saturated water content (porosity) [-]
@@ -17,6 +19,8 @@
     θᵣ::Vector{T} | "-"
     # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface
     kv₀::Vector{T}
+    # Vertical hydraulic conductivity [mm Δt⁻¹] per soil layer
+    kv::Vector{SVector{N,T}} | "-"
     # Muliplication factor [-] applied to kv_z (vertical flow)
     kvfrac::Vector{SVector{N,T}} | "-"
     # Air entry pressure [cm] of soil (Brooks-Corey)
@@ -55,6 +59,10 @@
     throughfall::Vector{T}
     # A scaling parameter [mm⁻¹] (controls exponential decline of kv₀)
     f::Vector{T} | "mm-1"
+    # Depth [mm] from soil surface for which exponential decline of kv₀ is valid
+    z_exp::Vector{T} | "mm"
+    # Depth [mm] from soil surface for which layered profile is valid
+    z_layered::Vector{T} | "mm"
     # Amount of water in the unsaturated store, per layer [mm]
     ustorelayerdepth::Vector{SVector{N,T}} | "mm"
     # Saturated store [mm]
@@ -254,6 +262,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
 
     Δt = Second(config.timestepsecs)
     config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
+    ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
     if length(config_thicknesslayers) > 0
         thicknesslayers = SVector(Tuple(push!(Float.(config_thicknesslayers), mv)))
         sumlayers = pushfirst(cumsum(thicknesslayers), 0.0)
@@ -463,9 +472,15 @@ function initialize_sbm(nc, config, riverfrac, inds)
 
     cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
 
+    θₑ = θₛ .- θᵣ
+    soilwatercapacity = soilthickness .* θₑ
+    satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth
+    zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi
+
     # these are filled in the loop below
     # TODO see if we can replace this approach
     nlayers = zeros(Int, n)
+    n_unsatlayers = zeros(Int, n)
     act_thickl = zeros(Float, maxlayers, n)
     s_layers = zeros(Float, maxlayers + 1, n)
 
@@ -478,39 +493,93 @@ function initialize_sbm(nc, config, riverfrac, inds)
             nlayers[i] = nlayers_
             act_thickl[:, i] = act_thickl_
             s_layers[:, i] = s_layers_
+            _, n_unsatlayers[i] = set_layerthickness(zi[i], sumlayers, thicknesslayers)
         else
             nlayers[i] = 1
             act_thickl[:, i] = SVector(soilthickness[i])
             s_layers[:, i] = pushfirst(cumsum(SVector(soilthickness[i])), 0.0)
         end
     end
-
-    # needed for derived parameters below
     act_thickl = svectorscopy(act_thickl, Val{maxlayers}())
-    θₑ = θₛ .- θᵣ
-    soilwatercapacity = soilthickness .* θₑ
-    satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth
-    zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi
 
     # copied to array of sarray below
     vwc = fill(mv, maxlayers, n)
     vwc_perc = fill(mv, maxlayers, n)
+    sumlayers = svectorscopy(s_layers, Val{maxlayers + 1}())
+
+    # ksat profiles
+    if ksat_profile == "exponential"
+        z_exp = soilthickness
+        z_layered = fill(mv, n)
+        kv = fill(mv, (maxlayers, n))
+        nlayers_kv = fill(0, n)
+    elseif ksat_profile == "exponential_constant"
+        z_exp =
+            ncread(nc, config, "vertical.z_exp"; optional = false, sel = inds, type = Float)
+        z_layered = fill(mv, n)
+        kv = fill(mv, (maxlayers, n))
+        nlayers_kv = fill(0, n)
+    elseif ksat_profile == "layered" || ksat_profile == "layered_exponential"
+        z_exp = fill(mv, n)
+        kv =
+            ncread(
+                nc,
+                config,
+                "vertical.kv";
+                sel = inds,
+                defaults = 1000.0,
+                type = Float,
+                dimname = :layer,
+            ) .* (Δt / basetimestep)
+        if size(kv, 1) != maxlayers
+            parname = param(config.input.vertical, "kv")
+            size1 = size(kv, 1)
+            error("$parname needs a layer dimension of size $maxlayers, but is $size1")
+        end
+        if ksat_profile == "layered"
+            z_layered = soilthickness
+            nlayers_kv = nlayers
+        else
+            z_layered = ncread(
+                nc,
+                config,
+                "vertical.z_layered";
+                optional = false,
+                sel = inds,
+                type = Float,
+            )
+            nlayers_kv = fill(0, n)
+            for i in eachindex(nlayers_kv)
+                layers = @view sumlayers[i][2:nlayers[i]]
+                _, k = findmin(abs.(z_layered[i] .- layers))
+                nlayers_kv[i] = k
+                z_layered[i] = layers[k]
+            end
+        end
+    else
+        error("""An unknown "ksat_profile" is specified in the TOML file ($ksat_profile).
+              This should be "exponential", "exponential_constant", "layered" or
+              "layered_exponential".
+              """)
+    end
 
     sbm = SBM{Float,maxlayers,maxlayers + 1}(
         Δt = tosecond(Δt),
         maxlayers = maxlayers,
         n = n,
         nlayers = nlayers,
-        n_unsatlayers = fill(0, n),
+        n_unsatlayers = n_unsatlayers,
+        nlayers_kv = nlayers_kv,
         riverfrac = riverfrac,
         θₛ = θₛ,
         θᵣ = θᵣ,
         kv₀ = kv₀,
+        kv = svectorscopy(kv, Val{maxlayers}()),
         kvfrac = svectorscopy(kvfrac, Val{maxlayers}()),
         hb = hb,
         soilthickness = soilthickness,
         act_thickl = act_thickl,
-        sumlayers = svectorscopy(s_layers, Val{maxlayers + 1}()),
+        sumlayers = sumlayers,
         infiltcappath = infiltcappath,
         infiltcapsoil = infiltcapsoil,
         maxleakage = maxleakage,
@@ -525,6 +594,8 @@ function initialize_sbm(nc, config, riverfrac, inds)
         stemflow = fill(mv, n),
         throughfall = fill(mv, n),
         f = f,
+        z_exp = z_exp,
+        z_layered = z_layered,
         ustorelayerdepth = zero(act_thickl),
         satwaterdepth = satwaterdepth,
         zi = zi,
@@ -692,6 +763,7 @@ function update_until_recharge(sbm::SBM, config)
     modelsnow = get(config.model, "snow", false)::Bool
     transfermethod = get(config.model, "transfermethod", false)::Bool
     ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code?
+    ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
 
     threaded_foreach(1:sbm.n, basesize = 250) do i
         if modelsnow
@@ -775,7 +847,7 @@ function update_until_recharge(sbm::SBM, config)
             # (ast) and the updated unsaturated storage for each soil layer.
             if transfermethod && sbm.maxlayers == 1
                 ustorelayerdepth = sbm.ustorelayerdepth[i][1] + infiltsoilpath
-                kv_z = sbm.kvfrac[i][1] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i])
+                kv_z = hydraulic_conductivity_at_depth(sbm, sbm.zi[i], i, 1, ksat_profile)
                 ustorelayerdepth, ast = unsatzone_flow_sbm(
                     ustorelayerdepth,
                     sbm.soilwatercapacity[i],
@@ -789,7 +861,7 @@ function update_until_recharge(sbm::SBM, config)
             else
                 for m = 1:n_usl
                     l_sat = usl[m] * (sbm.θₛ[i] - sbm.θᵣ[i])
-                    kv_z = sbm.kvfrac[i][m] * sbm.kv₀[i] * exp(-sbm.f[i] * z[m])
+                    kv_z = hydraulic_conductivity_at_depth(sbm, z[m], i, m, ksat_profile)
                     ustorelayerdepth =
                         m == 1 ? sbm.ustorelayerdepth[i][m] + infiltsoilpath :
                         sbm.ustorelayerdepth[i][m] + ast
@@ -901,7 +973,7 @@ function update_until_recharge(sbm::SBM, config)
 
         actcapflux = 0.0
         if n_usl > 0
-            ksat = sbm.kvfrac[i][n_usl] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.zi[i])
+            ksat = hydraulic_conductivity_at_depth(sbm, sbm.zi[i], i, n_usl, ksat_profile)
             ustorecapacity =
                 sbm.soilwatercapacity[i] - satwaterdepth - sum(@view usld[1:sbm.nlayers[i]])
             maxcapflux = max(0.0, min(ksat, actevapustore, ustorecapacity, satwaterdepth))
@@ -925,7 +997,13 @@ function update_until_recharge(sbm::SBM, config)
                 actcapflux = actcapflux + toadd
             end
         end
-        deepksat = sbm.kvfrac[i][end] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.soilthickness[i])
+        deepksat = hydraulic_conductivity_at_depth(
+            sbm,
+            sbm.soilthickness[i],
+            i,
+            sbm.nlayers[i],
+            ksat_profile,
+        )
         deeptransfer = min(satwaterdepth, deepksat)
         actleakage = max(0.0, min(sbm.maxleakage[i], deeptransfer))
 
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 085fca0d4..cdac2694b 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -126,11 +126,15 @@ function initialize_sbm_model(config::Config)
         f = sbm.f .* 1000.0
         zi = sbm.zi .* 0.001
         soilthickness = sbm.soilthickness .* 0.001
+        z_exp = sbm.z_exp .* 0.001
 
         ssf = LateralSSF{Float}(
             kh₀ = kh₀,
             f = f,
+            kh = fill(mv, n),
+            khfrac = khfrac,
             zi = zi,
+            z_exp = z_exp,
             soilthickness = soilthickness,
             θₛ = sbm.θₛ,
             θᵣ = sbm.θᵣ,
@@ -140,11 +144,20 @@ function initialize_sbm_model(config::Config)
             dw = dw,
             exfiltwater = fill(mv, n),
             recharge = fill(mv, n),
-            ssf = ((kh₀ .* βₗ) ./ f) .* (exp.(-f .* zi) - exp.(-f .* soilthickness)) .* dw,
+            ssf = fill(mv, n),
             ssfin = fill(mv, n),
-            ssfmax = ((kh₀ .* βₗ) ./ f) .* (1.0 .- exp.(-f .* soilthickness)),
+            ssfmax = fill(mv, n),
             to_river = zeros(n),
         )
+        # update variables `ssf`, `ssfmax` and `kh` (layered profile) based on ksat_profile
+        ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
+        if ksat_profile == "exponential"
+            initialize_lateralssf_exp!(ssf::LateralSSF)
+        elseif ksat_profile == "exponential_constant"
+            initialize_lateralssf_exp_const!(ssf::LateralSSF)
+        elseif ksat_profile == "layered" || ksat_profile == "layered_exponential"
+            initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
+        end
     else
         # when the SBM model is coupled (BMI) to a groundwater model, the following
         # variables are expected to be exchanged from the groundwater model.
@@ -368,13 +381,21 @@ end
 function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
 
     @unpack lateral, vertical, network, clock, config = model
+    ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
+
     model = update_until_recharge(model)
     # exchange of recharge between vertical sbm concept and subsurface flow domain
     lateral.subsurface.recharge .= vertical.recharge ./ 1000.0
     lateral.subsurface.recharge .*= lateral.subsurface.dw
     lateral.subsurface.zi .= vertical.zi ./ 1000.0
     # update lateral subsurface flow domain (kinematic wave)
-    update(lateral.subsurface, network.land, network.frac_toriver)
+    if (ksat_profile == "layered") || (ksat_profile == "layered_exponential")
+        for i in eachindex(lateral.subsurface.kh)
+            lateral.subsurface.kh[i] =
+                kh_layered_profile(vertical, lateral.subsurface.khfrac[i], i, ksat_profile)
+        end
+    end
+    update(lateral.subsurface, network.land, network.frac_toriver, ksat_profile)
     model = update_after_subsurfaceflow(model)
     model = update_total_water_storage(model)
 end
diff --git a/src/utils.jl b/src/utils.jl
index b229be67d..11ac219c1 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -656,3 +656,173 @@ function threaded_foreach(f, x::AbstractArray; basesize::Integer)
     end
     return nothing
 end
+
+"""
+    hydraulic_conductivity_at_depth(sbm::SBM, z, i, ksat_profile)
+
+Return vertical hydraulic conductivity `kv_z` for soil layer `n` at depth `z` for vertical
+concept `SBM` (at index `i`) based on hydraulic conductivity profile `ksat_profile`.
+"""
+function hydraulic_conductivity_at_depth(sbm::SBM, z, i, n, ksat_profile)
+    if ksat_profile == "exponential"
+        kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z)
+    elseif ksat_profile == "exponential_constant"
+        if z < sbm.z_exp[i]
+            kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z)
+        else
+            kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.z_exp[i])
+        end
+    elseif ksat_profile == "layered"
+        kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n]
+    elseif ksat_profile == "layered_exponential"
+        if z < sbm.z_layered[i]
+            kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n]
+        else
+            n = sbm.nlayers_kv[i]
+            kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n] * exp(-sbm.f[i] * (z - sbm.z_layered[i]))
+        end
+    else
+        error("unknown ksat_profile")
+    end
+    return kv_z
+end
+
+"""
+    kh_layered_profile(sbm::SBM, khfrac, z, i, ksat_profile)
+
+Return equivalent horizontal hydraulic conductivity `kh` [m d⁻¹] for a layered soil profile
+of vertical concept `SBM` (at index `i`) based on multiplication factor `khfrac` [-], water
+table depth `z` [mm] and hydraulic conductivity profile `ksat_profile`.
+"""
+function kh_layered_profile(sbm::SBM, khfrac, i, ksat_profile)
+
+    m = sbm.nlayers[i]
+    t_factor = (tosecond(basetimestep) / sbm.Δt)
+    if (sbm.soilthickness[i] - sbm.zi[i]) > 0.0
+        transmissivity = 0.0
+        sumlayers = @view sbm.sumlayers[i][2:end]
+        n = max(sbm.n_unsatlayers[i], 1)
+
+        if ksat_profile == "layered"
+            transmissivity += (sumlayers[n] - sbm.zi[i]) * sbm.kv[i][n]
+        elseif ksat_profile == "layered_exponential"
+            if sbm.zi[i] >= sbm.z_layered[i]
+                zt = sbm.soilthickness[i] - sbm.z_layered[i]
+                j = sbm.nlayers_kv[i]
+                transmissivity +=
+                    sbm.kv[i][j] / sbm.f[i] *
+                    (exp(-sbm.f[i] * (sbm.zi[i] - sbm.z_layered[i])) - exp(-sbm.f[i] * zt))
+                n = m
+            else
+                transmissivity += (sumlayers[n] - sbm.zi[i]) * sbm.kv[i][n]
+            end
+        else
+            error("unknown ksat_profile, `layered` or `layered_exponential` are allowed")
+        end
+        n += 1
+        while n <= m
+            if ksat_profile == "layered"
+                transmissivity += sbm.act_thickl[i][n] * sbm.kv[i][n]
+            elseif ksat_profile == "layered_exponential"
+                if n > sbm.nlayers_kv[i]
+                    zt = sbm.soilthickness[i] - sbm.z_layered[i]
+                    j = sbm.nlayers_kv[i]
+                    transmissivity += sbm.kv[i][j] / sbm.f[i] * (1.0 - exp(-sbm.f[i] * zt))
+                    n = m
+                else
+                    transmissivity += sbm.act_thickl[i][n] * sbm.kv[i][n]
+                end
+            end
+            n += 1
+        end
+        # convert units for kh [m d⁻¹] computation (transmissivity [mm² Δt⁻¹], soilthickness
+        # [mm] and zi [mm])
+        kh =
+            0.001 *
+            (transmissivity / (sbm.soilthickness[i] - sbm.zi[i])) *
+            t_factor *
+            khfrac
+    else
+        if ksat_profile == "layered"
+            kh = 0.001 * sbm.kv[i][m] * t_factor * khfrac
+        elseif ksat_profile == "layered_exponential"
+            if sbm.zi[i] >= sbm.z_layered[i]
+                j = sbm.nlayers_kv[i]
+                kh =
+                    0.001 *
+                    sbm.kv[i][j] *
+                    exp(-sbm.f[i] * (sbm.zi[i] - sbm.z_layered[i])) *
+                    khfrac *
+                    t_factor
+            else
+                kh = 0.001 * sbm.kv[i][m] * t_factor * khfrac
+            end
+        else
+            error("unknown ksat_profile, `layered` or `layered_exponential` are allowed")
+        end
+    end
+    return kh
+end
+
+"Initialize lateral subsurface variables `ssf` and `ssfmax` with `ksat_profile` `exponential`"
+function initialize_lateralssf_exp!(ssf::LateralSSF)
+    for i in eachindex(ssf.ssf)
+        ssf.ssfmax[i] =
+            ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) *
+            (1.0 - exp(-ssf.f[i] * ssf.soilthickness[i]))
+        ssf.ssf[i] =
+            ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) *
+            (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.soilthickness[i])) *
+            ssf.dw[i]
+    end
+end
+
+"Initialize lateral subsurface variables `ssf` and `ssfmax` with `ksat_profile` `exponential_constant`"
+function initialize_lateralssf_exp_const!(ssf::LateralSSF)
+    ssf_constant = @. ssf.khfrac *
+       ssf.kh₀ *
+       exp(-ssf.f * ssf.z_exp) *
+       ssf.βₗ *
+       (ssf.soilthickness - ssf.z_exp)
+    for i in eachindex(ssf.ssf)
+        ssf.ssfmax[i] =
+            ((ssf.khfrac[i] * ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) *
+            (1.0 - exp(-ssf.f[i] * ssf.z_exp[i])) + ssf_constant[i]
+        if ssf.zi[i] < ssf.z_exp[i]
+            ssf.ssf[i] =
+                (
+                    ((ssf.kh₀[i] * ssf.βₗ[i]) / ssf.f[i]) *
+                    (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.z_exp[i])) +
+                    ssf_constant[i]
+                ) * ssf.dw[i]
+        else
+            ssf.ssf[i] =
+                ssf.kh₀[i] *
+                exp(-ssf.f[i] * ssf.zi[i]) *
+                ssf.βₗ[i] *
+                (ssf.soilthickness[i] - ssf.zi[i]) *
+                ssf.dw[i]
+        end
+    end
+end
+
+"Initialize lateral subsurface variables `ssf`, `ssfmax` and `kh` with ksat_profile` `layered` or `layered_exponential`"
+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.βₗ[i] * ssf.dw[i]
+        kh_max = 0.0
+        for j = 1:sbm.nlayers[i]
+            if j <= sbm.nlayers_kv[i]
+                kh_max += sbm.kv[i][j] * sbm.act_thickl[i][j]
+            else
+                zt = sbm.soilthickness[i] - sbm.z_layered[i]
+                k = max(j - 1, 1)
+                kh_max += sbm.kv[i][k] / sbm.f[i] * (1.0 - exp(-sbm.f[i] * zt))
+                break
+            end
+        end
+        kh_max = kh_max * ssf.khfrac[i] * 0.001 * 0.001
+        ssf.ssfmax[i] = kh_max * ssf.βₗ[i]
+    end
+end
diff --git a/test/bmi.jl b/test/bmi.jl
index 848393681..e972574fb 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -20,9 +20,9 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "model information functions" begin
             @test BMI.get_component_name(model) == "sbm"
-            @test BMI.get_input_item_count(model) == 181
-            @test BMI.get_output_item_count(model) == 181
-            @test BMI.get_input_var_names(model)[[1, 5, 151, 175]] == [
+            @test BMI.get_input_item_count(model) == 192
+            @test BMI.get_output_item_count(model) == 192
+            @test BMI.get_input_var_names(model)[[1, 6, 161, 186]] == [
                 "vertical.nlayers",
                 "vertical.θᵣ",
                 "lateral.river.q",
diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl
index ea7873af2..851429dd7 100644
--- a/test/horizontal_process.jl
+++ b/test/horizontal_process.jl
@@ -62,6 +62,8 @@ end
                 1697.05 * 1000.0,
                 1200.0 * 1000.0,
                 2443723.716252628,
+                1.0,
+                "exponential",
             ),
             (7.410313985168225e10, 1540.1496836278836, -0.0),
         ),
@@ -230,7 +232,7 @@ end
         sw_river.inwater[1] = 20.0
         h0 = mean(sw_river.h)
         Δt = Wflow.stable_timestep(sw_river)
-        Wflow.shallowwater_river_update(sw_river, network, Δt, nothing, 0.0, true)
+        Wflow.shallowwater_river_update(sw_river, network, Δt, 0.0, true)
         d = abs(h0 - mean(sw_river.h))
         if d <= ϵ
             break
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 5f22f35a0..7e2194976 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -404,3 +404,109 @@ model = Wflow.run_timestep(model)
     @test h[501] ≈ 0.05610231297517167f0
 end
 Wflow.close_files(model, delete_output = false)
+
+# test different ksat profiles
+@testset "ksat profiles (SBM)" begin
+    i = 100
+    tomlpath = joinpath(@__DIR__, "sbm_config.toml")
+    config = Wflow.Config(tomlpath)
+    config.input.vertical.kv = "kv"
+    config.input.vertical.z_exp = Dict("value" => 400.0)
+    config.input.vertical.z_layered = Dict("value" => 400.0)
+
+    @testset "exponential profile" begin
+        model = Wflow.initialize_sbm_model(config)
+        @unpack vertical = model
+        z = vertical.zi[i]
+        kv_z = Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential")
+        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z)
+        @test vertical.z_exp == vertical.soilthickness
+        @test_throws ErrorException Wflow.kh_layered_profile(
+            vertical,
+            100.0,
+            i,
+            "exponential",
+        )
+        @test all(isnan.(vertical.z_layered))
+        @test all(isnan.(vertical.kv[i]))
+        @test all(vertical.nlayers_kv .== 0)
+    end
+
+    @testset "exponential constant profile" begin
+        config.input.vertical.ksat_profile = "exponential_constant"
+        model = Wflow.initialize_sbm_model(config)
+        @unpack vertical = model
+        z = vertical.zi[i]
+        kv_z =
+            Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential_constant")
+        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z)
+        kv_400 = Wflow.hydraulic_conductivity_at_depth(
+            vertical,
+            400.0,
+            i,
+            2,
+            "exponential_constant",
+        )
+        kv_1000 = Wflow.hydraulic_conductivity_at_depth(
+            vertical,
+            1000.0,
+            i,
+            3,
+            "exponential_constant",
+        )
+        @test kv_400 ≈ kv_1000
+        @test_throws ErrorException Wflow.kh_layered_profile(
+            vertical,
+            100.0,
+            i,
+            "exponential_constant",
+        )
+        @test all(isnan.(vertical.z_layered))
+        @test all(isnan.(vertical.kv[i]))
+        @test all(vertical.nlayers_kv .== 0)
+        @test all(vertical.z_exp .== 400.0)
+    end
+
+    @testset "layered profile" begin
+        config.input.vertical.ksat_profile = "layered"
+        model = Wflow.initialize_sbm_model(config)
+        @unpack vertical = model
+        z = vertical.zi[i]
+        @test Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "layered") ≈
+              vertical.kv[100][2]
+        @test Wflow.kh_layered_profile(vertical, 100.0, i, "layered") ≈ 47.508932674632355f0
+        @test vertical.nlayers_kv[i] == 4
+        @test vertical.z_layered == vertical.soilthickness
+        @test all(isnan.(vertical.z_exp))
+    end
+
+    @testset "layered exponential profile" begin
+        config.input.vertical.ksat_profile = "layered_exponential"
+        model = Wflow.initialize_sbm_model(config)
+        @unpack vertical = model
+        z = vertical.zi[i]
+        @test Wflow.hydraulic_conductivity_at_depth(
+            vertical,
+            z,
+            i,
+            2,
+            "layered_exponential",
+        ) ≈ vertical.kv[i][2]
+        @test vertical.nlayers_kv[i] == 2
+        @test Wflow.kh_layered_profile(vertical, 100.0, i, "layered_exponential") ≈
+              33.76026208801769f0
+        @test all(vertical.z_layered[1:10] .== 400.0)
+        @test all(isnan.(vertical.z_exp))
+    end
+
+    model = Wflow.run_timestep(model)
+    model = Wflow.run_timestep(model)
+    @testset "river flow layered exponential profile" begin
+        q = model.lateral.river.q_av
+        @test sum(q) ≈ 3118.8690178033266f0
+        @test q[1622] ≈ 0.000548447582354063f0
+        @test q[43] ≈ 9.860543811678328f0
+    end
+
+    Wflow.close_files(model, delete_output = false)
+end