From 8c98969ea5012fb1e6122341b5bb2de503e02501 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Thu, 18 Jan 2024 09:54:07 +0100
Subject: [PATCH 1/9] Limit lake evaporation Based on total available volume
 (storage, inflow and precipitation).

---
 src/flow.jl            |  3 +++
 src/reservoir_lake.jl  | 27 +++++++++++++--------------
 test/reservoir_lake.jl |  3 +++
 3 files changed, 19 insertions(+), 14 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index af77daa89..056fce145 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -270,6 +270,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy)
     if !isnothing(sf.lake)
         sf.lake.inflow .= 0.0
         sf.lake.totaloutflow .= 0.0
+        sf.lake.actevap .= 0.0
     end
 
     Δt, its = stable_timestep(sf)
@@ -902,6 +903,7 @@ function update(
     if !isnothing(sw.lake)
         sw.lake.inflow .= 0.0
         sw.lake.totaloutflow .= 0.0
+        sw.lake.actevap .= 0.0
     end
     if !isnothing(sw.floodplain)
         sw.floodplain.q_av .= 0.0
@@ -1142,6 +1144,7 @@ function update(
     if !isnothing(swr.lake)
         swr.lake.inflow .= 0.0
         swr.lake.totaloutflow .= 0.0
+        swr.lake.actevap .= 0.0
     end
     swr.q_av .= 0.0
     swr.h_av .= 0.0
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 5abd0506a..2599bbb68 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -222,6 +222,7 @@ end
     totaloutflow::Vector{T} | "m3"              # total outflow lake [m³]
     precipitation::Vector{T}                    # average precipitation for lake area [mm Δt⁻¹]
     evaporation::Vector{T}                      # average evaporation for lake area [mm Δt⁻¹]
+    actevap::Vector{T}                          # average actual evapotranspiration for lake area [mm Δt⁻¹]
 
     function Lake{T}(args...) where {T}
         equal_size_vectors(args)
@@ -415,6 +416,7 @@ function initialize_lake(config, nc, inds_riv, nriv, pits, Δt)
         totaloutflow = fill(mv, n),
         precipitation = fill(mv, n),
         evaporation = fill(mv, n),
+        actevap = fill(mv, n),
     )
 
     return lakes,
@@ -495,22 +497,21 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
     lo = lake.lowerlake_ind[i]
     has_lowerlake = lo != 0
 
+    # limit lake evaporation based on total available volume [m³]
+    precipitation = 0.001 * lake.precipitation[i] * (timestepsecs / lake.Δt) * lake.area[i]
+    available_volume = lake.storage[i] + inflow * timestepsecs + precipitation
+    evap = 0.001 * lake.evaporation[i] * (timestepsecs / lake.Δt) * lake.area[i]
+    actevap = min(available_volume, evap) # [m³/timestepsecs]
+
     ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ###
     # outflowfunc = 3
     # Calculate lake factor and SI parameter
     if lake.outflowfunc[i] == 3
         lakefactor = lake.area[i] / (timestepsecs * pow(lake.b[i], 0.5))
-        si_factor =
-            (
-                lake.storage[i] +
-                (lake.precipitation[i] - lake.evaporation[i]) *
-                (timestepsecs / lake.Δt) *
-                lake.area[i] / 1000.0
-            ) / timestepsecs + inflow
-        #Adjust SIFactor for ResThreshold != 0
+        si_factor = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow
+        # Adjust SIFactor for ResThreshold != 0
         si_factor_adj = si_factor - lake.area[i] * lake.threshold[i] / timestepsecs
-        #Calculate the new lake outflow/waterlevel/storage
-
+        # Calculate the new lake outflow/waterlevel/storage
         if si_factor_adj > 0.0
             outflow = pow(
                 0.5 *
@@ -553,10 +554,7 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
         end
 
         storage =
-            lake.storage[i] +
-            inflow * timestepsecs +
-            (lake.precipitation[i] / 1000.0) * (timestepsecs / lake.Δt) * lake.area[i] -
-            (lake.evaporation[i] / 1000.0) * (timestepsecs / lake.Δt) * lake.area[i] -
+            lake.storage[i] + inflow * timestepsecs + precipitation - actevap -
             outflow * timestepsecs
 
         # update storage and outflow for lake with rating curve of type 1.
@@ -596,6 +594,7 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
     lake.inflow[i] += inflow * timestepsecs
     lake.totaloutflow[i] += outflow * timestepsecs
     lake.storage[i] = storage
+    lake.actevap[i] += 1000.0 * (actevap / lake.area[i])
 
     return lake
 end
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index 66d251fe6..6094e6478 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -46,6 +46,7 @@ end
         waterlevel = [18.5],
         precipitation = [20.0],
         evaporation = [3.2],
+        actevap = [3.2],
         outflow = [NaN],
     )
 
@@ -90,6 +91,7 @@ sh = [
         waterlevel = [395.03027, 394.87833],
         precipitation = [10.0, 10.0],
         evaporation = [2.0, 2.0],
+        actevap = [2.0, 2.0],
         outflow = [NaN, NaN],
         storage = Wflow.initialize_storage(
             [2, 2],
@@ -139,6 +141,7 @@ end
         waterlevel = [397.75],
         precipitation = [10.0],
         evaporation = [2.0],
+        actevap = [2.0],
         outflow = [NaN],
         storage = [410_760_000],
     )

From b12570786f86e632e743b26f204c7b16ed0635da Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Thu, 18 Jan 2024 11:24:17 +0100
Subject: [PATCH 2/9] Limit reservoir evaporation Based on total available
 volume (storage, inflow and precipitation).

---
 src/flow.jl            |  3 +++
 src/reservoir_lake.jl  | 19 ++++++++++---------
 test/bmi.jl            |  4 ++--
 test/reservoir_lake.jl |  1 +
 4 files changed, 16 insertions(+), 11 deletions(-)

diff --git a/src/flow.jl b/src/flow.jl
index 056fce145..ab583963a 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -266,6 +266,7 @@ function update(sf::SurfaceFlowRiver, network, inflow_wb, doy)
     if !isnothing(sf.reservoir)
         sf.reservoir.inflow .= 0.0
         sf.reservoir.totaloutflow .= 0.0
+        sf.reservoir.actevap .= 0.0
     end
     if !isnothing(sf.lake)
         sf.lake.inflow .= 0.0
@@ -899,6 +900,7 @@ function update(
     if !isnothing(sw.reservoir)
         sw.reservoir.inflow .= 0.0
         sw.reservoir.totaloutflow .= 0.0
+        sw.reservoir.actevap .= 0.0
     end
     if !isnothing(sw.lake)
         sw.lake.inflow .= 0.0
@@ -1140,6 +1142,7 @@ function update(
     if !isnothing(swr.reservoir)
         swr.reservoir.inflow .= 0.0
         swr.reservoir.totaloutflow .= 0.0
+        swr.reservoir.actevap .= 0.0
     end
     if !isnothing(swr.lake)
         swr.lake.inflow .= 0.0
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 2599bbb68..a1a3e3147 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -14,6 +14,7 @@
     demandrelease::Vector{T} | "m3 s-1"                 # minimum (environmental) flow released from reservoir [m³ s⁻¹]
     precipitation::Vector{T}                            # average precipitation for reservoir area [mm Δt⁻¹]
     evaporation::Vector{T}                              # average evaporation for reservoir area [mm Δt⁻¹]
+    actevap::Vector{T}                                  # average actual evaporation for reservoir area [mm Δt⁻¹]
 
     function SimpleReservoir{T}(args...) where {T}
         equal_size_vectors(args)
@@ -148,6 +149,7 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, Δt)
         demandrelease = fill(mv, n),
         precipitation = fill(mv, n),
         evaporation = fill(mv, n),
+        actevap = fill(mv, n),
     )
 
     return reservoirs,
@@ -168,15 +170,13 @@ element rather than all at once.
 """
 function update(res::SimpleReservoir, i, inflow, timestepsecs)
 
-    vol = max(
-        0.0,
-        (
-            res.volume[i] +
-            (inflow * timestepsecs) +
-            (res.precipitation[i] * (timestepsecs / res.Δt) / 1000.0) * res.area[i] -
-            (res.evaporation[i] * (timestepsecs / res.Δt) / 1000.0) * res.area[i]
-        ),
-    )
+    # limit lake evaporation based on total available volume [m³]
+    precipitation = 0.001 * res.precipitation[i] * (timestepsecs / res.Δt) * res.area[i]
+    available_volume = res.volume[i] + inflow * timestepsecs + precipitation
+    evap = 0.001 * res.evaporation[i] * (timestepsecs / res.Δt) * res.area[i]
+    actevap = min(available_volume, evap) # [m³/timestepsecs]
+
+    vol = res.volume[i] + (inflow * timestepsecs) + precipitation - actevap
 
     percfull = vol / res.maxvolume[i]
     # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level
@@ -199,6 +199,7 @@ function update(res::SimpleReservoir, i, inflow, timestepsecs)
     res.demandrelease[i] = demandrelease / timestepsecs
     res.percfull[i] = percfull
     res.volume[i] = vol
+    res.actevap[i] += 1000.0 * (actevap / res.area[i])
 
     return res
 end
diff --git a/test/bmi.jl b/test/bmi.jl
index 35d147914..991161be9 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -20,8 +20,8 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "model information functions" begin
             @test BMI.get_component_name(model) == "sbm"
-            @test BMI.get_input_item_count(model) == 183
-            @test BMI.get_output_item_count(model) == 183
+            @test BMI.get_input_item_count(model) == 184
+            @test BMI.get_output_item_count(model) == 184
             @test BMI.get_input_var_names(model)[[1, 5, 150, 174]] == [
                 "vertical.nlayers",
                 "vertical.θᵣ",
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index 6094e6478..b26b74247 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -12,6 +12,7 @@
         targetminfrac = [0.2425554726620697],
         precipitation = [4.2],
         evaporation = [1.5],
+        actevap = [1.5],
         outflow = [NaN],
         percfull = [NaN],
         demandrelease = [NaN],

From 38b06a8da03857df169c7dd25c48aafb8c8a6c46 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Thu, 18 Jan 2024 11:37:53 +0100
Subject: [PATCH 3/9] Update tests (lake and reservoir)

---
 test/reservoir_lake.jl | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index b26b74247..0354f400d 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -12,7 +12,7 @@
         targetminfrac = [0.2425554726620697],
         precipitation = [4.2],
         evaporation = [1.5],
-        actevap = [1.5],
+        actevap = [0.0],
         outflow = [NaN],
         percfull = [NaN],
         demandrelease = [NaN],
@@ -26,6 +26,7 @@
     @test res.demandrelease[1] ≈ 52.5229994727611
     @test res.precipitation[1] ≈ 4.2
     @test res.evaporation[1] ≈ 1.5
+    @test res.actevap[1] ≈ 1.5
 end
 
 @testset "lake" begin
@@ -47,7 +48,7 @@ end
         waterlevel = [18.5],
         precipitation = [20.0],
         evaporation = [3.2],
-        actevap = [3.2],
+        actevap = [0.0],
         outflow = [NaN],
     )
 
@@ -58,6 +59,7 @@ end
     @test lake.waterlevel[1] ≈ 19.672653848925634
     @test lake.precipitation[1] ≈ 20.0
     @test lake.evaporation[1] ≈ 3.2
+    @test lake.actevap[1] ≈ 3.2
 end
 
 datadir = joinpath(@__DIR__, "data")
@@ -92,7 +94,7 @@ sh = [
         waterlevel = [395.03027, 394.87833],
         precipitation = [10.0, 10.0],
         evaporation = [2.0, 2.0],
-        actevap = [2.0, 2.0],
+        actevap = [0.0, 0.0],
         outflow = [NaN, NaN],
         storage = Wflow.initialize_storage(
             [2, 2],
@@ -108,6 +110,7 @@ sh = [
     @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3
     @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4
     @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2
+    lake.actevap .= 0.0
     lake.totaloutflow .= 0.0
     lake.inflow .= 0.0
     Wflow.update(lake, 1, 500.0, 15, 86400.0)
@@ -116,6 +119,7 @@ sh = [
     @test lake.totaloutflow ≈ [-2.2446764487487033e7, 2.070723775102806e7] atol = 1e3
     @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4
     @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2
+    @test lake.actevap ≈ [2.0, 2.0]
 end
 
 @testset "overflowing lake with sh and hq" begin
@@ -142,7 +146,7 @@ end
         waterlevel = [397.75],
         precipitation = [10.0],
         evaporation = [2.0],
-        actevap = [2.0],
+        actevap = [0.0],
         outflow = [NaN],
         storage = [410_760_000],
     )

From 2c154514885392099dfdd47ab3307e6174fffcaa Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Thu, 18 Jan 2024 21:11:33 +0100
Subject: [PATCH 4/9] Limit lake outflow and solution Modified Puls approach
 Limit lake outflow by lake storage, and check if solution term of quadratic
 equations is above zero.

---
 src/reservoir_lake.jl  | 26 +++++++++++++++++---------
 test/reservoir_lake.jl |  2 +-
 2 files changed, 18 insertions(+), 10 deletions(-)

diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index a1a3e3147..ee48c7fa8 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -177,6 +177,7 @@ function update(res::SimpleReservoir, i, inflow, timestepsecs)
     actevap = min(available_volume, evap) # [m³/timestepsecs]
 
     vol = res.volume[i] + (inflow * timestepsecs) + precipitation - actevap
+    vol = max(vol, 0.0)
 
     percfull = vol / res.maxvolume[i]
     # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level
@@ -514,14 +515,17 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
         si_factor_adj = si_factor - lake.area[i] * lake.threshold[i] / timestepsecs
         # Calculate the new lake outflow/waterlevel/storage
         if si_factor_adj > 0.0
-            outflow = pow(
-                0.5 *
-                (-lakefactor + pow((pow(lakefactor, 2.0) + 4.0 * si_factor_adj), 0.5)),
-                2.0,
-            )
+            quadratic_sol_term =
+                -lakefactor + pow((pow(lakefactor, 2.0) + 4.0 * si_factor_adj), 0.5)
+            if quadratic_sol_term > 0.0
+                outflow = pow(0.5 * quadratic_sol_term, 2.0)
+            else
+                outflow = 0.0
+            end
         else
             outflow = 0.0
         end
+        outflow = min(outflow, si_factor)
         storage = (si_factor - outflow) * timestepsecs
         waterlevel = storage / lake.area[i]
     end
@@ -553,10 +557,13 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
                 end
             end
         end
-
-        storage =
-            lake.storage[i] + inflow * timestepsecs + precipitation - actevap -
-            outflow * timestepsecs
+        storage_input = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow
+        if diff_wl >= 0.0
+            outflow = min(outflow, storage_input)
+        else
+            outflow = max(outflow, -lake.storage[lo])
+        end
+        storage = (storage_input - outflow) * timestepsecs
 
         # update storage and outflow for lake with rating curve of type 1.
         if lake.outflowfunc[i] == 1
@@ -583,6 +590,7 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
 
             # update values for the lower lake in place
             lake.outflow[lo] = -outflow
+            lake.totaloutflow[lo] += -outflow * timestepsecs
             lake.storage[lo] = lowerlake_storage
             lake.waterlevel[lo] = lowerlake_waterlevel
         end
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index 0354f400d..6e17dfe77 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -116,7 +116,7 @@ sh = [
     Wflow.update(lake, 1, 500.0, 15, 86400.0)
     Wflow.update(lake, 2, 500.0, 15, 86400.0)
     @test lake.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2
-    @test lake.totaloutflow ≈ [-2.2446764487487033e7, 2.070723775102806e7] atol = 1e3
+    @test lake.totaloutflow ≈ [-2.2446764487487033e7, 4.3154002238515094e7] atol = 1e3
     @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4
     @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2
     @test lake.actevap ≈ [2.0, 2.0]

From 25b6f331f81838a18e1b7585dae9ed34c7cef233 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Fri, 19 Jan 2024 08:17:04 +0100
Subject: [PATCH 5/9] Add parameter `actevap` to docs Lake and reservoir
 module.

---
 docs/src/model_docs/params_lateral.md | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 50326452e..236502508 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -106,6 +106,7 @@ locs = "wflow_reservoirlocs"
 | `percfull`             | fraction full (of max storage) | - | - |
 | `precipitation`             | average precipitation for reservoir area | mm Δt⁻¹ | - |
 | `evaporation`             | average evaporation for reservoir area | mm Δt⁻¹ | - |
+| `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
 
 ### [Lakes](@id lake_params)
 The Table below shows the parameters (fields) of struct `Lake`, including a description of
@@ -147,7 +148,8 @@ between parentheses.
 | `outflow` | outflow lake | m``^3`` s``^{-1}``  | - |
 | `totaloutflow` | total outflow lake | m``^3``  | - |
 | `precipitation` | average precipitation for lake area | mm Δt⁻¹  | - |
-| `evaporation` | average precipitation for lake area | mm Δt⁻¹  | - |
+| `evaporation` | average evaporation for lake area | mm Δt⁻¹  | - |
+| `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
 
 ### [Lateral subsurface flow](@id params_ssf)
 The Table below shows the parameters (fields) of struct `LateralSSF`, including a

From 02ea9ce19b3be7ef4cd147d6003b976557ff141c Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Fri, 19 Jan 2024 08:43:41 +0100
Subject: [PATCH 6/9] Small refactoring (limit lake outflow)

---
 src/reservoir_lake.jl | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index ee48c7fa8..5e0a2d0d5 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -535,14 +535,17 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
 
         diff_wl = has_lowerlake ? lake.waterlevel[i] - lake.waterlevel[lo] : 0.0
 
+        storage_input = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow
         if lake.outflowfunc[i] == 1
             outflow =
                 interpolate_linear(lake.waterlevel[i], lake.hq[i].H, lake.hq[i].Q[:, doy])
+            outflow = min(outflow, storage_input)
         else
             if diff_wl >= 0.0
                 if lake.waterlevel[i] > lake.threshold[i]
                     outflow =
                         lake.b[i] * pow((lake.waterlevel[i] - lake.threshold[i]), lake.e[i])
+                    outflow = min(outflow, storage_input)
                 else
                     outflow = Float(0)
                 end
@@ -552,17 +555,12 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
                         -1.0 *
                         lake.b[i] *
                         pow((lake.waterlevel[lo] - lake.threshold[i]), lake.e[i])
+                    outflow = max(outflow, -lake.storage[lo])
                 else
                     outflow = Float(0)
                 end
             end
         end
-        storage_input = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow
-        if diff_wl >= 0.0
-            outflow = min(outflow, storage_input)
-        else
-            outflow = max(outflow, -lake.storage[lo])
-        end
         storage = (storage_input - outflow) * timestepsecs
 
         # update storage and outflow for lake with rating curve of type 1.

From 3fa3e166ceccbbef9c1c9c380aed85cf36f53671 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Fri, 19 Jan 2024 09:42:52 +0100
Subject: [PATCH 7/9] Update changelog

---
 docs/src/changelog.md | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index ad8633198..c93b8101f 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 - Added missing BMI function `get_grid_size`, it is used for unstructured grids, for example
   to get the length of arrays returned by BMI functions `get_grid_x` and `get_grid_y`.
+- Added a check for the solution of the quadratic equation as part of the Modified Puls
+  approach for lake outflow. Lower limit should be zero (very small negative values can
+  occur).
+- Limit lake evaporation (added variable `actevap`) and lake outflow to prevent negative
+  lake storage. The variable `actevap` has also been added to the reservoir module.
 
 ## v0.7.3 - 2024-01-12
 

From 162b769a1a898a7f2cbfdf3457db3b12d528da26 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Fri, 19 Jan 2024 10:07:26 +0100
Subject: [PATCH 8/9] Fix limit outflow for lakes with threshold

---
 src/reservoir_lake.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 5e0a2d0d5..6b9d27049 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -543,19 +543,19 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
         else
             if diff_wl >= 0.0
                 if lake.waterlevel[i] > lake.threshold[i]
-                    outflow =
-                        lake.b[i] * pow((lake.waterlevel[i] - lake.threshold[i]), lake.e[i])
-                    outflow = min(outflow, storage_input)
+                    Δh = lake.waterlevel[i] - lake.threshold[i]
+                    outflow = lake.b[i] * pow(Δh, lake.e[i])
+                    maxflow = (Δh * lake.area[i]) / timestepsecs
+                    outflow = min(outflow, maxflow)
                 else
                     outflow = Float(0)
                 end
             else
                 if lake.waterlevel[lo] > lake.threshold[i]
-                    outflow =
-                        -1.0 *
-                        lake.b[i] *
-                        pow((lake.waterlevel[lo] - lake.threshold[i]), lake.e[i])
-                    outflow = max(outflow, -lake.storage[lo])
+                    Δh = lake.waterlevel[lo] - lake.threshold[i]
+                    outflow = -1.0 * lake.b[i] * pow(Δh, lake.e[i])
+                    maxflow = (Δh * lake.area[lo]) / timestepsecs
+                    outflow = max(outflow, -maxflow)
                 else
                     outflow = Float(0)
                 end

From b24bde100222600eb614ce60f06bbb06775c8a08 Mon Sep 17 00:00:00 2001
From: Willem van Verseveld <willem.vanverseveld@gmail.com>
Date: Tue, 30 Jan 2024 20:08:56 +0100
Subject: [PATCH 9/9] Address review comments To clarify parameter
 descriptions.

---
 docs/src/model_docs/params_lateral.md | 4 ++--
 src/reservoir_lake.jl                 | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 236502508..e9bd012bc 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -105,7 +105,7 @@ locs = "wflow_reservoirlocs"
 | `totaloutflow`        | total outflow into reservoir | m``^3`` | - |
 | `percfull`             | fraction full (of max storage) | - | - |
 | `precipitation`             | average precipitation for reservoir area | mm Δt⁻¹ | - |
-| `evaporation`             | average evaporation for reservoir area | mm Δt⁻¹ | - |
+| `evaporation`             | average potential evaporation for reservoir area | mm Δt⁻¹ | - |
 | `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
 
 ### [Lakes](@id lake_params)
@@ -148,7 +148,7 @@ between parentheses.
 | `outflow` | outflow lake | m``^3`` s``^{-1}``  | - |
 | `totaloutflow` | total outflow lake | m``^3``  | - |
 | `precipitation` | average precipitation for lake area | mm Δt⁻¹  | - |
-| `evaporation` | average evaporation for lake area | mm Δt⁻¹  | - |
+| `evaporation` | average potential evaporation for lake area | mm Δt⁻¹  | - |
 | `actevap` | average actual evaporation for lake area | mm Δt⁻¹  | - |
 
 ### [Lateral subsurface flow](@id params_ssf)
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 6b9d27049..110d2edff 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -13,7 +13,7 @@
     percfull::Vector{T} | "-"                           # fraction full (of max storage) [-]
     demandrelease::Vector{T} | "m3 s-1"                 # minimum (environmental) flow released from reservoir [m³ s⁻¹]
     precipitation::Vector{T}                            # average precipitation for reservoir area [mm Δt⁻¹]
-    evaporation::Vector{T}                              # average evaporation for reservoir area [mm Δt⁻¹]
+    evaporation::Vector{T}                              # average potential evaporation for reservoir area [mm Δt⁻¹]
     actevap::Vector{T}                                  # average actual evaporation for reservoir area [mm Δt⁻¹]
 
     function SimpleReservoir{T}(args...) where {T}
@@ -223,7 +223,7 @@ end
     outflow::Vector{T} | "m3 s-1"               # outflow lake [m³ s⁻¹]
     totaloutflow::Vector{T} | "m3"              # total outflow lake [m³]
     precipitation::Vector{T}                    # average precipitation for lake area [mm Δt⁻¹]
-    evaporation::Vector{T}                      # average evaporation for lake area [mm Δt⁻¹]
+    evaporation::Vector{T}                      # average potential evaporation for lake area [mm Δt⁻¹]
     actevap::Vector{T}                          # average actual evapotranspiration for lake area [mm Δt⁻¹]
 
     function Lake{T}(args...) where {T}