Repository: Deltares/Wflow.jl
Pull request #114 (2021-03-19 07:22:04)

Check params
Review and changes of struct parameters to get a bit more consistency. An example: changing the Manning's n parameter when the model was already initialized did not result in any changes, because the update of the alpha parameter was only part of the initialization and not the update of the model. Changing parameters after initialization remains tricky, and requires detailed knowledge of the model.

visr: Looks mostly good, great to see the reduction in code.
Looking at `altitude`, is it correct that for `sbm_gwf` altitude is now under `input` in the toml, whereas for other models it is under `input.vertical`? What is the reason for this difference? This change here:

https://github.com/Deltares/Wflow.jl/pull/114/commits/4028bd18357edf0109f09ad16a1481858dc7f97d#diff-5091fa59158833690f71ac46441bdc5de6c5a40b1fd099edebbe8ab21b8054dd

is technically a breaking change to the config file input. If we want that, we should think a bit more about how to approach breaking changes.

Patch:
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 35370c3f1..a70173c56 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -14,6 +14,7 @@ jobs:
       matrix:
         version:
           - '1.5'
+          - '1'
         os:
           - ubuntu-latest
           - macOS-latest
diff --git a/Project.toml b/Project.toml
index 3ac0c4be3..ab895503d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Wflow"
 uuid = "d48b7d99-76e7-47ae-b1d5-ff0c1cf9a818"
 authors = ["Deltares and contributors"]
-version = "0.1.2"
+version = "0.2.0"
 
 [deps]
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
diff --git a/docs/src/lateral/kinwave.md b/docs/src/lateral/kinwave.md
index e13266958..1cc9b1a59 100644
--- a/docs/src/lateral/kinwave.md
+++ b/docs/src/lateral/kinwave.md
@@ -43,13 +43,8 @@ unit width of slope ``w`` [mm] according to:
 where ``\beta`` is element slope angle [deg.], ``q`` is subsurface flow [mm``^{2}``/t],
 ``K_{0}`` is the saturated hydraulic conductivity at the soil surface [mm/t], ``z_{i}`` is
 the water table depth [mm],``z_{t}`` is total soil depth [mm], and ``f`` is a scaling
-parameter [mm``^{-1}``]:
-```math
-    f=\frac{\theta_{s}-\theta_{r}}{M},\,
-```
-where ``\theta_{s}`` is saturated water content [mm/mm] and ``\theta_{r}`` is residual water
-content [mm/mm] and ``M`` represents a model parameter [mm], that determines the decrease of
-vertical saturated conductivity with depth.
+parameter [mm``^{-1}``], that controls the decrease of vertical saturated conductivity with 
+depth.
 
 Combining with the following continuity equation:
 ```math
diff --git a/docs/src/vertical/sbm.md b/docs/src/vertical/sbm.md
index df4865758..3852ca345 100644
--- a/docs/src/vertical/sbm.md
+++ b/docs/src/vertical/sbm.md
@@ -125,13 +125,7 @@ to:
     K_{sat}=K_{0}e^{(-fz)}
 ```
 where ``K_{0}`` is the saturated conductivity at the soil surface and ``f`` is a scaling
-parameter [mm``^{-1}``]. The scaling parameter :math:`f` is defined by:
-
-```math
-f=\frac{\theta_{s}-\theta_{r}}{M} 
-```
-with ``\theta_{s}`` and ``\theta_{r}`` as defined previously and ``M`` [mm] representing a
-model parameter.
+parameter [mm``^{-1}``].
 
 ## Transpiration and soil evaporation
 The potential evaporation left over after interception and open water evaporation (rivers
diff --git a/src/flow.jl b/src/flow.jl
index e1a206d92..2cb312ca4 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -12,7 +12,6 @@
     Δt::T | "s"                             # Model time step [s]
     its::Int | "-"                          # Number of fixed iterations
     width::Vector{T} | "m"                  # Flow width [m]
-    alpha_term::Vector{T} | "s3/5 m-1/5"    # Constant part of α
     alpha_pow::T | "-"                      # Used in the power part of α
     α::Vector{T} | "s3/5 m1/5"              # Constant in momentum equation A = αQᵝ, based on Manning's equation
     eps::T | "s3/5 m1/5"                    # Maximum allowed change in α, if exceeded cross sectional area and h is recalculated
@@ -69,6 +68,9 @@ function update(
     # sub time step
     adt = sf.Δt / its
 
+    alpha_term = pow.(sf.n ./ sqrt.(sf.sl),sf.β)
+    sf.α .= alpha_term .* pow.(sf.width .+ 2.0 .* sf.h, sf.alpha_pow)
+
     q_sum = zeros(n)
     h_sum = zeros(n)
     sf.to_river .= 0.0
@@ -135,7 +137,7 @@ function update(
             sf.h[v] = crossarea / sf.width[v]
             wetper = sf.width[v] + (2.0 * sf.h[v]) # wetted perimeter
             α = sf.α[v]
-            sf.α[v] = sf.alpha_term[v] * pow(wetper, sf.alpha_pow)
+            sf.α[v] = alpha_term[v] * pow(wetper, sf.alpha_pow)
 
             if abs(α - sf.α[v]) > sf.eps
                 crossarea = sf.α[v] * pow(sf.q[v], sf.β)
@@ -156,7 +158,8 @@ end
     kh₀::Vector{T} | "mm Δt-1"              # Horizontal hydraulic conductivity at soil surface [mm Δt⁻¹]
     f::Vector{T} | "mm-1"                   # A scaling parameter [mm⁻¹] (controls exponential decline of kh₀)
     soilthickness::Vector{T} | "mm"         # Soil thickness [mm]
-    θₑ::Vector{T} | "mm mm-1"               # Effective porosity [-]
+    θₛ::Vector{T} | "mm mm-1"               # Saturated water content (porosity) [-]
+    θᵣ::Vector{T} | "mm mm-1"               # Residual water content [-]
     t::T | "Δt s"                           # time step [Δt s]
     Δt::T | "s"                             # model time step [s]
     βₗ::Vector{T} | "m m-1"                  # Slope [m m⁻¹]
@@ -205,7 +208,7 @@ function update(ssf::LateralSSF, network, frac_toriver)
             ssf.recharge[v],
             ssf.kh₀[v],
             ssf.βₗ[v],
-            ssf.θₑ[v],
+            ssf.θₛ[v] - ssf.θᵣ[v],
             ssf.f[v],
             ssf.soilthickness[v],
             ssf.t,
diff --git a/src/hbv.jl b/src/hbv.jl
index bb0beac87..3843fccd5 100644
--- a/src/hbv.jl
+++ b/src/hbv.jl
@@ -1,9 +1,6 @@
 @get_units @with_kw struct HBV{T}
     Δt::T | "s"                     # Model time step [s]
     n::Int | "-"                    # Number of cells
-    yl::Vector{T} | "m"             # Length of cells in y direction [m]
-    xl::Vector{T} | "m"             # Length of cells in x direction [m]
-    altitude::Vector{T} | "m"       # Vertical elevation [m]
     fc::Vector{T}                   # Field capacity [mm]
     betaseepage::Vector{T} | "-"    # Exponent in soil runoff generation equation [-]
     lp::Vector{T} | "-"             # Fraction of field capacity below which actual evaporation=potential evaporation [-]
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index 71579f923..f0ffa3575 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -288,9 +288,6 @@ function initialize_hbv_model(config::Config)
     hbv = HBV{Float64}(
         Δt = tosecond(Δt),
         n = n,
-        yl = yl,
-        xl = xl,
-        altitude = altitude,
         fc = fc,
         betaseepage = betaseepage,
         lp = lp,
@@ -408,8 +405,6 @@ function initialize_hbv_model(config::Config)
 
     alpha_pow = (2.0 / 3.0) * 0.6
     β = 0.6
-    h = fill(0.0, n)
-    alpha_term = pow.(n_land ./ sqrt.(βₗ), β)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
@@ -419,15 +414,14 @@ function initialize_hbv_model(config::Config)
         qin = fill(0.0, n),
         q_av = fill(0.0, n),
         qlat = fill(0.0, n),
-        h = h,
+        h = fill(0.0, n),
         h_av = fill(0.0, n),
         Δt = tosecond(Δt),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = dw,
         wb_pit = pits[inds],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(dw .+ 2.0 .* h, alpha_pow),
+        α = fill(mv, n),
         eps = 1e-03,
         cel = fill(0.0, n),
         to_river = fill(0.0, n),
@@ -467,8 +461,6 @@ function initialize_hbv_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)
 
-    h_river = fill(0.0, nriv)
-    alpha_term = pow.(n_river ./ sqrt.(riverslope), β)
     rf = SurfaceFlow(
         β = β,
         sl = riverslope,
@@ -478,16 +470,15 @@ function initialize_hbv_model(config::Config)
         qin = fill(0.0, nriv),
         q_av = fill(0.0, nriv),
         qlat = fill(0.0, nriv),
-        h = h_river,
+        h = fill(0.0, nriv),
         h_av = fill(0.0, nriv),
         Δt = tosecond(Δt),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(riverwidth .+ 2.0 .* h_river, alpha_pow),
+        α = fill(mv, nriv),
         eps = 1e-03,
         cel = fill(0.0, nriv),
         to_river = fill(0.0, nriv),
@@ -517,6 +508,8 @@ function initialize_hbv_model(config::Config)
         order = topological_sort_by_dfs(graph),
         indices = inds,
         reverse_indices = rev_inds,
+        xl = xl,
+        yl = yl,
     )
     river = (
         graph = graph_riv,
@@ -583,7 +576,7 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:HBV,R,W}
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # hbv concept
     lateral.land.qlat .=
-        (vertical.runoff .* vertical.xl .* vertical.yl .* 0.001) ./ lateral.land.Δt ./
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt ./
         lateral.land.dl
 
     # run kinematic wave for overland flow
diff --git a/src/sbm.jl b/src/sbm.jl
index ca4c8fdbb..4f3bf6d9e 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -9,24 +9,16 @@
     nlayers::Vector{Int} | "-"
     # Number of unsaturated soil layers
     n_unsatlayers::Vector{Int} | "-"
-    # length of cells in y direction [m]
-    yl::Vector{T} | "m"
-    # length of cells in x direction [m]
-    xl::Vector{T} | "m"
     # Fraction of river [-]
     riverfrac::Vector{T} | "-"
     # Saturated water content (porosity) [mm mm⁻¹]
     θₛ::Vector{T} | "mm mm-1"
     # Residual water content [mm mm⁻¹]
     θᵣ::Vector{T} | "mm mm-1"
-    # Effictive porosity [mm mm⁻¹]
-    θₑ::Vector{T} | "mm mm-1"
     # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface
     kv₀::Vector{T} | "mm Δt-1"
     # Muliplication factor [-] applied to kv_z (vertical flow)
     kvfrac::Vector{SVector{N,T}} | "-"
-    # Parameter [mm] controlling f
-    m::Vector{T}
     # Air entry pressure [cm] of soil (Brooks-Corey)
     hb::Vector{T} | "cm"
     # Soil thickness [mm]
@@ -45,8 +37,6 @@
     waterfrac::Vector{T} | "-"
     # Fraction of compacted area  [-]
     pathfrac::Vector{T} | "-"
-    # Vertical elevation [m]
-    altitude::Vector{T} | "m"
     # Rooting depth [mm]
     rootingdepth::Vector{T}
     # Controls how roots are linked to water table [-]
@@ -92,8 +82,6 @@
     transpiration::Vector{T}
     # Actual evaporation from unsaturated store [mm]
     ae_ustore::Vector{T}
-    # Actual evaporation from saturated store [mm]
-    ae_sat::Vector{T}
     # Interception [mm]
     interception::Vector{T}
     # Soil evaporation from unsaturated store [mm]
@@ -102,7 +90,7 @@
     soilevapsat::Vector{T}
     # Actual capillary rise [mm]
     actcapflux::Vector{T}
-    # Actual evaporation from saturated store (transpiration and soil evaporation) [mm]
+    # Actual transpiration from saturated store [mm]
     actevapsat::Vector{T}
     # Total actual evapotranspiration [mm]
     actevap::Vector{T}
@@ -266,7 +254,7 @@ function initialize_canopy(nc, config, inds)
     return cmax, e_r, canopygapfraction, sl, swood, kext
 end
 
-function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
+function initialize_sbm(nc, config, riverfrac, inds)
 
     Δt = Second(config.timestepsecs)
     config_thicknesslayers = get(config.model, "thicknesslayers", Float64[])
@@ -280,9 +268,6 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
 
     n = length(inds)
 
-    altitude =
-        ncread(nc, param(config, "input.vertical.altitude"); sel = inds, type = Float64)
-
     cfmax =
         ncread(
             nc,
@@ -399,11 +384,11 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
             defaults = 3000.0,
             type = Float64,
         ) .* (Δt / basetimestep)
-    m = ncread(
+    f = ncread(
         nc,
-        param(config, "input.vertical.m", nothing);
+        param(config, "input.vertical.f", nothing);
         sel = inds,
-        defaults = 300.0,
+        defaults = 0.001,
         type = Float64,
     )
     hb = ncread(
@@ -559,15 +544,11 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
         n = n,
         nlayers = nlayers,
         n_unsatlayers = fill(0, n),
-        yl = yl,
-        xl = xl,
         riverfrac = riverfrac,
         θₛ = θₛ,
         θᵣ = θᵣ,
-        θₑ = θₑ,
         kv₀ = kv₀,
         kvfrac = svectorscopy(kvfrac, Val{maxlayers}()),
-        m = m,
         hb = hb,
         soilthickness = soilthickness,
         act_thickl = act_thickl,
@@ -577,7 +558,6 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
         maxleakage = maxleakage,
         waterfrac = max.(waterfrac .- riverfrac, 0.0),
         pathfrac = pathfrac,
-        altitude = altitude,
         rootingdepth = rootingdepth,
         rootdistpar = rootdistpar,
         capscale = capscale,
@@ -585,7 +565,7 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
         c = svectorscopy(c, Val{maxlayers}()),
         stemflow = fill(mv, n),
         throughfall = fill(mv, n),
-        f = θₑ ./ m,
+        f = f,
         ustorelayerdepth = zero(act_thickl),
         satwaterdepth = satwaterdepth,
         zi = zi,
@@ -600,7 +580,6 @@ function initialize_sbm(nc, config, riverfrac, xl, yl, inds)
         pottrans_soil = fill(mv, n),
         transpiration = fill(mv, n),
         ae_ustore = fill(mv, n),
-        ae_sat = fill(mv, n),
         interception = fill(mv, n),
         soilevap = fill(mv, n),
         soilevapsat = fill(mv, n),
@@ -670,7 +649,6 @@ function update_until_snow(sbm::SBM, config)
     do_lai = haskey(config.input.vertical, "leaf_area_index")
     modelglacier = get(config.model, "glacier", false)::Bool
     modelsnow = get(config.model, "snow", false)::Bool
-    Δt = Second(config.timestepsecs)
 
     for i = 1:sbm.n
         if do_lai
@@ -692,7 +670,7 @@ function update_until_snow(sbm::SBM, config)
         # should we include tempcor in SBM?
         # potential_evaporation = PotenEvap #??
 
-        if Δt >= Hour(23)
+        if Second(sbm.Δt) >= Hour(23)
             throughfall, interception, stemflow, canopystorage = rainfall_interception_gash(
                 cmax,
                 e_r,
@@ -756,7 +734,6 @@ 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?
-    Δt = Second(config.timestepsecs)
 
     for i = 1:sbm.n
         if modelsnow
@@ -774,7 +751,7 @@ function update_until_recharge(sbm::SBM, config)
                     sbm.g_tt[i],
                     sbm.g_cfmax[i],
                     sbm.g_sifrac[i],
-                    Δt,
+                    Second(sbm.Δt),
                 )
                 # Convert to mm per grid cell and add to snowmelt
                 glaciermelt = glaciermelt * sbm.glacierfrac[i]
@@ -975,7 +952,7 @@ function update_until_recharge(sbm::SBM, config)
             if sbm.zi[i] > rootingdepth
                 capfluxscale =
                     sbm.capscale[i] / (sbm.capscale[i] + sbm.zi[i] - rootingdepth) *
-                    tosecond(Δt) / tosecond(basetimestep)
+                    (sbm.Δt / tosecond(basetimestep))
             else
                 capfluxscale = 0.0
             end
@@ -1017,6 +994,7 @@ function update_until_recharge(sbm::SBM, config)
         sbm.runoff_river[i] = runoff_river
         sbm.actevapsat[i] = actevapsat
         sbm.actevap[i] = actevap
+        sbm.ae_ustore[i] = actevapustore
         sbm.ustorelayerdepth[i] = usld
         sbm.transfer[i] = transfer
         sbm.actcapflux[i] = actcapflux
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index da3874ca0..6859c5e0c 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -59,6 +59,8 @@ function initialize_sbm_gwf_model(config::Config)
         ncread(nc, param(config, "input.lateral.river.length"); type = Float64, fill = 0)
     riverlength = riverlength_2d[inds]
 
+    altitude =
+        ncread(nc, param(config, "input.altitude"); sel = inds, type = Float64)
     # read x, y coordinates and calculate cell length [m]
     y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
     x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
@@ -81,7 +83,7 @@ function initialize_sbm_gwf_model(config::Config)
     end
 
     # initialize vertical SBM concept
-    sbm = initialize_sbm(nc, config, riverfrac, xl, yl, inds)
+    sbm = initialize_sbm(nc, config, riverfrac, inds)
 
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
     nriv = length(inds_riv)
@@ -128,8 +130,6 @@ function initialize_sbm_gwf_model(config::Config)
 
     alpha_pow = (2.0 / 3.0) * 0.6
     β = 0.6
-    h = fill(0.0, n)
-    alpha_term = pow.(n_land ./ sqrt.(βₗ), β)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
@@ -139,15 +139,14 @@ function initialize_sbm_gwf_model(config::Config)
         qin = fill(0.0, n),
         q_av = fill(0.0, n),
         qlat = fill(0.0, n),
-        h = h,
+        h = fill(0.0, n),
         h_av = fill(0.0, n),
         Δt = tosecond(Δt),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = sw,
         wb_pit = pits[inds],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(sw .+ 2.0 .* h, alpha_pow),
+        α = fill(mv, n),
         eps = 1e-03,
         cel = fill(0.0, n),
         to_river = fill(0.0, n),
@@ -185,8 +184,6 @@ function initialize_sbm_gwf_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)
 
-    h_river = fill(0.0, nriv)
-    alpha_term = pow.(n_river ./ sqrt.(riverslope), β)
     rf = SurfaceFlow(
         β = β,
         sl = riverslope,
@@ -196,16 +193,15 @@ function initialize_sbm_gwf_model(config::Config)
         q_av = fill(0.0, nriv),
         qin = fill(0.0, nriv),
         qlat = fill(0.0, nriv),
-        h = h_river,
+        h = fill(0.0, nriv),
         h_av = fill(0.0, nriv),
         Δt = tosecond(Δt),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(riverwidth .+ 2.0 .* h_river, alpha_pow),
+        α = fill(mv, nriv),
         eps = 1e-03,
         cel = fill(0.0, nriv),
         to_river = fill(0.0, nriv),
@@ -242,9 +238,9 @@ function initialize_sbm_gwf_model(config::Config)
         type = Float64,
     )
 
-    connectivity = Connectivity(inds, rev_inds, sbm.xl, sbm.yl)
-    initial_head = sbm.altitude .- 0.10 # cold state for groundwater head
-    initial_head[index_river] = sbm.altitude[index_river]
+    connectivity = Connectivity(inds, rev_inds, xl, yl)
+    initial_head = altitude .- 0.10 # cold state for groundwater head
+    initial_head[index_river] = altitude[index_river]
 
     if do_constanthead
         initial_head[constant_head.index] = constant_head.head
@@ -253,16 +249,16 @@ function initialize_sbm_gwf_model(config::Config)
     aquifer = UnconfinedAquifer(
         initial_head,
         conductivity,
-        sbm.altitude,
-        sbm.altitude .- sbm.soilthickness ./ 1000.0,
+        altitude,
+        altitude .- sbm.soilthickness ./ 1000.0,
         xl .* yl,
         specific_yield,
         fill(0.0, connectivity.nconnection),  # conductance
     )
 
     # reset zi and satwaterdepth with groundwater head from unconfined aquifer 
-    sbm.zi .= (sbm.altitude .- initial_head) .* 1000.0
-    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* sbm.θₑ
+    sbm.zi .= (altitude .- initial_head) .* 1000.0
+    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* (sbm.θₛ .- sbm.θᵣ)
 
     # river boundary of unconfined aquifer
     infiltration_conductance = ncread(
@@ -385,6 +381,9 @@ function initialize_sbm_gwf_model(config::Config)
         order = topological_sort_by_dfs(graph),
         indices = inds,
         reverse_indices = rev_inds,
+        xl = xl,
+        yl = yl,
+        altitude = altitude, 
     )
     river = (
         graph = graph_riv,
@@ -484,14 +483,14 @@ function update_sbm_gwf(model)
     # update vertical sbm concept (runoff, ustorelayerdepth and satwaterdepth)
     update_after_lateralflow(
         vertical,
-        (vertical.altitude .- lateral.subsurface.flow.aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
+        (network.land.altitude .- lateral.subsurface.flow.aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM
         exfiltwater .* 1000.0,
     )
 
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # sbm concept
     lateral.land.qlat .=
-        (vertical.runoff .* vertical.xl .* vertical.yl .* 0.001) ./ lateral.land.Δt ./
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt ./
         lateral.land.dl
     # run kinematic wave for overland flow
     update(
@@ -504,8 +503,8 @@ function update_sbm_gwf(model)
     # determine net runoff from vertical sbm concept in river cells
     net_runoff_river =
         (
-            vertical.net_runoff_river[inds_riv] .* vertical.xl[inds_riv] .*
-            vertical.yl[inds_riv] .* 0.001
+            vertical.net_runoff_river[inds_riv] .* network.land.xl[inds_riv] .*
+            network.land.yl[inds_riv] .* 0.001
         ) ./ vertical.Δt
 
     # flux from groundwater domain to river (Q to river from drains (optional) and groundwater)
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 814a4b1af..7bdb5daaa 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -70,7 +70,7 @@ function initialize_sbm_model(config::Config)
             river[i] ? min((riverlength[i] * riverwidth[i]) / (xl[i] * yl[i]), 1.0) : 0.0
     end
 
-    sbm = initialize_sbm(nc, config, riverfrac, xl, yl, inds)
+    sbm = initialize_sbm(nc, config, riverfrac, inds)
 
     inds_riv, rev_inds_riv = active_indices(river_2d, 0)
     nriv = length(inds_riv)
@@ -126,7 +126,8 @@ function initialize_sbm_model(config::Config)
         f = sbm.f,
         zi = sbm.zi,
         soilthickness = sbm.soilthickness,
-        θₑ = sbm.θₛ .- sbm.θᵣ,
+        θₛ = sbm.θₛ,
+        θᵣ = sbm.θᵣ,
         Δt = tosecond(Δt),
         t = 1.0,
         βₗ = βₗ,
@@ -153,8 +154,6 @@ function initialize_sbm_model(config::Config)
 
     alpha_pow = (2.0 / 3.0) * 0.6
     β = 0.6
-    h = fill(0.0, n)
-    alpha_term = pow.(n_land ./ sqrt.(βₗ), β)
     olf = SurfaceFlow(
         β = β,
         sl = βₗ,
@@ -164,15 +163,14 @@ function initialize_sbm_model(config::Config)
         qin = fill(0.0, n),
         q_av = fill(0.0, n),
         qlat = fill(0.0, n),
-        h = h,
+        h = fill(0.0, n),
         h_av = fill(0.0, n),
         Δt = tosecond(Δt),
         its = kw_land_tstep > 0 ? Int(cld(tosecond(Δt), kw_land_tstep)) : kw_land_tstep,
         width = sw,
         wb_pit = pits[inds],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(sw .+ 2.0 .* h, alpha_pow),
+        α = fill(mv, n),
         eps = 1e-03,
         cel = fill(0.0, n),
         to_river = fill(0.0, n),
@@ -213,8 +211,6 @@ function initialize_sbm_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)
 
-    h_river = fill(0.0, nriv)
-    alpha_term = pow.(n_river ./ sqrt.(riverslope), β)
     rf = SurfaceFlow(
         β = β,
         sl = riverslope,
@@ -224,16 +220,15 @@ function initialize_sbm_model(config::Config)
         qin = fill(0.0, nriv),
         q_av = fill(0.0, nriv),
         qlat = fill(0.0, nriv),
-        h = h_river,
+        h = fill(0.0, nriv),
         h_av = fill(0.0, nriv),
         Δt = tosecond(Δt),
         its = kw_river_tstep > 0 ? ceil(Int(tosecond(Δt) / kw_river_tstep)) :
               kw_river_tstep,
         width = riverwidth,
         wb_pit = pits[inds_riv],
-        alpha_term = alpha_term,
         alpha_pow = alpha_pow,
-        α = alpha_term .* pow.(riverwidth .+ 2.0 .* h_river, alpha_pow),
+        α = fill(mv, nriv),
         eps = 1e-03,
         cel = fill(0.0, nriv),
         to_river = fill(0.0, nriv),
@@ -274,6 +269,8 @@ function initialize_sbm_model(config::Config)
         order = topological_sort_by_dfs(graph),
         indices = inds,
         reverse_indices = rev_inds,
+        xl = xl,
+        yl = yl,
     )
     river = (
         graph = graph_riv,
@@ -300,14 +297,8 @@ function initialize_sbm_model(config::Config)
         set_states(instate_path, model, state_ncnames; type = Float64)
         @unpack lateral, vertical = model
         # update zi for vertical sbm and α for river and overland flow
-        zi = max.(0.0, vertical.soilthickness .- vertical.satwaterdepth ./ vertical.θₑ)
+        zi = max.(0.0, vertical.soilthickness .- vertical.satwaterdepth ./ (vertical.θₛ .- vertical.θᵣ))
         vertical.zi .= zi
-        lateral.river.α .=
-            lateral.river.alpha_term .*
-            pow.(lateral.river.width .+ 2.0 .* lateral.river.h, lateral.river.alpha_pow)
-        lateral.land.α .=
-            lateral.land.alpha_term .*
-            pow.(lateral.land.width .+ 2.0 .* lateral.land.h, lateral.land.alpha_pow)
     end
 
     # make sure the forcing is already loaded
@@ -371,7 +362,7 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:SBM,R,W}
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # sbm concept
     lateral.land.qlat .=
-        (vertical.runoff .* vertical.xl .* vertical.yl .* 0.001) ./ lateral.land.Δt ./
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt ./
         lateral.land.dl
     # run kinematic wave for overland flow
     update(
@@ -385,8 +376,8 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:SBM,R,W}
     # overland flow lateral subsurface flow and net runoff to the river cells
     net_runoff_river =
         (
-            vertical.net_runoff_river[inds_riv] .* vertical.xl[inds_riv] .*
-            vertical.yl[inds_riv] .* 0.001
+            vertical.net_runoff_river[inds_riv] .* network.land.xl[inds_riv] .*
+            network.land.yl[inds_riv] .* 0.001
         ) ./ vertical.Δt
     lateral.river.qlat .=
         (


Test patch:
diff --git a/test/bmi.jl b/test/bmi.jl
index d965738ad..e30e73e9b 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -15,13 +15,13 @@ 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) == 182
-        @test BMI.get_output_item_count(model) == 182
-        @test BMI.get_input_var_names(model)[[1, 5, 120, 179]] == [
+        @test BMI.get_input_item_count(model) == 175
+        @test BMI.get_output_item_count(model) == 175
+        @test BMI.get_input_var_names(model)[[1, 5, 120, 175]] == [
             "vertical.Δt",
             "vertical.n_unsatlayers",
-            "lateral.land.sl",
-            "lateral.river.reservoir.percfull",
+            "lateral.land.q_av",
+            "lateral.river.reservoir.evaporation",
         ]
     end
 
@@ -37,7 +37,6 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         @test BMI.get_var_units(model, "vertical.θₛ") == "mm mm-1"
         @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == 8
         @test BMI.get_var_nbytes(model, "vertical.n") == 8
-        @test BMI.get_var_nbytes(model, "vertical.xl") == 400560
         @test BMI.get_var_nbytes(model, "lateral.river.q") == 45240
         @test BMI.get_var_location(model, "lateral.river.q") == "node"
     end
diff --git a/test/run_hbv.jl b/test/run_hbv.jl
index ec1496c55..26735695c 100644
--- a/test/run_hbv.jl
+++ b/test/run_hbv.jl
@@ -28,7 +28,6 @@ end
     hbv = model.vertical
     @test hbv.tt[1] ≈ 0.0
     @test model.clock.iteration == 2
-    @test hbv.altitude[1] == 484.2149963378906
     @test hbv.soilmoisture[1] == 134.35299682617188
     @test hbv.runoff[1] == 7.406898120121746
     @test hbv.soilevap[1] == 0.0
@@ -40,7 +39,6 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     hbv = model.vertical
-    @test hbv.altitude[1] == 484.2149963378906
     @test hbv.soilmoisture[1] == 134.35299682617188
     @test hbv.runoff[1] == 4.353345933656929
     @test hbv.soilevap[1] == 0.0
@@ -57,10 +55,10 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 51941.417623244706
-    @test q[500] ≈ 5.6934140250393055
-    @test q[88] ≈ 8.839325386868788
-    @test q[network.river.order[end]] ≈ 363.528899708456
+    @test sum(q) ≈ 51943.24072462076
+    @test q[500] ≈ 5.693405508134953
+    @test q[88] ≈ 8.839306661545436
+    @test q[network.river.order[end]] ≈ 363.5647479049738
 end
 
 benchmark = @benchmark Wflow.update(model)
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index af843452d..b3beb124a 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -18,25 +18,25 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     @test row.volume ≈ 4.364251626782245e7
     @test row.temp_bycoord ≈ 2.3279826641082764
     @test row.temp_byindex ≈ 2.3279826641082764
-    @test row.Q_6336050 ≈ 0.023884907014593577
-    @test row.Q_6336510 ≈ 0.012411069390755603
+    @test row.Q_6336050 ≈ 0.02388490515704142
+    @test row.Q_6336510 ≈ 0.012411069754330572
     @test row.Q_6836100 ≈ 0.004848609260578462
     @test row.Q_6336500 ≈ 0.011474801875419624
-    @test row.Q_6836190 ≈ 0.0005262528887021589
+    @test row.Q_6836190 ≈ 0.0005262529349807199
     @test row.Q_6336800 ≈ 0.013467698526109831
-    @test row.Q_6336900 ≈ 0.0034082809012198974
+    @test row.Q_6336900 ≈ 0.003408280699908697
     @test row.Q_6336930 ≈ 0.09773275295006544
     @test row.Q_6336910 ≈ 0.002147610203630289
-    @test row.Q_6336920 ≈ 0.002649393554348607
-    @test row.Q_6136100 ≈ 0.0008708128976810045
+    @test row.Q_6336920 ≈ 0.002649393436607859
+    @test row.Q_6136100 ≈ 0.0008708128761381345
     @test row.Q_6136500 ≈ 0.000729148906480041
     @test row.Q_6136520 ≈ 0.002155395279410574
-    @test row.Q_6136150 ≈ 0.0022298330217711665
-    @test row.Q_6136151 ≈ 0.0031045029537777745
+    @test row.Q_6136150 ≈ 0.0022298329229975987
+    @test row.Q_6136151 ≈ 0.0031045027476524524
     @test row.Q_6136160 ≈ 3.3423894540713786
     @test row.Q_6136200 ≈ 1.358270076420503
     @test row.Q_6136201 ≈ 5.942330218662147
-    @test row.Q_6136202 ≈ 1.6809973032577248
+    @test row.Q_6136202 ≈ 1.680997341533662
     @test row.recharge_1 ≈ -0.027398093386017383
 end
 
@@ -77,7 +77,6 @@ end
 
     @test model.clock.iteration == 2
 
-    @test sbm.altitude[1] ≈ 643.5469970703125
     @test sbm.θₛ[1] ≈ 0.48343977332115173
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
@@ -89,7 +88,6 @@ model = Wflow.update(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.altitude[1] ≈ 643.5469970703125
     @test sbm.θₛ[1] ≈ 0.48343977332115173
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] ≈ 0.005865651540305367
@@ -100,8 +98,8 @@ end
     ssf = model.lateral.subsurface.ssf
     @test sum(ssf) ≈ 6.368140761295825e16
     @test ssf[network.land.order[1]] ≈ 3.0449782003445332e13
-    @test ssf[network.land.order[end-100]] ≈ 7.855716706804882e11
-    @test ssf[network.land.order[end]] ≈ 2.1612469198596365e11
+    @test ssf[network.land.order[end-100]] ≈ 7.855716879739626e11
+    @test ssf[network.land.order[end]] ≈ 2.161246841709492e11
 end
 
 @testset "overland flow" begin
@@ -114,16 +112,16 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 2807.935527751757
+    @test sum(q) ≈ 2807.884971209105
     @test q[4061] ≈ 0.0016288040314320486
-    @test q[5617] ≈ 7.338204361667292
+    @test q[5617] ≈ 7.338169165884175
     @test q[network.river.order[end]] ≈ 0.00610520650626283
 end
 
 @testset "reservoir simple" begin
     res = model.lateral.river.reservoir
     @test res.outflow[2] ≈ 0.2174998592483153
-    @test res.inflow[2] ≈ 50.170880948187055
+    @test res.inflow[2] ≈ 50.170880189190626
     @test res.volume[2] ≈ 2.776162917050312e7
     @test res.precipitation[2] ≈ 0.1765228509902954
     @test res.evaporation[2] ≈ 0.5372688174247742
@@ -157,7 +155,7 @@ end
 # downstream from pit is at river index 1739, CartesianIndex(142, 85)
 @testset "river flow at and downstream of pit" begin
     q = model.lateral.river.q_av
-    @test q[1765] ≈ 8.133435755782642
+    @test q[1765] ≈ 8.133435566601927
     @test q[1739] ≈ 0.008996045895155833
 end
 
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 963f71860..4596cbf99 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -21,7 +21,6 @@ end
     sbm = model.vertical
 
     @test model.clock.iteration == 2
-    @test sbm.altitude[1] == 1.937000036239624
     @test sbm.θₛ[1] == 0.44999998807907104
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
@@ -33,7 +32,6 @@ model = Wflow.update_sbm_gwf(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.altitude[1] == 1.937000036239624
     @test sbm.θₛ[1] == 0.44999998807907104
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
@@ -47,10 +45,10 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 0.03515256649675796
-    @test q[6] ≈ 0.00795266793733648
-    @test q[13] ≈ 0.0005980549864213564
-    @test q[network.river.order[end]] ≈ 0.008482646477950898
+    @test sum(q) ≈ 0.03515257228971982
+    @test q[6] ≈ 0.007952670041402076
+    @test q[13] ≈ 0.0005980549969548235
+    @test q[network.river.order[end]] ≈ 0.008482648782941154
 end
 
 @testset "groundwater" begin
diff --git a/test/sbm_config.toml b/test/sbm_config.toml
index 48f4bdff0..76bd598c2 100644
--- a/test/sbm_config.toml
+++ b/test/sbm_config.toml
@@ -61,24 +61,22 @@ forcing = [
 cyclic = ["vertical.leaf_area_index"]
 
 [input.vertical]
-altitude = "wflow_dem"
 c = "c"
 cf_soil = "cf_soil"
 cfmax = "Cfmax"
 e_r = "EoverR"
+f = "f"
 infiltcappath = "InfiltCapPath"
 infiltcapsoil = "InfiltCapSoil"
 kext = "Kext"
 "kv₀" = "KsatVer"
 leaf_area_index = "LAI"
-m = "M"
 maxleakage = "MaxLeakage"
 pathfrac = "PathFrac"
 potential_evaporation = "PET"
 precipitation = "P"
 rootdistpar = "rootdistpar"
 rootingdepth = "RootingDepth"
-soilminthickness = "SoilMinThickness"
 soilthickness = "SoilThickness"
 specific_leaf = "Sl"
 storage_wood = "Swood"
@@ -86,7 +84,6 @@ temperature = "TEMP"
 tt = "TT"
 tti = "TTI"
 ttm = "TTM"
-w_soil = "wflow_soil"
 water_holding_capacity = "WHC"
 waterfrac = "WaterFrac"
 "θᵣ" = "thetaR"
diff --git a/test/sbm_gwf_config.toml b/test/sbm_gwf_config.toml
index 2fb4ace4c..2f82850c0 100644
--- a/test/sbm_gwf_config.toml
+++ b/test/sbm_gwf_config.toml
@@ -40,6 +40,7 @@ path_static = "data/staticmaps-sbm-groundwater.nc"
 ldd = "wflow_ldd"
 river_location = "wflow_river"
 subcatchment = "wflow_subcatch"
+altitude = "wflow_dem"
 
 # specify the internal IDs of the parameters which vary over time
 # the external name mapping needs to be below together with the other mappings
@@ -49,7 +50,6 @@ forcing = [
 ]
 
 [input.vertical]
-altitude = "wflow_dem"
 canopygapfraction = "canopygapfraction"
 cmax = "cmax"
 e_r = "EoverR"