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

Issue:
`SBM` evaporation and fractions (open water, vegetation, soil surface)
### Feature type

Changing existing functionality

### Improvement Description

For evaporation in `SBM` we use potential evaporation with an optional correction factor `et_reftopot`. Then different processes are computed in the following order:  1) interception, 2) open water evaporation (river and land), 3) soil evaporation and 4) transpiration. The remaining potential evaporation is computed for step 1 and 2, and split based on the canopy gap fraction into remaining potential evaporation for step 3 and 4. 

An improvement is to base the potential evaporation on the different fractions (which is now done partly) of open water (`riverfrac` and `waterfrac`), vegetation and bare soil for each process, instead of working with remaining potential evaporation (bookkeeping). We can then also make use of different factors (e.g. a soil factor for bare soil evaporation, crop factor for transpiration), and keep better track of the different fractions.  

### Implementation Description

_No response_

### Additional Context

_No response_


Discussion / hints:

Fix patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 65f774665..ec2e38872 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -22,6 +22,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `to_river` variable from overland flow and lateral subsurface flow was not added to the
   inflow of these locations.
 - Close netCDF `NCDataset` with state variables in extended BMI function `save_state`.
+- For the computation of Gash interception model parameter `e_r` multiply the precipitation
+  input with the canopy fraction (this was only done for the potential evapotranspiration
+  input).
+
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was not
@@ -46,6 +50,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   When in the Python version of wflow\_sbm the kinematic wave for surface flow was split
   into a river and land component, the lower limit was removed for river runoff
   (`net_runoff_river`), but not for land runoff.
+- Always use fractions for the computation of potential evapotranspiration (interception and
+  transpiration) and potential evaporation (open water and soil). Replaced variable
+  `et_reftopot` by crop coefficient `kc` (use of `et_reftopot` has been deprecated).
 
 ### Added
 - Total water storage as an export variable for `SBM` concept. This is the total water stored
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 214fe3c17..840eee979 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -64,7 +64,7 @@ profile `kv` is used and `z_layered` is required as input.
 | **`rootdistpar`** | controls how roots are linked to water table | - | -500.0  |
 | **`cap_hmax`** | water depth beyond which capillary flux ceases | mm | 2000.0  |
 | **`cap_n`** | coefficient controlling capillary rise | - | 2.0  |
-| **`et_reftopot`** | multiplication factor to correct reference evaporation | - | 1.0  |
+| **`kc`** | crop coefficient Kc | - | 1.0  |
 | **`sl`** (`specific_leaf`) | specific leaf storage  | mm | - |
 | **`swood`** (`storage_wood`) | storage woody part of vegetation | mm | - |
 | **`kext`** | extinction coefficient (to calculate canopy gap fraction) | - | - |
@@ -89,8 +89,8 @@ profile `kv` is used and `z_layered` is required as input.
 | `canopystorage`|  canopy storage | mm | - |
 |**`precipitation`** | precipitation | mm Δt``^{-1}``|  - |
 | **`temperature`** | temperature | ᵒC | - |
-| **`potential_evaporation`** | potential evaporation | mm Δt``^{-1}`` | - |
-| `pottrans_soil` | interception subtracted from potential evaporation) | mm Δt``^{-1}`` | - |
+| **`potential_evaporation`** | potential reference evapotranspiration | mm Δt``^{-1}`` | - |
+| `pottrans` | interception subtracted from potential evapotranspiration | mm Δt``^{-1}`` | - |
 | `transpiration` | transpiration | mm Δt``^{-1}`` | - |
 | `ae_ustore` | actual evaporation from unsaturated store | mm Δt``^{-1}`` | - |
 | `interception` | interception loss by evaporation | mm Δt``^{-1}`` | - |
diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md
index 159258e57..6a7744339 100644
--- a/docs/src/model_docs/vertical/sbm.md
+++ b/docs/src/model_docs/vertical/sbm.md
@@ -134,22 +134,32 @@ The extinction coefficient `kext` can be related to land cover.
 
 ## [Evaporation](@id evap)
 
-The wflow\_sbm model assumes the input to be potential evaporation. A multiplication factor
-(`et_reftopot`, set to 1 by default) is present to correct the input evaporation if
-required.
-
-The potential evaporation left over after interception and open water evaporation (rivers
-and water bodies) is split in potential soil evaporation and potential transpiration based
-on the canopy gap fraction (assumed to be identical to the amount of bare soil).
+The wflow\_sbm model assumes the input to be potential reference evapotranspiration. A crop
+coefficient (`kc`, set to 1 by default) is used to convert the potential evapotranspiration
+rate of a reference crop fully covering the soil to the potential evapotranspiration rate of
+vegetation (natural and agricultural) fully covering the soil. The crop coefficient `kc` of
+wflow\_sbm is used for a surface completely covered by vegetation, and does not include the
+effect of growing stages of vegetation and soil cover. These effects are handled separately
+through the use of the canopy gap fraction. 
+
+It is assumed that the potential evaporation rate of intercepted water by vegetation is
+equal to the potential evapotranspiration rate of vegetation (fully covering the soil)
+multiplied by the canopy fraction. The potential evapotranspiration rate left over after
+interception is available for transpiration. For potential open water evaporation (river and
+water bodies) the potential reference evapotranspiration rate is used (multipled by the
+river fraction `riverfrac`, and open water fraction `waterfrac`). Also for potential soil
+evaporation the potential reference evapotranspiration rate is used, multiplied by the
+canopy gap fraction corrected by the sum of total water fraction (`riverfrac` and
+`waterfrac`) and the fraction covered by a glacier (`glacierfrac`).
 
 ### Bare soil evaporation
 
 If there is only one soil layer present in the wflow\_sbm model, the bare soil evaporation
 is scaled according to the wetness of the soil layer. The fraction of bare soil is assumed
-to be equal to the fraction not covered by the canopy (`canopygapfraction`). When the soil
-is fully saturated, evaporation is set to equal the potential evaporation. When the soil is
-not fully saturated, actual evaporation decreases linearly with decreasing soil moisture
-values, as indicated by the figure below.
+to be equal to the fraction not covered by the canopy (`canopygapfraction`) corrected by the
+total water fraction. When the soil is fully saturated, evaporation is set to equal the
+potential reference evaporation. When the soil is not fully saturated, actual evaporation
+decreases linearly with decreasing soil moisture values, as indicated by the figure below.
 
 ![soil_evap](../../images/soil_evap.png)
 
@@ -182,10 +192,10 @@ index of the vector that contains all active cells within the spatial model doma
 
 ```julia
     # transpiration from saturated store
-    wetroots = scurve(sbm.zi[i], rootingdepth, 1.0, sbm.rootdistpar[i])
-    actevapsat = min(pottrans * wetroots, satwaterdepth)
+    wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
+    actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
     satwaterdepth = satwaterdepth - actevapsat
-    restpottrans = pottrans - actevapsat
+    restpottrans = sbm.pottrans[i] - actevapsat
 ```
 
 ![soil_wetroots](../../images/soil_wetroots.png)
diff --git a/src/sbm.jl b/src/sbm.jl
index ba508ddb5..e0282cda3 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -49,8 +49,8 @@
     cap_hmax::Vector{T} | "mm"
     # Coefficient [-] controlling capillary rise
     cap_n::Vector{T} | "-"
-    # Multiplication factor [-] to correct
-    et_reftopot::Vector{T} | "-"
+    # Crop coefficient Kc [-]
+    kc::Vector{T} | "-"
     # Brooks-Corey power coefﬁcient [-] for each soil layer
     c::Vector{SVector{N,T}} | "-"
     # Stemflow [mm Δt⁻¹]
@@ -84,10 +84,10 @@
     precipitation::Vector{T}
     # Temperature [ᵒC]
     temperature::Vector{T} | "°C"
-    # Potential evapotranspiration [mm Δt⁻¹]
+    # Potential reference evapotranspiration [mm Δt⁻¹]
     potential_evaporation::Vector{T}
-    # Potential transpiration, open water, river and soil evaporation (after subtracting interception from potential_evaporation)
-    pottrans_soil::Vector{T}
+    # Potential transpiration (after subtracting interception from potential_evaporation)
+    pottrans::Vector{T}
     # Transpiration [mm Δt⁻¹]
     transpiration::Vector{T}
     # Actual evaporation from unsaturated store [mm Δt⁻¹]
@@ -469,8 +469,21 @@ function initialize_sbm(nc, config, riverfrac, inds)
     cap_hmax =
         ncread(nc, config, "vertical.cap_hmax"; sel = inds, defaults = 2000.0, type = Float)
     cap_n = ncread(nc, config, "vertical.cap_n"; sel = inds, defaults = 2.0, type = Float)
-    et_reftopot =
-        ncread(nc, config, "vertical.et_reftopot"; sel = inds, defaults = 1.0, type = Float)
+    kc = ncread(
+        nc,
+        config,
+        "vertical.kc";
+        alias = "vertical.et_reftopot",
+        sel = inds,
+        defaults = 1.0,
+        type = Float,
+    )
+    if haskey(config.input.vertical, "et_reftopot")
+        @warn string(
+            "The `et_reftopot` key in `[input.vertical]` is now called ",
+            "`kc`. Please update your TOML file.",
+        )
+    end
 
     cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
 
@@ -591,7 +604,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         rootdistpar = rootdistpar,
         cap_hmax = cap_hmax,
         cap_n = cap_n,
-        et_reftopot = et_reftopot,
+        kc = kc,
         c = svectorscopy(c, Val{maxlayers}()),
         stemflow = fill(mv, n),
         throughfall = fill(mv, n),
@@ -609,7 +622,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         precipitation = fill(mv, n),
         temperature = fill(mv, n),
         potential_evaporation = fill(mv, n),
-        pottrans_soil = fill(mv, n),
+        pottrans = fill(mv, n),
         transpiration = fill(mv, n),
         ae_ustore = fill(mv, n),
         interception = fill(mv, n),
@@ -688,22 +701,19 @@ function update_until_snow(sbm::SBM, config)
         if do_lai
             cmax = sbm.sl[i] * sbm.leaf_area_index[i] + sbm.swood[i]
             canopygapfraction = exp(-sbm.kext[i] * sbm.leaf_area_index[i])
-            ewet =
-                (1.0 - exp(-sbm.kext[i] * sbm.leaf_area_index[i])) *
-                sbm.potential_evaporation[i]
+            canopyfraction = 1.0 - canopygapfraction
+            ewet = canopyfraction * sbm.potential_evaporation[i] * sbm.kc[i]
             e_r =
                 sbm.precipitation[i] > 0.0 ?
-                min(0.25, ewet / max(0.0001, sbm.precipitation[i])) : 0.0
+                min(0.25, ewet / max(0.0001, canopyfraction * sbm.precipitation[i])) : 0.0
         else
             cmax = sbm.cmax[i]
             canopygapfraction = sbm.canopygapfraction[i]
             e_r = sbm.e_r[i]
         end
 
-        potential_evaporation = sbm.potential_evaporation[i] * sbm.et_reftopot[i]
-        # should we include tempcor in SBM?
-        # potential_evaporation = PotenEvap #??
-
+        canopy_potevap =
+            sbm.kc[i] * sbm.potential_evaporation[i] * (1.0 - canopygapfraction)
         if Second(sbm.Δt) >= Hour(23)
             throughfall, interception, stemflow, canopystorage = rainfall_interception_gash(
                 cmax,
@@ -711,19 +721,19 @@ function update_until_snow(sbm::SBM, config)
                 canopygapfraction,
                 sbm.precipitation[i],
                 sbm.canopystorage[i],
-                potential_evaporation,
+                canopy_potevap,
             )
-            pottrans_soil = max(0.0, potential_evaporation - interception) # now in mm
+            pottrans = max(0.0, canopy_potevap - interception) # now in mm
         else
             netinterception, throughfall, stemflow, leftover, interception, canopystorage =
                 rainfall_interception_modrut(
                     sbm.precipitation[i],
-                    potential_evaporation,
+                    canopy_potevap,
                     sbm.canopystorage[i],
                     canopygapfraction,
                     cmax,
                 )
-            pottrans_soil = max(0.0, leftover)  # now in mm
+            pottrans = max(0.0, leftover)  # now in mm
         end
 
         if modelsnow
@@ -748,7 +758,7 @@ function update_until_snow(sbm::SBM, config)
         sbm.interception[i] = interception
         sbm.stemflow[i] = stemflow
         sbm.throughfall[i] = throughfall
-        sbm.pottrans_soil[i] = pottrans_soil
+        sbm.pottrans[i] = pottrans
         if modelsnow
             sbm.snow[i] = snow
             sbm.snowwater[i] = snowwater
@@ -806,18 +816,20 @@ function update_until_recharge(sbm::SBM, config)
 
         ae_openw_r = min(
             sbm.waterlevel_river[i] * sbm.riverfrac[i],
-            sbm.riverfrac[i] * sbm.pottrans_soil[i],
+            sbm.riverfrac[i] * sbm.potential_evaporation[i],
         )
         ae_openw_l = min(
             sbm.waterlevel_land[i] * sbm.waterfrac[i],
-            sbm.waterfrac[i] * sbm.pottrans_soil[i],
+            sbm.waterfrac[i] * sbm.potential_evaporation[i],
         )
 
-        restevap = sbm.pottrans_soil[i] - ae_openw_r - ae_openw_l
-
-        # evap available for soil evaporation and transpiration
-        potsoilevap = restevap * sbm.canopygapfraction[i]
-        pottrans = restevap * (1.0 - sbm.canopygapfraction[i])
+        # evap available for soil evaporation
+        soilevap_fraction = max(
+            sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.waterfrac[i] -
+            sbm.glacierfrac[i],
+            0.0,
+        )
+        potsoilevap = soilevap_fraction * sbm.potential_evaporation[i]
 
         # Calculate the initial capacity of the unsaturated store
         ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
@@ -926,9 +938,9 @@ function update_until_recharge(sbm::SBM, config)
 
         # transpiration from saturated store
         wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
-        actevapsat = min(pottrans * wetroots, satwaterdepth)
+        actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
         satwaterdepth = satwaterdepth - actevapsat
-        restpottrans = pottrans - actevapsat
+        restpottrans = sbm.pottrans[i] - actevapsat
 
         # actual transpiration from ustore
         actevapustore = 0.0


Test patch:
diff --git a/server/test/client.jl b/server/test/client.jl
index 45051e81f..9c06e66fd 100644
--- a/server/test/client.jl
+++ b/server/test/client.jl
@@ -15,8 +15,7 @@ function request(message)
 end
 
 @testset "initialization and time functions" begin
-    msg =
-        (fn = "initialize", config_file = joinpath(@__DIR__, "sbm_config.toml"))
+    msg = (fn = "initialize", config_file = joinpath(@__DIR__, "sbm_config.toml"))
     @test request(msg) == Dict("status" => "OK")
     @test request((fn = "get_end_time",)) == Dict("end_time" => 2678400)
     @test request((fn = "get_start_time",)) == Dict("start_time" => 0)
@@ -67,7 +66,7 @@ vwc_1_size = 0
     vwc_1_size = Int(vwc_1_nbytes / vwc_1_itemsize)
     @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
-    @test mean(request(msg)["value"]) ≈ 278.1510965581235
+    @test mean(request(msg)["value"]) ≈ 277.8362161218515
     msg = (fn = "get_value_ptr", name = "vertical.θₛ")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
@@ -77,7 +76,7 @@ vwc_1_size = 0
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.1901434445889123, 2.6778265820894545, 3.470059871798756]
+          [2.196562540141252f0, 2.6853997213000866f0, 3.4814142040324985f0]
     msg = (fn = "set_value", name = "vertical.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
@@ -97,7 +96,7 @@ vwc_1_size = 0
     )
     @test request(msg)["value_at_indices"] == [250.0, 350.0, 300.0]
     msg = (fn = "get_value", name = "vertical.vwc[1]", dest = fill(0.0, vwc_1_size))
-    @test mean(request(msg)["value"]) ≈ 0.1845159140308566f0
+    @test mean(request(msg)["value"]) ≈ 0.18600013563085036f0
     msg = (
         fn = "get_value_at_indices",
         name = "vertical.vwc[1]",
@@ -105,7 +104,7 @@ vwc_1_size = 0
         inds = [1, 2, 3],
     )
     @test request(msg)["value_at_indices"] ≈
-          [0.12089607119560242f0, 0.11967185322648062f0, 0.14503555864288548f0]
+          [0.12089607119560242f0, 0.11968416924304527f0, 0.14602328618707333f0]
     msg = (fn = "set_value", name = "vertical.vwc[1]", src = fill(0.3, vwc_1_size))
     @test request(msg) == Dict("status" => "OK")
     msg = (fn = "get_value", name = "vertical.vwc[1]", dest = fill(0.0, vwc_1_size))
diff --git a/test/bmi.jl b/test/bmi.jl
index 66a97ad46..75470c4d5 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -57,7 +57,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test_throws ErrorException BMI.get_value_ptr(model, "vertical.")
             dest = zeros(Float, size(model.vertical.zi))
             BMI.get_value(model, "vertical.zi", dest)
-            @test mean(dest) ≈ 276.3767651555451
+            @test mean(dest) ≈ 276.16325589542333
             @test BMI.get_value_at_indices(
                 model,
                 "vertical.vwc[1]",
@@ -75,7 +75,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
                 "lateral.river.q",
                 zeros(Float, 3),
                 [1, 100, 5617],
-            ) ≈ [0.6211503865184697, 5.219305686635002, 0.026163746306482282]
+            ) ≈ [0.623325399343309, 5.227139951657074, 0.02794287432778194]
             BMI.set_value(model, "vertical.zi", fill(300.0, length(model.vertical.zi)))
             @test mean(
                 BMI.get_value(model, "vertical.zi", zeros(Float, size(model.vertical.zi))),
@@ -150,10 +150,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "recharge part of SBM" begin
             sbm = model.vertical
-            @test sbm.interception[1] ≈ 0.6299999952316284f0
+            @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.1912317735997524f0
-            @test sbm.recharge[5] ≈ -0.0727941579914808f0
+            @test sbm.snow[1] ≈ 3.484789961176288f0
+            @test sbm.recharge[5] ≈ -0.0f0
             @test sbm.zi[5] ≈ 300.0f0
         end
 
@@ -174,10 +174,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         @testset "SBM after subsurface flow" begin
             sbm = model.vertical
             sub = model.lateral.subsurface
-            @test sbm.interception[1] ≈ 0.6299999952316284
+            @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.1912317735997524f0
-            @test sbm.recharge[5] ≈ -0.0727941579914808f0
+            @test sbm.snow[1] ≈ 3.484789961176288f0
+            @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 250.0f0
             @test sub.zi[5] ≈ 0.25f0
             @test sub.exfiltwater[1] ≈ 1.0f-5
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 148e524b6..e9dcd6ec3 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -14,53 +14,53 @@ flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
     row = csv_first_row(model.writer.csv_path)
 
     @test row.time == DateTime("2000-01-02T00:00:00")
-    @test row.Q ≈ 8.1396354649197f0
+    @test row.Q ≈ 8.15299947254324f0
     @test row.volume ≈ 2.7535003939625636f7
     @test row.temp_bycoord ≈ 2.390000104904175f0
-    @test row.vwc_layer2_bycoord ≈ 0.25901943991019094f0
+    @test row.vwc_layer2_bycoord ≈ 0.25938809638672006f0
     @test row.temp_byindex ≈ 2.390000104904175f0
-    @test row.Q_6336050 ≈ 0.006160282939850074f0
-    @test row.Q_6336510 ≈ 0.029177309556712334f0
-    @test row.Q_6836100 ≈ 0.19621282450614713f0
-    @test row.Q_6336500 ≈ 0.006089112638001381f0
+    @test row.Q_6336050 ≈ 0.006583064321841488f0
+    @test row.Q_6336510 ≈ 0.029864230092642642f0
+    @test row.Q_6836100 ≈ 0.19995488963854305f0
+    @test row.Q_6336500 ≈ 0.006277726622788425f0
     @test row.Q_6836190 ≈ 0.0031262850749354237f0
-    @test row.Q_6336800 ≈ 0.007770868657277307f0
-    @test row.Q_6336900 ≈ 0.006403194169947582f0
-    @test row.Q_6336930 ≈ 0.08888787154163148f0
-    @test row.Q_6336910 ≈ 0.007071851236520184f0
-    @test row.Q_6136500 ≈ 0.0016367337487926633f0
-    @test row.Q_6136520 ≈ 0.002084670434294102f0
-    @test row.Q_6136150 ≈ 0.006095549758915344f0
-    @test row.Q_6136151 ≈ 0.007643634432992056f0
-    @test row.Q_6136160 ≈ 3.9199531493174726f0
-    @test row.Q_6136202 ≈ 1.4125847550988493f0
-    @test row.recharge_1 ≈ -0.05653226176238641f0
+    @test row.Q_6336800 ≈ 0.008278375560053742f0
+    @test row.Q_6336900 ≈ 0.0066141980189014385f0
+    @test row.Q_6336930 ≈ 0.09141703511009937f0
+    @test row.Q_6336910 ≈ 0.007475453481320056f0
+    @test row.Q_6136500 ≈ 0.001834989281902289f0
+    @test row.Q_6136520 ≈ 0.0022266031120691397f0
+    @test row.Q_6136150 ≈ 0.006310361139139334f0
+    @test row.Q_6136151 ≈ 0.007946301730645885f0
+    @test row.Q_6136160 ≈ 3.927719795530719f0
+    @test row.Q_6136202 ≈ 1.4162246003743886f0
+    @test row.recharge_1 ≈ -0.002257181032501202f0
 end
 
 @testset "NetCDF scalar output" begin
     ds = model.writer.dataset_scalar
     @test ds["time"][1] == DateTime("2000-01-02T00:00:00")
     @test ds["Q"][:][1:20] ≈ [
-        0.73993874f0,
-        1.4125848f0,
-        1.4389194f0,
-        1.4036233f0,
-        5.7276225f0,
-        2.7586424f0,
-        2.1080604f0,
-        4.1026044f0,
-        0.008365918f0,
-        3.919953f0,
-        4.0615087f0,
-        0.006050462f0,
-        0.0076436345f0,
-        0.00780229f0,
-        0.00346879f0,
-        0.7066179f0,
-        0.0022446192f0,
-        1.3323202f0,
-        3.8271446f0,
-        1.6729931f0,
+        0.7425387f0,
+        1.4162246f0,
+        1.4425076f0,
+        1.4044669f0,
+        5.738109f0,
+        2.7616737f0,
+        2.1128905f0,
+        4.105428f0,
+        0.008651769f0,
+        3.9277198f0,
+        4.069447f0,
+        0.006356805f0,
+        0.007946302f0,
+        0.008135906f0,
+        0.0037393502f0,
+        0.70888275f0,
+        0.0024000728f0,
+        1.3347782f0,
+        3.8374817f0,
+        1.676597f0,
     ]
     @test ds["Q_gauges"].attrib["cf_role"] == "timeseries_id"
     @test ds["temp_index"][:] ≈ [2.39f0]
@@ -77,13 +77,13 @@ end
 
     @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
     @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
-    @test mean(sbm.runoff) ≈ 0.0416027711985709f0
-    @test mean(sbm.soilevap) ≈ 0.00616346765710303f0
-    @test mean(sbm.actevap) ≈ 0.5393128846181061f0
-    @test mean(sbm.actinfilt) ≈ 1.4860556930083888f0
-    @test sbm.snow[5] ≈ 3.592840840467347f0
-    @test mean(sbm.snow) ≈ 0.035197394854698104f0
-    @test sbm.total_storage[50063] ≈ 559.70849973344f0
+    @test mean(sbm.runoff) ≈ 0.04177459898728149f0
+    @test mean(sbm.soilevap) ≈ 0.02122698830889417f0
+    @test mean(sbm.actevap) ≈ 0.33545623834952154f0
+    @test mean(sbm.actinfilt) ≈ 1.6444774688444848f0
+    @test sbm.snow[5] ≈ 3.768513390588815f0
+    @test mean(sbm.snow) ≈ 0.038019723676094325f0
+    @test sbm.total_storage[50063] ≈ 559.9035608052374f0
     @test sbm.total_storage[429] ≈ 597.4578475404879f0 # river cell
 end
 
@@ -94,28 +94,28 @@ model = Wflow.run_timestep(model)
     sbm = model.vertical
     @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
     @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
-    @test mean(sbm.net_runoff) ≈ 0.23657364759232993f0
-    @test mean(sbm.runoff) ≈ 0.2369017651524159f0
-    @test mean(sbm.soilevap) ≈ 0.015275916086112940
-    @test mean(sbm.actevap) ≈ 0.29716147685525746f0
-    @test mean(sbm.actinfilt) ≈ 0.08829600284076534f0
-    @test sbm.snow[5] ≈ 3.667748983774868f0
-    @test mean(sbm.snow) ≈ 0.03209680141455693f0
-    @test sbm.total_storage[50063] ≈ 559.7935411649405f0
-    @test sbm.total_storage[429] ≈ 617.0062092646873f0 # river cell
+    @test mean(sbm.net_runoff) ≈ 0.23710522236947756f0
+    @test mean(sbm.runoff) ≈ 0.23747368431143517f0
+    @test mean(sbm.soilevap) ≈ 0.018750808322054897f0
+    @test mean(sbm.actevap) ≈ 0.27379330211866343f0
+    @test mean(sbm.actinfilt) ≈ 0.08863102527394363f0
+    @test sbm.snow[5] ≈ 3.843412524052313f0
+    @test mean(sbm.snow) ≈ 0.03461317061870949f0
+    @test sbm.total_storage[50063] ≈ 560.0152135062889f0
+    @test sbm.total_storage[429] ≈ 617.2238533241972f0 # river cell
 end
 
 @testset "subsurface flow" begin
     ssf = model.lateral.subsurface.ssf
     @test sum(ssf) ≈ 6.370399148012509f7
     @test ssf[network.land.order[1]] ≈ 7.169036749244327f2
-    @test ssf[network.land.order[end-100]] ≈ 2335.2465707069578f0
+    @test ssf[network.land.order[end-100]] ≈ 2333.801056570759f0
     @test ssf[network.land.order[end]] ≈ 288.19428729403944f0
 end
 
 @testset "overland flow" begin
     q = model.lateral.land.q_av
-    @test sum(q) ≈ 290.5520014030802f0
+    @test sum(q) ≈ 291.27639107427285f0
     @test q[26625] ≈ 0.0
     @test q[39308] ≈ 0.0
     @test q[network.land.order[end]] ≈ 1.0f-30
@@ -123,10 +123,10 @@ end
 
 @testset "river flow" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3615.2084127447865f0
-    @test q[1622] ≈ 0.0005986272622329333f0
-    @test q[43] ≈ 12.036342425160155f0
-    @test q[network.river.order[end]] ≈ 0.03835913312643948f0
+    @test sum(q) ≈ 3622.7369292570543f0
+    @test q[1622] ≈ 0.0006497468064774366f0
+    @test q[43] ≈ 12.05767242907667f0
+    @test q[network.river.order[end]] ≈ 0.039200124520463835f0
 end
 
 @testset "reservoir simple" begin
@@ -166,10 +166,10 @@ end
 
 @testset "river flow at basin outlets and downstream of one pit" begin
     q = model.lateral.river.q_av
-    @test q[4009] ≈ 8.60480399680283f0 # pit/ outlet, CartesianIndex(141, 228)
+    @test q[4009] ≈ 8.51907041734622f0 # pit/ outlet, CartesianIndex(141, 228)
     @test q[4020] ≈ 0.006779014715290862f0 # downstream of pit 4009, CartesianIndex(141, 229)
-    @test q[2508] ≈ 150.15321264134985f0 # pit/ outlet
-    @test q[5808] ≈ 0.12625654862968252f0 # pit/ outlet
+    @test q[2508] ≈ 150.28398167251638f0 # pit/ outlet
+    @test q[5808] ≈ 0.12419895007970105f0 # pit/ outlet
 end
 
 # test changing forcing and cyclic LAI parameter
@@ -212,7 +212,7 @@ model = Wflow.run_timestep(model)
 
 @testset "river inflow (cyclic)" begin
     @test model.lateral.river.inflow[44] ≈ 0.75
-    @test model.lateral.river.q_av[44] ≈ 10.698591283662008
+    @test model.lateral.river.q_av[44] ≈ 10.71846407068599
 end
 
 # test fixed forcing (precipitation = 2.5)
@@ -248,14 +248,14 @@ model = Wflow.run_timestep(model)
 
 @testset "river flow and depth (local inertial)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3910.2095742376546f0
-    @test q[1622] ≈ 6.0094181857060604f-5
-    @test q[43] ≈ 11.900372477232786f0
-    @test q[501] ≈ 3.536628093804679f0
+    @test sum(q) ≈ 3919.6025219496014f0
+    @test q[1622] ≈ 7.31010246736994f-5
+    @test q[43] ≈ 11.92153120707289f0
+    @test q[501] ≈ 3.5736389982451895f0
     h = model.lateral.river.h_av
-    @test h[1622] ≈ 0.001809965063947279f0
-    @test h[43] ≈ 0.43627044208669874f0
-    @test h[501] ≈ 0.05669956233680719f0
+    @test h[1622] ≈ 0.001987887644883981f0
+    @test h[43] ≈ 0.4366415244811759f0
+    @test h[501] ≈ 0.057265962518284294f0
     q_channel = model.lateral.river.q_channel_av
     @test q ≈ q_channel
 end
@@ -271,21 +271,21 @@ model = Wflow.run_timestep(model)
 
 @testset "river and overland flow and depth (local inertial)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 2375.1186861861243f0
-    @test q[1622] ≈ 6.011407534125278f-5
-    @test q[43] ≈ 5.358152280519331f0
-    @test q[501] ≈ 1.5878151534724314f0
+    @test sum(q) ≈ 2380.64389229669f0
+    @test q[1622] ≈ 7.322956970529551f-5
+    @test q[43] ≈ 5.361283165612762f0
+    @test q[501] ≈ 1.6021771576366957f0
     h = model.lateral.river.h_av
-    @test h[1622] ≈ 0.0018101771426228585f0
+    @test h[1622] ≈ 0.0019891342000364796f0
     @test h[43] ≈ 0.3003008810153667f0
-    @test h[501] ≈ 0.03162351626631113f0
+    @test h[501] ≈ 0.031925992442532f0
     qx = model.lateral.land.qx
     qy = model.lateral.land.qy
-    @test qx[[26, 35, 631]] ≈ [0.18613687016733824f0, 0.0004519163131931592f0, 0.0f0]
-    @test qy[[26, 35, 631]] ≈ [0.12681702046955443f0, 1.7210193779889194f0, 0.0f0]
+    @test qx[[26, 35, 631]] ≈ [0.18614776104106373f0, 0.029502872625766417f0, 0.0f0]
+    @test qy[[26, 35, 631]] ≈ [0.12757214437549858f0, 1.7212079599401755f0, 0.0f0]
     h = model.lateral.land.h
     @test h[[26, 35, 631]] ≈
-          [0.07341443653334193f0, 0.009152294150993293f0, 0.0006875940563996746f0]
+          [0.07361854999908582f0, 0.009155393111676267f0, 0.0007258741013439351f0]
 end
 Wflow.close_files(model, delete_output = false)
 
@@ -401,19 +401,19 @@ model = Wflow.run_timestep(model)
 
 @testset "river flow (local inertial) with floodplain schematization simulation" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3898.719057830299f0
-    @test q[1622] ≈ 6.0094627478450016f-5
-    @test q[43] ≈ 11.900372477232796f0
-    @test q[501] ≈ 3.470259878228359f0
-    @test q[5808] ≈ 0.002199141435542889f0
+    @test sum(q) ≈ 3908.039208613999f0
+    @test q[1622] ≈ 7.310102468091527f-5
+    @test q[43] ≈ 11.921531207072922f0
+    @test q[501] ≈ 3.5061516913374717f0
+    @test q[5808] ≈ 0.0022258226497210145f0
     h = model.lateral.river.h_av
-    @test h[1622] ≈ 0.0018099697988149294f0
-    @test h[43] ≈ 0.4362704420867342f0
-    @test h[501] ≈ 0.05610231297517167f0
-    @test h[5808] ≈ 0.005901095907129176f0
+    @test h[1622] ≈ 0.001987887580593841f0
+    @test h[43] ≈ 0.436641524481545f0
+    @test h[501] ≈ 0.05665942153713204f0
+	@test h[5808] ≈ 0.005933025055544241f0
 end
 
-# set boundary condition local inertial routing from NetCDF file
+# set boundary condition local inertial routing from netCDF file
 config.input.lateral.river.riverlength_bc = "riverlength_bc"
 config.input.lateral.river.riverdepth_bc = "riverdepth_bc"
 model = Wflow.initialize_sbm_model(config)
@@ -422,15 +422,15 @@ model = Wflow.run_timestep(model)
 
 @testset "change boundary condition for local inertial routing (including floodplain)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 3899.01831081148f0
-    @test q[1622] ≈ 6.008529317397127f-5
-    @test q[43] ≈ 11.900394994314782f0
-    @test q[501] ≈ 3.4702817507735366f0
-    @test q[5808] ≈ 0.060081151400640784f0
+    @test sum(q) ≈ 3908.2493947503826f0
+    @test q[1622] ≈ 7.310183576829848f-5
+    @test q[43] ≈ 11.921531207085014f0
+    @test q[501] ≈ 3.5061513868779763f0
+    @test q[5808] ≈ 0.060518234525259465f0
     h = model.lateral.river.h_av
-    @test h[1622] ≈ 0.001809870969211623f0
-    @test h[43] ≈ 0.4362709200147433f0
-    @test h[501] ≈ 0.05610440991535934f0
+    @test h[1622] ≈ 0.0019878952928530183f0
+    @test h[43] ≈ 0.4366415249636809f0
+    @test h[501] ≈ 0.05665929962422606f0
     @test h[5808] ≈ 2.0000006940603936f0
 end
 Wflow.close_files(model, delete_output = false)
@@ -533,9 +533,9 @@ Wflow.close_files(model, delete_output = false)
     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
+        @test sum(q) ≈ 3158.1868156296678f0
+        @test q[1622] ≈ 0.0005972577112819149f0
+        @test q[43] ≈ 10.013363662625276f0
     end
 
     Wflow.close_files(model, delete_output = false)
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 629e38389..9d2ba53ac 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -46,12 +46,12 @@ end
 @testset "river domain (kinematic wave)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.034846861707851576f0
-    @test q[6] ≈ 0.007866587223009643f0
-    @test river.volume[6] ≈ 4.474556355925399f0
-    @test river.inwater[6] ≈ 0.00036392604840276924f0
-    @test q[13] ≈ 0.0005952292000597283f0
-    @test q[network.river.order[end]] ≈ 0.008387430132453135f0
+    @test sum(q) ≈ 0.034904078613089404f0
+    @test q[6] ≈ 0.007881734464420193f0
+    @test river.volume[6] ≈ 4.479866566743959f0
+    @test river.inwater[6] ≈ 0.00036640965826798314f0
+    @test q[13] ≈ 0.0005958036525897934f0
+    @test q[network.river.order[end]] ≈ 0.008405059891856978f0
 end
 
 @testset "groundwater" begin
@@ -64,7 +64,7 @@ end
         1.803682542252168f0,
         1.4049539807162532f0,
     ]
-    @test gw.river.flux[1] ≈ -50.46515313302901f0
+    @test gw.river.flux[1] ≈ -50.568389286368145f0
     @test gw.drain.flux[1] ≈ 0.0
     @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
 end
@@ -93,12 +93,12 @@ model = Wflow.run_timestep(model)
 @testset "river domain (local inertial)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
-    @test sum(q) ≈ 0.026752212882189867f0
-    @test q[6] ≈ 0.005962226036202787f0
-    @test river.volume[6] ≈ 7.498085951460259f0
-    @test river.inwater[6] ≈ 0.00017741362380231486f0
-    @test q[13] ≈ 0.0004586401133311585f0
-    @test q[5] ≈ 0.006312330297863098f0
+    @test sum(q) ≈ 0.026802342283209452f0
+    @test q[6] ≈ 0.005975747422547475f0
+    @test river.volume[6] ≈ 7.509678393246221f0
+    @test river.inwater[6] ≈ 0.00017989723366752878f0
+    @test q[13] ≈ 0.0004590884299495001f0
+    @test q[5] ≈ 0.006328157455390906f0
 end
 Wflow.close_files(model, delete_output = false)
 
@@ -123,14 +123,14 @@ model = Wflow.run_timestep(model)
 
 @testset "river and land domain (local inertial)" begin
     q = model.lateral.river.q_av
-    @test sum(q) ≈ 0.026759512926399318f0
-    @test q[6] ≈ 0.005964050519182144f0
-    @test q[13] ≈ 0.00045874863609566906f0
-    @test q[5] ≈ 0.006314336504100363f0
+    @test sum(q) ≈ 0.026809644054064f0
+    @test q[6] ≈ 0.005977572310302902f0
+    @test q[13] ≈ 0.0004591969811223254f0
+    @test q[5] ≈ 0.006330164124991123f0
     h = model.lateral.river.h_av
-    @test h[6] ≈ 0.08023321813680799f0
-    @test h[5] ≈ 0.07761589477563689f0
-    @test h[13] ≈ 0.08226917799750776f0
+    @test h[6] ≈ 0.08033159784710718f0
+    @test h[5] ≈ 0.0777121612660625f0
+    @test h[13] ≈ 0.08236630657766124f0
     qx = model.lateral.land.qx
     qy = model.lateral.land.qy
     @test all(qx .== 0.0f0)