From 49478ee3cefc489a1be9d768cdabfb294c1e0917 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Sat, 30 Nov 2024 17:28:20 +0100
Subject: [PATCH] Remove model parameter `dt` From structs `LakeParameters`,
 `ReservoirParameters`, `LateralSsfParameters` and
 `GroundwaterExchangeParameters`.

---
 src/routing/lake.jl                   | 46 ++++++++++++---------------
 src/routing/reservoir.jl              | 32 ++++++++-----------
 src/routing/subsurface.jl             | 22 +++----------
 src/routing/surface_kinwave.jl        |  8 ++---
 src/routing/surface_local_inertial.jl |  9 +++---
 src/sbm_gwf_model.jl                  |  4 +--
 src/sbm_model.jl                      | 10 +++---
 test/reservoir_lake.jl                | 18 ++++-------
 test/routing_process.jl               |  2 +-
 9 files changed, 63 insertions(+), 88 deletions(-)

diff --git a/src/routing/lake.jl b/src/routing/lake.jl
index 842842f30..682e040da 100644
--- a/src/routing/lake.jl
+++ b/src/routing/lake.jl
@@ -1,6 +1,5 @@
 "Struct for storing lake model parameters"
 @get_units @grid_loc @with_kw struct LakeParameters{T}
-    dt::T                                           # Model time step [s]
     lowerlake_ind::Vector{Int} | "-"                # Index of lower lake (linked lakes)
     area::Vector{T} | "m2"                          # lake area [m²]
     maxstorage::Vector{Union{T, Missing}} | "m3"    # lake maximum storage from rating curve 1 [m³]
@@ -14,7 +13,7 @@
 end
 
 "Initialize lake model parameters"
-function LakeParameters(config, dataset, inds_riv, nriv, pits, dt)
+function LakeParameters(config, dataset, inds_riv, nriv, pits)
     # read only lake data if lakes true
     # allow lakes only in river cells
     # note that these locations are only the lake outlet pixels
@@ -181,7 +180,6 @@ function LakeParameters(config, dataset, inds_riv, nriv, pits, dt)
         end
     end
     parameters = LakeParameters{Float}(;
-        dt = dt,
         lowerlake_ind = lowerlake_ind,
         area = lakearea,
         maxstorage = maximum_storage(lake_storfunc, lake_outflowfunc, lakearea, sh, hq),
@@ -249,9 +247,9 @@ end
 end
 
 "Initialize lake model"
-function Lake(dataset, config, indices_river, n_river_cells, pits, dt)
+function Lake(dataset, config, indices_river, n_river_cells, pits)
     parameters, lake_network, inds_lake_map2river, lake_waterlevel, pits =
-        LakeParameters(dataset, config, indices_river, n_river_cells, pits, dt)
+        LakeParameters(dataset, config, indices_river, n_river_cells, pits)
 
     n_lakes = length(parameters.area)
     variables = LakeVariables(n_lakes, lake_waterlevel)
@@ -324,7 +322,7 @@ Update a single lake at position `i`.
 This is called from within the kinematic wave loop, therefore updating only for a single
 element rather than all at once.
 """
-function update!(model::Lake, i, inflow, doy, timestepsecs)
+function update!(model::Lake, i, inflow, doy, dt, dt_forcing)
     lake_bc = model.boundary_conditions
     lake_p = model.parameters
     lake_v = model.variables
@@ -333,20 +331,19 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
     has_lowerlake = lo != 0
 
     # limit lake evaporation based on total available volume [m³]
-    precipitation =
-        0.001 * lake_bc.precipitation[i] * (timestepsecs / lake_p.dt) * lake_p.area[i]
-    available_volume = lake_v.storage[i] + inflow * timestepsecs + precipitation
-    evap = 0.001 * lake_bc.evaporation[i] * (timestepsecs / lake_p.dt) * lake_p.area[i]
-    actevap = min(available_volume, evap) # [m³/timestepsecs]
+    precipitation = 0.001 * lake_bc.precipitation[i] * (dt / dt_forcing) * lake_p.area[i]
+    available_volume = lake_v.storage[i] + inflow * dt + precipitation
+    evap = 0.001 * lake_bc.evaporation[i] * (dt / dt_forcing) * lake_p.area[i]
+    actevap = min(available_volume, evap) # [m³/dt]
 
     ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ###
     # outflowfunc = 3
     # Calculate lake factor and SI parameter
     if lake_p.outflowfunc[i] == 3
-        lakefactor = lake_p.area[i] / (timestepsecs * pow(lake_p.b[i], 0.5))
-        si_factor = (lake_v.storage[i] + precipitation - actevap) / timestepsecs + inflow
+        lakefactor = lake_p.area[i] / (dt * pow(lake_p.b[i], 0.5))
+        si_factor = (lake_v.storage[i] + precipitation - actevap) / dt + inflow
         # Adjust SIFactor for ResThreshold != 0
-        si_factor_adj = si_factor - lake_p.area[i] * lake_p.threshold[i] / timestepsecs
+        si_factor_adj = si_factor - lake_p.area[i] * lake_p.threshold[i] / dt
         # Calculate the new lake outflow/waterlevel/storage
         if si_factor_adj > 0.0
             quadratic_sol_term =
@@ -360,7 +357,7 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
             outflow = 0.0
         end
         outflow = min(outflow, si_factor)
-        storage = (si_factor - outflow) * timestepsecs
+        storage = (si_factor - outflow) * dt
         waterlevel = storage / lake_p.area[i]
     end
 
@@ -368,8 +365,7 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
     if lake_p.outflowfunc[i] == 1 || lake_p.outflowfunc[i] == 2
         diff_wl = has_lowerlake ? lake_v.waterlevel[i] - lake_v.waterlevel[lo] : 0.0
 
-        storage_input =
-            (lake_v.storage[i] + precipitation - actevap) / timestepsecs + inflow
+        storage_input = (lake_v.storage[i] + precipitation - actevap) / dt + inflow
         if lake_p.outflowfunc[i] == 1
             outflow = interpolate_linear(
                 lake_v.waterlevel[i],
@@ -382,7 +378,7 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
                 if lake_v.waterlevel[i] > lake_p.threshold[i]
                     dh = lake_v.waterlevel[i] - lake_p.threshold[i]
                     outflow = lake_p.b[i] * pow(dh, lake_p.e[i])
-                    maxflow = (dh * lake_p.area[i]) / timestepsecs
+                    maxflow = (dh * lake_p.area[i]) / dt
                     outflow = min(outflow, maxflow)
                 else
                     outflow = Float(0)
@@ -391,18 +387,18 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
                 if lake_v.waterlevel[lo] > lake_p.threshold[i]
                     dh = lake_v.waterlevel[lo] - lake_p.threshold[i]
                     outflow = -1.0 * lake_p.b[i] * pow(dh, lake_p.e[i])
-                    maxflow = (dh * lake_p.area[lo]) / timestepsecs
+                    maxflow = (dh * lake_p.area[lo]) / dt
                     outflow = max(outflow, -maxflow)
                 else
                     outflow = Float(0)
                 end
             end
         end
-        storage = (storage_input - outflow) * timestepsecs
+        storage = (storage_input - outflow) * dt
 
         # update storage and outflow for lake with rating curve of type 1.
         if lake_p.outflowfunc[i] == 1
-            overflow = max(0.0, (storage - lake_p.maxstorage[i]) / timestepsecs)
+            overflow = max(0.0, (storage - lake_p.maxstorage[i]) / dt)
             storage = min(storage, lake_p.maxstorage[i])
             outflow = outflow + overflow
         end
@@ -415,7 +411,7 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
 
         # update lower lake (linked lakes) in case flow from lower lake to upper lake occurs
         if diff_wl < 0.0
-            lowerlake_storage = lake_v.storage[lo] + outflow * timestepsecs
+            lowerlake_storage = lake_v.storage[lo] + outflow * dt
 
             lowerlake_waterlevel = if lake_p.storfunc[lo] == 1
                 lake_v.waterlevel[lo] +
@@ -426,7 +422,7 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
 
             # update values for the lower lake in place
             lake_v.outflow[lo] = -outflow
-            lake_v.totaloutflow[lo] += -outflow * timestepsecs
+            lake_v.totaloutflow[lo] += -outflow * dt
             lake_v.storage[lo] = lowerlake_storage
             lake_v.waterlevel[lo] = lowerlake_waterlevel
         end
@@ -435,8 +431,8 @@ function update!(model::Lake, i, inflow, doy, timestepsecs)
     # update values in place
     lake_v.outflow[i] = max(outflow, 0.0) # for a linked lake flow can be negative
     lake_v.waterlevel[i] = waterlevel
-    lake_bc.inflow[i] += inflow * timestepsecs
-    lake_v.totaloutflow[i] += outflow * timestepsecs
+    lake_bc.inflow[i] += inflow * dt
+    lake_v.totaloutflow[i] += outflow * dt
     lake_v.storage[i] = storage
     lake_v.actevap[i] += 1000.0 * (actevap / lake_p.area[i])
 
diff --git a/src/routing/reservoir.jl b/src/routing/reservoir.jl
index 6fea790f4..33a053249 100644
--- a/src/routing/reservoir.jl
+++ b/src/routing/reservoir.jl
@@ -1,6 +1,5 @@
 "Struct for storing reservoir model parameters"
 @get_units @grid_loc @with_kw struct ReservoirParameters{T}
-    dt::T                                   # Model time step [s]
     maxvolume::Vector{T} | "m3"             # maximum storage (above which water is spilled) [m³]
     area::Vector{T} | "m2"                  # reservoir area [m²]
     maxrelease::Vector{T} | "m3 s-1"        # maximum amount that can be released if below spillway [m³ s⁻¹]
@@ -10,7 +9,7 @@
 end
 
 "Initialize reservoir model parameters"
-function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt)
+function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits)
     # read only reservoir data if reservoirs true
     # allow reservoirs only in river cells
     # note that these locations are only the reservoir outlet pixels
@@ -125,7 +124,6 @@ function ReservoirParameters(dataset, config, indices_river, n_river_cells, pits
     )
 
     parameters = ReservoirParameters{Float}(;
-        dt = dt,
         demand = resdemand,
         maxrelease = resmaxrelease,
         maxvolume = resmaxvolume,
@@ -186,9 +184,9 @@ end
 end
 
 "Initialize reservoir model `SimpleReservoir`"
-function SimpleReservoir(dataset, config, indices_river, n_river_cells, pits, dt)
+function SimpleReservoir(dataset, config, indices_river, n_river_cells, pits)
     parameters, reservoir_network, inds_reservoir_map2river, pits =
-        ReservoirParameters(dataset, config, indices_river, n_river_cells, pits, dt)
+        ReservoirParameters(dataset, config, indices_river, n_river_cells, pits)
 
     n_reservoirs = length(parameters.area)
     @info "Read `$n_reservoirs` reservoir locations."
@@ -206,41 +204,39 @@ Update a single reservoir at position `i`.
 This is called from within the kinematic wave loop, therefore updating only for a single
 element rather than all at once.
 """
-function update!(model::SimpleReservoir, i, inflow, timestepsecs)
+function update!(model::SimpleReservoir, i, inflow, dt, dt_forcing)
     res_bc = model.boundary_conditions
     res_p = model.parameters
     res_v = model.variables
 
     # limit lake evaporation based on total available volume [m³]
-    precipitation =
-        0.001 * res_bc.precipitation[i] * (timestepsecs / res_p.dt) * res_p.area[i]
-    available_volume = res_v.volume[i] + inflow * timestepsecs + precipitation
-    evap = 0.001 * res_bc.evaporation[i] * (timestepsecs / res_p.dt) * res_p.area[i]
-    actevap = min(available_volume, evap) # [m³/timestepsecs]
+    precipitation = 0.001 * res_bc.precipitation[i] * (dt / dt_forcing) * res_p.area[i]
+    available_volume = res_v.volume[i] + inflow * dt + precipitation
+    evap = 0.001 * res_bc.evaporation[i] * (dt / dt_forcing) * res_p.area[i]
+    actevap = min(available_volume, evap) # [m³/dt]
 
-    vol = res_v.volume[i] + (inflow * timestepsecs) + precipitation - actevap
+    vol = res_v.volume[i] + (inflow * dt) + precipitation - actevap
     vol = max(vol, 0.0)
 
     percfull = vol / res_p.maxvolume[i]
     # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level
     fac = scurve(percfull, res_p.targetminfrac[i], Float(1.0), Float(30.0))
-    demandrelease = min(fac * res_p.demand[i] * timestepsecs, vol)
+    demandrelease = min(fac * res_p.demand[i] * dt, vol)
     vol = vol - demandrelease
 
     wantrel = max(0.0, vol - (res_p.maxvolume[i] * res_p.targetfullfrac[i]))
     # Assume extra maximum Q if spilling
     overflow_q = max((vol - res_p.maxvolume[i]), 0.0)
-    torelease =
-        min(wantrel, overflow_q + res_p.maxrelease[i] * timestepsecs - demandrelease)
+    torelease = min(wantrel, overflow_q + res_p.maxrelease[i] * dt - demandrelease)
     vol = vol - torelease
     outflow = torelease + demandrelease
     percfull = vol / res_p.maxvolume[i]
 
     # update values in place
-    res_v.outflow[i] = outflow / timestepsecs
-    res_bc.inflow[i] += inflow * timestepsecs
+    res_v.outflow[i] = outflow / dt
+    res_bc.inflow[i] += inflow * dt
     res_v.totaloutflow[i] += outflow
-    res_v.demandrelease[i] = demandrelease / timestepsecs
+    res_v.demandrelease[i] = demandrelease / dt
     res_v.percfull[i] = percfull
     res_v.volume[i] = vol
     res_v.actevap[i] += 1000.0 * (actevap / res_p.area[i])
diff --git a/src/routing/subsurface.jl b/src/routing/subsurface.jl
index 4a4a116f5..3bf73f642 100644
--- a/src/routing/subsurface.jl
+++ b/src/routing/subsurface.jl
@@ -29,7 +29,6 @@ end
     soilthickness::Vector{T} | "m"         # Soil thickness [m]
     theta_s::Vector{T} | "-"               # Saturated water content (porosity) [-]
     theta_r::Vector{T} | "-"               # Residual water content [-]   
-    dt::T                                  # model time step [d]
     slope::Vector{T} | "m m-1"             # Slope [m m⁻¹]
     flow_length::Vector{T} | "m"           # Flow length [m]
     flow_width::Vector{T} | "m"            # Flow width [m]
@@ -44,7 +43,6 @@ function LateralSsfParameters(
     slope,
     flow_length,
     flow_width,
-    dt,
 )
     khfrac = ncread(
         dataset,
@@ -60,6 +58,7 @@ function LateralSsfParameters(
     soilthickness = soilthickness .* 0.001
 
     kh_profile_type = get(config.input.vertical, "ksat_profile", "exponential")::String
+    dt = Second(config.timestepsecs) / basetimestep
     if kh_profile_type == "exponential"
         (; kv_0, f) = soil.kv_profile
         kh_0 = khfrac .* kv_0 .* 0.001 .* dt
@@ -79,7 +78,6 @@ function LateralSsfParameters(
         soilthickness,
         theta_s,
         theta_r,
-        dt,
         slope,
         flow_length,
         flow_width,
@@ -139,7 +137,6 @@ function LateralSSF(
     flow_width,
     x_length,
     y_length,
-    dt,
 )
     parameters = LateralSsfParameters(
         dataset,
@@ -149,7 +146,6 @@ function LateralSSF(
         slope,
         flow_length,
         flow_width,
-        dt,
     )
     zi = 0.001 * soil.variables.zi
     variables = LateralSsfVariables(parameters, zi, x_length, y_length)
@@ -159,7 +155,7 @@ function LateralSSF(
 end
 
 "Update lateral subsurface model for a single timestep"
-function update!(model::LateralSSF, network)
+function update!(model::LateralSSF, network, dt)
     (;
         order_of_subdomains,
         order_subdomain,
@@ -171,7 +167,7 @@ function update!(model::LateralSSF, network)
 
     (; recharge) = model.boundary_conditions
     (; ssfin, ssf, ssfmax, to_river, zi, exfiltwater, volume) = model.variables
-    (; slope, theta_s, theta_r, soilthickness, flow_length, flow_width, dt, kh_profile) =
+    (; slope, theta_s, theta_r, soilthickness, flow_length, flow_width, kh_profile) =
         model.parameters
 
     ns = length(order_of_subdomains)
@@ -237,23 +233,15 @@ function GroundwaterExchangeVariables(n)
     return variables
 end
 
-"Struct for storing groundwater exchange parameters for coupling with an external groundwater
-model."
-@with_kw struct GroundwaterExchangeParameters{T}
-    dt::T       # model time step [d]
-end
-
 "Groundwater exchange"
 @with_kw struct GroundwaterExchange{T} <: SubsurfaceFlow
-    parameters::GroundwaterExchangeParameters{T}
     variables::GroundwaterExchangeVariables{T}
 end
 
 "Initialize groundwater exchange"
-function GroundwaterExchange(n, dt)
-    parameters = GroundwaterExchangeParameters{Float}(; dt = dt / basetimestep)
+function GroundwaterExchange(n)
     variables = GroundwaterExchangeVariables(n)
-    ssf = GroundwaterExchange{Float}(; parameters, variables)
+    ssf = GroundwaterExchange{Float}(; variables)
     return ssf
 end
 
diff --git a/src/routing/surface_kinwave.jl b/src/routing/surface_kinwave.jl
index 6b1ae9045..94203e9bb 100644
--- a/src/routing/surface_kinwave.jl
+++ b/src/routing/surface_kinwave.jl
@@ -340,7 +340,7 @@ function update!(model::KinWaveOverlandFlow, network, dt)
 end
 
 "Update river flow model `KinWaveRiverFlow` for a single timestep"
-function surfaceflow_river_update!(model::KinWaveRiverFlow, network, doy, dt)
+function surfaceflow_river_update!(model::KinWaveRiverFlow, network, doy, dt, dt_forcing)
     (;
         graph,
         order_of_subdomains,
@@ -390,7 +390,7 @@ function surfaceflow_river_update!(model::KinWaveRiverFlow, network, doy, dt)
                     # run reservoir model and copy reservoir outflow to inflow (qin) of
                     # downstream river cell
                     i = reservoir_indices[v]
-                    update!(reservoir, i, q[v] + inflow_waterbody[v], dt)
+                    update!(reservoir, i, q[v] + inflow_waterbody[v], dt, dt_forcing)
 
                     downstream_nodes = outneighbors(graph, v)
                     n_downstream = length(downstream_nodes)
@@ -411,7 +411,7 @@ function surfaceflow_river_update!(model::KinWaveRiverFlow, network, doy, dt)
                     # run lake model and copy lake outflow to inflow (qin) of downstream river
                     # cell
                     i = lake_indices[v]
-                    update!(lake, i, q[v] + inflow_waterbody[v], doy, dt)
+                    update!(lake, i, q[v] + inflow_waterbody[v], doy, dt, dt_forcing)
 
                     downstream_nodes = outneighbors(graph, v)
                     n_downstream = length(downstream_nodes)
@@ -484,7 +484,7 @@ function update!(model::KinWaveRiverFlow, network, doy, dt)
     while t < dt
         dt_s = adaptive ? stable_timestep(model, 0.05) : model.timestepping.dt_fixed
         dt_s = check_timestepsize(dt_s, t, dt)
-        surfaceflow_river_update!(model, network, doy, dt_s)
+        surfaceflow_river_update!(model, network, doy, dt_s, dt)
         t = t + dt_s
     end
     q_av ./= dt
diff --git a/src/routing/surface_local_inertial.jl b/src/routing/surface_local_inertial.jl
index 55af544f8..33abdeae7 100644
--- a/src/routing/surface_local_inertial.jl
+++ b/src/routing/surface_local_inertial.jl
@@ -303,6 +303,7 @@ function local_inertial_river_update!(
     model::LocalInertialRiverFlow,
     network,
     dt,
+    dt_forcing,
     doy,
     update_h,
 )
@@ -461,7 +462,7 @@ function local_inertial_river_update!(
         i = inds_reservoir[v]
 
         q_in = get_inflow_waterbody(model, edges_at_node.src[i])
-        update!(reservoir, v, q_in + inflow_waterbody[i], dt)
+        update!(reservoir, v, q_in + inflow_waterbody[i], dt, dt_forcing)
         river_v.q[i] = reservoir.variables.outflow[v]
         river_v.q_av[i] += river_v.q[i] * dt
     end
@@ -471,7 +472,7 @@ function local_inertial_river_update!(
         i = inds_lake[v]
 
         q_in = get_inflow_waterbody(model, edges_at_node.src[i])
-        update!(lake, v, q_in + inflow_waterbody[i], doy, dt)
+        update!(lake, v, q_in + inflow_waterbody[i], doy, dt, dt_forcing)
         river_v.q[i] = lake.variables.outflow[v]
         river_v.q_av[i] += river_v.q[i] * dt
     end
@@ -568,7 +569,7 @@ function update!(
         if t + dt_s > dt
             dt_s = dt - t
         end
-        local_inertial_river_update!(model, network, dt_s, doy, update_h)
+        local_inertial_river_update!(model, network, dt_s, dt, doy, update_h)
         t = t + dt_s
     end
     model.variables.q_av ./= dt
@@ -929,7 +930,7 @@ function update!(
         if t + dt_s > dt
             dt_s = dt - t
         end
-        local_inertial_river_update!(river, network, dt_s, doy, update_h)
+        local_inertial_river_update!(river, network, dt_s, dt, doy, update_h)
         local_inertial_update!(land, river, network, dt_s)
         t = t + dt_s
     end
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 62f86e046..31d037e5f 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -100,7 +100,7 @@ function initialize_sbm_gwf_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoir, reservoir_network, inds_reservoir_map2river, pits =
-            SimpleReservoir(dataset, config, inds_river, n_river_cells, pits, tosecond(dt))
+            SimpleReservoir(dataset, config, inds_river, n_river_cells, pits)
     else
         reservoir_network = (river_indices = [],)
         inds_reservoir_map2river = fill(0, n_river_cells)
@@ -110,7 +110,7 @@ function initialize_sbm_gwf_model(config::Config)
     # lakes
     if do_lakes
         lake, lake_network, inds_lake_map2river, pits =
-            Lake(dataset, config, inds_river, n_river_cells, pits, tosecond(dt))
+            Lake(dataset, config, inds_river, n_river_cells, pits)
     else
         lake_network = (river_indices = [],)
         inds_lake_map2river = fill(0, n_river_cells)
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index cd11e45f6..e03d48cbf 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -94,7 +94,7 @@ function initialize_sbm_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoir, reservoir_network, inds_reservoir_map2river, pits =
-            SimpleReservoir(dataset, config, inds_river, n_river_cells, pits, tosecond(dt))
+            SimpleReservoir(dataset, config, inds_river, n_river_cells, pits)
     else
         reservoir_network = (river_indices = [],)
         inds_reservoir_map2river = fill(0, n_river_cells)
@@ -104,7 +104,7 @@ function initialize_sbm_model(config::Config)
     # lakes
     if do_lakes
         lake, lake_network, inds_lake_map2river, pits =
-            Lake(dataset, config, inds_river, n_river_cells, pits, tosecond(dt))
+            Lake(dataset, config, inds_river, n_river_cells, pits)
     else
         lake_network = (river_indices = [],)
         inds_lake_map2river = fill(0, n_river_cells)
@@ -135,7 +135,6 @@ function initialize_sbm_model(config::Config)
     # to another groundwater model, this component is not defined in the TOML file.
     do_lateral_ssf = haskey(config.input.lateral, "subsurface")
     if do_lateral_ssf
-        dt_ssf = dt / basetimestep
         subsurface_flow = LateralSSF(
             dataset,
             config,
@@ -146,7 +145,6 @@ function initialize_sbm_model(config::Config)
             flow_width,
             x_length,
             y_length,
-            dt = dt_ssf,
         )
         # update variables `ssf`, `ssfmax` and `kh` (layered profile) based on ksat_profile
         kh_profile_type = get(config.input.vertical, "ksat_profile", "exponential")::String
@@ -164,7 +162,7 @@ function initialize_sbm_model(config::Config)
     else
         # when the SBM model is coupled (BMI) to a groundwater model, the following
         # variables are expected to be exchanged from the groundwater model.
-        subsurface_flow = GroundwaterExchange(n_land_cells, dt)
+        subsurface_flow = GroundwaterExchange(n_land_cells)
     end
 
     graph = flowgraph(ldd, indices, pcr_dir)
@@ -444,7 +442,7 @@ function update!(model::Model{N, L, V, R, W, T}) where {N, L, V, R, W, T <: SbmM
     lateral.subsurface.variables.zi .= vertical.soil.variables.zi ./ 1000.0
     # update lateral subsurface flow domain (kinematic wave)
     kh_layered_profile!(vertical.soil, lateral.subsurface, kv_profile, dt)
-    update!(lateral.subsurface, network.land)
+    update!(lateral.subsurface, network.land, clock.dt / basetimestep)
     update_after_subsurfaceflow!(model)
     update_total_water_storage!(model)
     return nothing
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index 02d569341..928ee01a7 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -2,7 +2,6 @@
 res_bc =
     Wflow.ReservoirBC{Float}(; inflow = [0.0], precipitation = [4.2], evaporation = [1.5])
 res_params = Wflow.ReservoirParameters{Float64}(;
-    dt = 86400.0,
     demand = [52.523],
     maxrelease = [420.184],
     maxvolume = [25_000_000.0],
@@ -25,7 +24,7 @@ res = Wflow.SimpleReservoir{Float64}(;
     variables = res_vars,
 )
 @testset "Update reservoir simple" begin
-    Wflow.update!(res, 1, 100.0, 86400.0)
+    Wflow.update!(res, 1, 100.0, 86400.0, 86400.0)
     @test res.variables.outflow[1] ≈ 91.3783714867453
     @test res.variables.totaloutflow[1] ≈ 7.895091296454794e6
     @test res.variables.volume[1] ≈ 2.0e7
@@ -38,7 +37,6 @@ end
 
 lake_bc = Wflow.LakeBC{Float}(; inflow = [0.0], precipitation = [20.0], evaporation = [3.2])
 lake_params = Wflow.LakeParameters{Float}(;
-    dt = 86400.0,
     lowerlake_ind = [0],
     area = [180510409.0],
     maxstorage = Wflow.maximum_storage([1], [3], [180510409.0], [missing], [missing]),
@@ -67,7 +65,7 @@ lake = Wflow.Lake{Float64}(;
     lake_p = lake.parameters
     lake_v = lake.variables
     lake_bc = lake.boundary_conditions
-    Wflow.update!(lake, 1, 2500.0, 181, 86400.0)
+    Wflow.update!(lake, 1, 2500.0, 181, 86400.0, 86400.0)
     @test Wflow.waterlevel(lake_p.storfunc, lake_p.area, lake_v.storage, lake_p.sh)[1] ≈
           19.672653848925634
     @test lake_v.outflow[1] ≈ 85.14292808113598
@@ -89,7 +87,6 @@ sh = [
     @test typeof(values(sh[1])) == Tuple{Vector{Float}, Vector{Float}}
 
     lake_params = Wflow.LakeParameters{Float}(;
-        dt = 86400.0,
         lowerlake_ind = [2, 0],
         area = [472461536.0, 60851088.0],
         maxstorage = Wflow.maximum_storage(
@@ -131,8 +128,8 @@ sh = [
         variables = lake_vars,
     )
 
-    Wflow.update!(lake, 1, 500.0, 15, 86400.0)
-    Wflow.update!(lake, 2, 500.0, 15, 86400.0)
+    Wflow.update!(lake, 1, 500.0, 15, 86400.0, 86400.0)
+    Wflow.update!(lake, 2, 500.0, 15, 86400.0, 86400.0)
     lake_v = lake.variables
     lake_bc = lake.boundary_conditions
     @test lake_v.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2
@@ -142,8 +139,8 @@ sh = [
     lake_v.actevap .= 0.0
     lake_v.totaloutflow .= 0.0
     lake_bc.inflow .= 0.0
-    Wflow.update!(lake, 1, 500.0, 15, 86400.0)
-    Wflow.update!(lake, 2, 500.0, 15, 86400.0)
+    Wflow.update!(lake, 1, 500.0, 15, 86400.0, 86400.0)
+    Wflow.update!(lake, 2, 500.0, 15, 86400.0, 86400.0)
     @test lake_v.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2
     @test lake_v.totaloutflow ≈ [-2.2446764487487033e7, 4.3154002238515094e7] atol = 1e3
     @test lake_v.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4
@@ -155,7 +152,6 @@ end
     lake_bc =
         Wflow.LakeBC{Float}(; inflow = [0.00], precipitation = [10.0], evaporation = [2.0])
     lake_params = Wflow.LakeParameters{Float}(;
-        dt = 86400.0,
         lowerlake_ind = [0],
         area = [200_000_000],
         maxstorage = Wflow.maximum_storage(
@@ -186,7 +182,7 @@ end
         variables = lake_vars,
     )
 
-    Wflow.update!(lake, 1, 1500.0, 15, 86400.0)
+    Wflow.update!(lake, 1, 1500.0, 15, 86400.0, 86400.0)
     lake_p = lake.parameters
     lake_v = lake.variables
     @test Wflow.waterlevel(lake_p.storfunc, lake_p.area, lake_v.storage, lake_p.sh) ≈
diff --git a/test/routing_process.jl b/test/routing_process.jl
index 2a8447d4c..931060fbe 100644
--- a/test/routing_process.jl
+++ b/test/routing_process.jl
@@ -244,7 +244,7 @@ end
         sw_river.boundary_conditions.inwater[1] = 20.0
         h0 = mean(sw_river.variables.h)
         dt = Wflow.stable_timestep(sw_river)
-        Wflow.local_inertial_river_update!(sw_river, network, dt, 0.0, true)
+        Wflow.local_inertial_river_update!(sw_river, network, dt, 86400.0, 0.0, true)
         d = abs(h0 - mean(sw_river.variables.h))
         if d <= epsilon
             break