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

Issue:
Support different vertical hydraulic conductivity (`kv`) profiles in `SBM`
### Feature type

Adding new functionality

### Improvement Description

The vertical hydraulic conductivity at the surface `kv₀` declines with depth exponentially, controlled by the model parameter `f`. This vertical hydraulic conductivity profile is often based on data available for soil layers to a depth of ~2 m. For below 2 m an exponential decline of `kv₀` may not always be the best choice, and a user may want to use available vertical hydraulic conductivity data of soil layers directly in the model. Because of this, an improvement is to support the following `kv` profiles:
- `exponential`: Exponential decline (default).
- `exponential_constant`: Exponential decline (for example until soil layers with measurements (e.g. 2 m)) and constant (rest of the soil column).
- `layered`: `kv` value for each soil layer.
- `layered_exponential`: Layered (for example until soil layers with measurements (e.g. 2 m))  and exponential (rest of the soil column).

### Implementation Description

Add different profiles to initialization of wflow_sbm, vertical unsaturated flow of `SBM` and lateral subsurface flow (kinematic wave).
For `layered` and `layered_exponential` profiles, an average horizontal hydraulic conductivity can used for lateral subsurface flow.

### Additional Context

_No response_
Support different vertical hydraulic conductivity (`kv`) profiles in `SBM`
### Feature type

Adding new functionality

### Improvement Description

The vertical hydraulic conductivity at the surface `kv₀` declines with depth exponentially, controlled by the model parameter `f`. This vertical hydraulic conductivity profile is often based on data available for soil layers to a depth of ~2 m. For below 2 m an exponential decline of `kv₀` may not always be the best choice, and a user may want to use available vertical hydraulic conductivity data of soil layers directly in the model. Because of this, an improvement is to support the following `kv` profiles:
- `exponential`: Exponential decline (default).
- `exponential_constant`: Exponential decline (for example until soil layers with measurements (e.g. 2 m)) and constant (rest of the soil column).
- `layered`: `kv` value for each soil layer.
- `layered_exponential`: Layered (for example until soil layers with measurements (e.g. 2 m))  and exponential (rest of the soil column).

### Implementation Description

Add different profiles to initialization of wflow_sbm, vertical unsaturated flow of `SBM` and lateral subsurface flow (kinematic wave).
For `layered` and `layered_exponential` profiles, an average horizontal hydraulic conductivity can used for lateral subsurface flow.

### Additional Context

_No response_


Discussion / hints:

Fix patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 0b1f22c54..318667d20 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -32,6 +32,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/src/flow.jl b/src/flow.jl
index c0c0a1515..220784593 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -391,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 [-]
@@ -399,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⁻¹]
@@ -413,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)
@@ -435,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
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


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index 848393681..af6554220 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) == 191
+            @test BMI.get_output_item_count(model) == 191
+            @test BMI.get_input_var_names(model)[[1, 6, 161, 185]] == [
                 "vertical.nlayers",
                 "vertical.θᵣ",
                 "lateral.river.q",
diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl
index ea7873af2..9e395a098 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),
         ),
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