Remove model parameter dt

From structs `LakeParameters`, `ReservoirParameters`, `LateralSsfParameters` and `GroundwaterExchangeParameters`.

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,16 +357,15 @@ 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
 
     ### Linearisation for specific storage/rating curves ###
     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])
 

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])

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
 

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