Units reservoir and lake, rename totaloutflow to outflow_av

Change unit (m3 per model timestep) `totaloutflow` (renamed to `outflow_av`) and `inflow` to m3/s . Use functions for setting reservoir or lake variables (`inflow`, `outflow_av` and `actevap`) at the start of each timestep to zero and to compute averages valid for the model timestep for variables `inflow` and `outflow_av`.

src/Wflow.jl

@@ -146,6 +146,8 @@ include("groundwater/aquifer.jl")
 include("groundwater/boundary_conditions.jl")
 include("routing/timestepping.jl")
 include("routing/subsurface.jl")
+include("routing/reservoir.jl")
+include("routing/lake.jl")
 include("routing/surface_kinwave.jl")
 include("routing/surface_local_inertial.jl")
 include("routing/surface_routing.jl")
@@ -162,8 +164,6 @@ include("soil/soil_process.jl")
 include("sbm.jl")
 include("demand/water_demand.jl")
 include("sediment.jl")
-include("routing/reservoir.jl")
-include("routing/lake.jl")
 include("sbm_model.jl")
 include("sediment_model.jl")
 include("sbm_gwf_model.jl")

src/routing/lake.jl

@@ -204,8 +204,8 @@ end
 @get_units @grid_loc @with_kw struct LakeVariables{T}
     waterlevel::Vector{T} | "m"                 # waterlevel H [m] of lake
     storage::Vector{T} | "m3"                   # storage lake [m³]
-    outflow::Vector{T} | "m3 s-1"               # outflow lake [m³ s⁻¹]
-    totaloutflow::Vector{T} | "m3"              # total outflow lake [m³]
+    outflow::Vector{T} | "m3 s-1"               # outflow of lake outlet [m³ s⁻¹]
+    outflow_av::Vector{T} | "m3 s-1"            # average outflow lake [m³ s⁻¹] for model timestep Δt (including flow from lower to upper lake)
     actevap::Vector{T}                          # average actual evapotranspiration for lake area [mm Δt⁻¹] 
 end
 
@@ -216,15 +216,15 @@ function LakeVariables(n, lake_waterlevel)
         inflow = fill(mv, n),
         storage = initialize_storage(lake_storfunc, lakearea, lake_waterlevel, sh),
         outflow = fill(mv, n),
-        totaloutflow = fill(mv, n),
+        outflow_av = fill(mv, n),
         actevap = fill(mv, n),
     )
     return variables
 end
 
 "Struct for storing lake model boundary conditions"
 @get_units @grid_loc @with_kw struct LakeBC{T}
-    inflow::Vector{T} | "m3"                    # inflow to the lake [m³]
+    inflow::Vector{T} | "m3 s-1"                # inflow to the lake [m³ s⁻¹] for model timestep Δt
     precipitation::Vector{T}                    # average precipitation for lake area [mm Δt⁻¹]
     evaporation::Vector{T}                      # average potential evaporation for lake area [mm Δt⁻¹]
 end
@@ -422,17 +422,17 @@ function update!(model::Lake, i, inflow, doy, dt, dt_forcing)
 
             # update values for the lower lake in place
             lake_v.outflow[lo] = -outflow
-            lake_v.totaloutflow[lo] += -outflow * dt
+            lake_v.outflow_av[lo] += -outflow * dt
             lake_v.storage[lo] = lowerlake_storage
             lake_v.waterlevel[lo] = lowerlake_waterlevel
         end
     end
 
     # update values in place
-    lake_v.outflow[i] = max(outflow, 0.0) # for a linked lake flow can be negative
+    lake_v.outflow[i] = outflow
     lake_v.waterlevel[i] = waterlevel
     lake_bc.inflow[i] += inflow * dt
-    lake_v.totaloutflow[i] += outflow * dt
+    lake_v.outflow_av[i] += outflow * dt
     lake_v.storage[i] = storage
     lake_v.actevap[i] += 1000.0 * (actevap / lake_p.area[i])
 

src/routing/reservoir.jl

@@ -139,7 +139,7 @@ end
 @get_units @grid_loc @with_kw struct ReservoirVariables{T}
     volume::Vector{T} | "m3"                # reservoir volume [m³]
     outflow::Vector{T} | "m3 s-1"           # outflow from reservoir [m³ s⁻¹]
-    totaloutflow::Vector{T} | "m3"          # total outflow from reservoir [m³]
+    outflow_av::Vector{T} | "m3 s-1"        # average outflow from reservoir [m³ s⁻¹] for model timestep Δt
     percfull::Vector{T} | "-"               # fraction full (of max storage) [-]
     demandrelease::Vector{T} | "m3 s-1"     # minimum (environmental) flow released from reservoir [m³ s⁻¹]
     actevap::Vector{T}                      # average actual evaporation for reservoir area [mm Δt⁻¹]
@@ -151,7 +151,7 @@ function ReservoirVariables(n, parameters)
     variables = ReservoirVariables{Float}(;
         volume = targetfullfrac .* maxvolume,
         outflow = fill(mv, n),
-        totaloutflow = fill(mv, n),
+        outflow_av = fill(mv, n),
         percfull = fill(mv, n),
         demandrelease = fill(mv, n),
         actevap = fill(mv, n),
@@ -161,7 +161,7 @@ end
 
 "Struct for storing reservoir model boundary conditions"
 @get_units @grid_loc @with_kw struct ReservoirBC{T}
-    inflow::Vector{T} | "m3"        # total inflow into reservoir [m³]
+    inflow::Vector{T} | "m3 s-1"    # inflow into reservoir [m³ s⁻¹] for model timestep Δt
     precipitation::Vector{T}        # average precipitation for reservoir area [mm Δt⁻¹]
     evaporation::Vector{T}          # average potential evaporation for reservoir area [mm Δt⁻¹]
 end
@@ -235,7 +235,7 @@ function update!(model::SimpleReservoir, i, inflow, dt, dt_forcing)
     # update values in place
     res_v.outflow[i] = outflow / dt
     res_bc.inflow[i] += inflow * dt
-    res_v.totaloutflow[i] += outflow
+    res_v.outflow_av[i] += outflow
     res_v.demandrelease[i] = demandrelease / dt
     res_v.percfull[i] = percfull
     res_v.volume[i] = vol

src/routing/surface_kinwave.jl

@@ -246,6 +246,29 @@ function KinWaveOverlandFlow(dataset, config, indices; slope, flow_length, flow_
     return sf_land
 end
 
+"""
+Helper function to set waterbody variables `inflow`,`outflow_av` and `actevap` to zero. This
+is done at the start of each simulation timestep, during the timestep the total (weighted)
+sum is computed from values at each sub timestep.
+"""
+function set_waterbody_vars!(waterbody::W) where {W <: Union{SimpleReservoir, Lake}}
+    waterbody.boundary_conditions.inflow .= 0.0
+    waterbody.variables.outflow_av .= 0.0
+    waterbody.variables.actevap .= 0.0
+    return nothing
+end
+set_waterbody_vars!(waterbody) = nothing
+
+"""
+Helper function to compute the average of waterbody variables `inflow` and `outflow_av`. This
+is done at the end of each simulation timestep.
+"""
+function average_waterbody_vars!(waterbody::W, dt) where {W <: Union{SimpleReservoir, Lake}}
+    waterbody.variables.outflow_av ./= dt
+    waterbody.boundary_conditions.inflow ./= dt
+end
+average_waterbody_vars!(waterbody, dt) = nothing
+
 "Update overland flow model `KinWaveOverlandFlow` for a single timestep"
 function kinwave_land_update!(model::KinWaveOverlandFlow, network, dt)
     (;
@@ -417,7 +440,7 @@ function kinwave_river_update!(model::KinWaveRiverFlow, network, doy, dt, dt_for
                     n_downstream = length(downstream_nodes)
                     if n_downstream == 1
                         j = only(downstream_nodes)
-                        qin[j] = lake.variables.outflow[i]
+                        qin[j] = max(lake.variables.outflow[i], 0.0)
                     elseif n_downstream == 0
                         error(
                             """A lake without a downstream river node is not supported.
@@ -467,18 +490,9 @@ function update!(model::KinWaveRiverFlow, network, doy, dt)
 
     q_av .= 0.0
     h_av .= 0.0
-    # because of possible iterations set reservoir and lake inflow and total outflow at
-    # start to zero, the total sum of inflow and outflow at each sub time step is calculated
-    if !isnothing(reservoir)
-        reservoir.boundary_conditions.inflow .= 0.0
-        reservoir.variables.totaloutflow .= 0.0
-        reservoir.variables.actevap .= 0.0
-    end
-    if !isnothing(lake)
-        lake.boundary_conditions.inflow .= 0.0
-        lake.variables.totaloutflow .= 0.0
-        lake.variables.actevap .= 0.0
-    end
+
+    set_waterbody_vars!(reservoir)
+    set_waterbody_vars!(lake)
 
     t = 0.0
     while t < dt
@@ -487,6 +501,10 @@ function update!(model::KinWaveRiverFlow, network, doy, dt)
         kinwave_river_update!(model, network, doy, dt_s, dt)
         t = t + dt_s
     end
+
+    average_waterbody_vars!(reservoir, dt)
+    average_waterbody_vars!(lake, dt)
+
     q_av ./= dt
     h_av ./= dt
     volume .= flow_length .* flow_width .* h

src/routing/surface_local_inertial.jl

@@ -473,7 +473,7 @@ function local_inertial_river_update!(
 
         q_in = get_inflow_waterbody(model, edges_at_node.src[i])
         update!(lake, v, q_in + inflow_waterbody[i], doy, dt, dt_forcing)
-        river_v.q[i] = lake.variables.outflow[v]
+        river_v.q[i] = max(lake.variables.outflow[v], 0.0)
         river_v.q_av[i] += river_v.q[i] * dt
     end
     if update_h
@@ -546,16 +546,10 @@ function update!(
     update_h = true,
 ) where {T}
     (; reservoir, lake) = model.boundary_conditions
-    if !isnothing(reservoir)
-        reservoir.boundary_conditions.inflow .= 0.0
-        reservoir.variables.totaloutflow .= 0.0
-        reservoir.variables.actevap .= 0.0
-    end
-    if !isnothing(lake)
-        lake.boundary_conditions.inflow .= 0.0
-        lake.variables.totaloutflow .= 0.0
-        lake.variables.actevap .= 0.0
-    end
+
+    set_waterbody_vars!(reservoir)
+    set_waterbody_vars!(lake)
+
     if !isnothing(model.floodplain)
         model.floodplain.variables.q_av .= 0.0
         model.floodplain.variables.h_av .= 0.0
@@ -575,6 +569,9 @@ function update!(
     model.variables.q_av ./= dt
     model.variables.h_av ./= dt
 
+    average_waterbody_vars!(reservoir, dt)
+    average_waterbody_vars!(lake, dt)
+
     if !isnothing(model.floodplain)
         model.floodplain.variables.q_av ./= dt
         model.floodplain.variables.h_av ./= dt

test/reservoir_lake.jl

@@ -11,7 +11,7 @@ res_params = Wflow.ReservoirParameters{Float64}(;
 )
 res_vars = Wflow.ReservoirVariables{Float64}(;
     volume = [1.925e7],
-    totaloutflow = [0.0],
+    outflow_av = [0.0],
     actevap = [0.0],
     outflow = [NaN],
     percfull = [NaN],
@@ -25,8 +25,9 @@ res = Wflow.SimpleReservoir{Float64}(;
 )
 @testset "Update reservoir simple" begin
     Wflow.update!(res, 1, 100.0, 86400.0, 86400.0)
+    res.variables.outflow_av ./= 86400.0
     @test res.variables.outflow[1] ≈ 91.3783714867453
-    @test res.variables.totaloutflow[1] ≈ 7.895091296454794e6
+    @test res.variables.outflow_av[1] == res.variables.outflow[1]
     @test res.variables.volume[1] ≈ 2.0e7
     @test res.variables.percfull[1] ≈ 0.80
     @test res.variables.demandrelease[1] ≈ 52.5229994727611
@@ -49,7 +50,7 @@ lake_params = Wflow.LakeParameters{Float}(;
     hq = [missing],
 )
 lake_vars = Wflow.LakeVariables{Float}(;
-    totaloutflow = [0.0],
+    outflow_av = [0.0],
     storage = Wflow.initialize_storage([1], [180510409.0], [18.5], [missing]),
     waterlevel = [18.5],
     actevap = [0.0],
@@ -66,10 +67,11 @@ lake = Wflow.Lake{Float64}(;
     lake_v = lake.variables
     lake_bc = lake.boundary_conditions
     Wflow.update!(lake, 1, 2500.0, 181, 86400.0, 86400.0)
+    lake_v.outflow_av ./= 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
-    @test lake_v.totaloutflow[1] ≈ 7.356348986210149e6
+    @test lake_v.outflow_av[1] ≈ lake_v.outflow[1]
     @test lake_v.storage[1] ≈ 3.55111879238499e9
     @test lake_v.waterlevel[1] ≈ 19.672653848925634
     @test lake_bc.precipitation[1] ≈ 20.0
@@ -105,7 +107,7 @@ sh = [
         hq = [missing, Wflow.read_hq_csv(joinpath(datadir, "input", "lake_hq_2.csv"))],
     )
     lake_vars = Wflow.LakeVariables{Float}(;
-        totaloutflow = [0.0, 0.0],
+        outflow_av = [0.0, 0.0],
         waterlevel = [395.03027, 394.87833],
         actevap = [0.0, 0.0],
         outflow = [NaN, NaN],
@@ -130,21 +132,22 @@ sh = [
 
     Wflow.update!(lake, 1, 500.0, 15, 86400.0, 86400.0)
     Wflow.update!(lake, 2, 500.0, 15, 86400.0, 86400.0)
+    lake.variables.outflow_av ./= 86400.0
     lake_v = lake.variables
     lake_bc = lake.boundary_conditions
-    @test lake_v.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2
-    @test lake_v.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3
-    @test lake_v.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4
-    @test lake_v.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2
+    @test lake_v.outflow ≈ [214.80170846121263, 236.83281600000214]
+    @test lake_v.outflow_av ≈ lake_v.outflow
+    @test lake_v.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8]
     lake_v.actevap .= 0.0
-    lake_v.totaloutflow .= 0.0
+    lake_v.outflow_av .= 0.0
     lake_bc.inflow .= 0.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
-    @test lake_v.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2
+    lake.variables.outflow_av ./= 86400.0
+    @test lake_v.outflow ≈ [-259.8005149014703, 239.66710359986183]
+    @test lake_v.outflow_av ≈ [-259.8005149014703, 499.4676185013321]
+    @test lake_v.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8]
+    @test lake_v.waterlevel ≈ [395.239782021054, 395.21771942667266]
     @test lake_v.actevap ≈ [2.0, 2.0]
 end
 
@@ -170,7 +173,7 @@ end
         hq = [Wflow.read_hq_csv(joinpath(datadir, "input", "lake_hq_2.csv"))],
     )
     lake_vars = Wflow.LakeVariables{Float}(;
-        totaloutflow = [0.0],
+        outflow_av = [0.0],
         waterlevel = [397.75],
         actevap = [0.0],
         outflow = [NaN],
@@ -183,12 +186,13 @@ end
     )
 
     Wflow.update!(lake, 1, 1500.0, 15, 86400.0, 86400.0)
+    lake.variables.outflow_av ./= 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) ≈
           [398.0] atol = 1e-2
-    @test lake_v.outflow ≈ [1303.67476852] atol = 1e-2
-    @test lake_v.totaloutflow ≈ [11.26375000e7] atol = 1e3
-    @test lake_v.storage ≈ [4.293225e8] atol = 1e4
-    @test lake_v.waterlevel ≈ [398.000000] atol = 1e-2
+    @test lake_v.outflow ≈ [1303.67476852]
+    @test lake_v.outflow_av ≈ lake_v.outflow
+    @test lake_v.storage ≈ [4.293225e8]
+    @test lake_v.waterlevel ≈ [398.000000]
 end

test/run_sbm.jl

@@ -133,7 +133,7 @@ end
 @testset "reservoir simple" begin
     res = model.lateral.river.boundary_conditions.reservoir
     @test res.variables.outflow[1] ≈ 0.21750000119148086f0
-    @test res.boundary_conditions.inflow[1] ≈ 44.312772458823474f0
+    @test res.boundary_conditions.inflow[1] ≈ 0.00051287944327482
     @test res.variables.volume[1] ≈ 2.751299001489657f7
     @test res.boundary_conditions.precipitation[1] ≈ 0.17999997735023499f0
     @test res.boundary_conditions.evaporation[1] ≈ 0.5400000810623169f0