Check lake balance

docs/src/changelog.md

@@ -10,13 +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`.
-
-### Changed
-
-### Added
-- Total water storage as an export variable for `SBM` concept. This is the total water stored
-  per grid cell in millimeters. Excluded from this variable are the floodplain, lakes and
-  reservoirs.
+- 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.
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was
@@ -25,6 +23,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   required (e.g. for model coupling) while code changes for these variables (including
   struct fields) are required.
 
+### Added
+- Total water storage as an export variable for `SBM` concept. This is the total water stored
+  per grid cell in millimeters. Excluded from this variable are the floodplain, lakes and
+  reservoirs.
+
 ## v0.7.3 - 2024-01-12
 
 ### Fixed

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

src/flow.jl

@@ -266,10 +266,12 @@ 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
         sf.lake.totaloutflow .= 0.0
+        sf.lake.actevap .= 0.0
     end
 
     Δt, its = stable_timestep(sf)
@@ -898,10 +900,12 @@ 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
         sw.lake.totaloutflow .= 0.0
+        sw.lake.actevap .= 0.0
     end
     if !isnothing(sw.floodplain)
         sw.floodplain.q_av .= 0.0
@@ -1138,10 +1142,12 @@ 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
         swr.lake.totaloutflow .= 0.0
+        swr.lake.actevap .= 0.0
     end
     swr.q_av .= 0.0
     swr.h_av .= 0.0

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,14 @@ 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
+    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
@@ -199,6 +200,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
@@ -222,6 +224,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 +418,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,31 +499,33 @@ 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 *
-                (-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
@@ -529,35 +535,33 @@ 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])
+                    Δ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])
+                    Δ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
             end
         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] -
-            outflow * timestepsecs
+        storage = (storage_input - outflow) * timestepsecs
 
         # update storage and outflow for lake with rating curve of type 1.
         if lake.outflowfunc[i] == 1
@@ -584,6 +588,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
@@ -596,6 +601,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

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) == 180
-            @test BMI.get_output_item_count(model) == 180
+            @test BMI.get_input_item_count(model) == 181
+            @test BMI.get_output_item_count(model) == 181
             @test BMI.get_input_var_names(model)[[1, 5, 151, 175]] == [
                 "vertical.nlayers",
                 "vertical.θᵣ",

test/reservoir_lake.jl

@@ -12,6 +12,7 @@
         targetminfrac = [0.2425554726620697],
         precipitation = [4.2],
         evaporation = [1.5],
+        actevap = [0.0],
         outflow = [NaN],
         percfull = [NaN],
         demandrelease = [NaN],
@@ -25,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
@@ -46,6 +48,7 @@ end
         waterlevel = [18.5],
         precipitation = [20.0],
         evaporation = [3.2],
+        actevap = [0.0],
         outflow = [NaN],
     )
 
@@ -56,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")
@@ -90,6 +94,7 @@ sh = [
         waterlevel = [395.03027, 394.87833],
         precipitation = [10.0, 10.0],
         evaporation = [2.0, 2.0],
+        actevap = [0.0, 0.0],
         outflow = [NaN, NaN],
         storage = Wflow.initialize_storage(
             [2, 2],
@@ -105,14 +110,16 @@ 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)
     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]
 end
 
 @testset "overflowing lake with sh and hq" begin
@@ -139,6 +146,7 @@ end
         waterlevel = [397.75],
         precipitation = [10.0],
         evaporation = [2.0],
+        actevap = [0.0],
         outflow = [NaN],
         storage = [410_760_000],
     )