Total Water Storage

docs/src/model_docs/params_vertical.md

@@ -116,7 +116,7 @@ specific_leaf = "Sl"
 | **`leaf_area_index`** | leaf area index | m``^2`` m``{-2}`` | - |
 | `waterlevel_land` | water level land | mm | - |
 | `waterlevel_river` | water level river | mm | - |
-| `total_storage` | total water storage (excluding floodplains) | mm | - |
+| `total_storage` | total water storage (excluding floodplains, lakes and reservoirs) | mm | - |
 
 
 ## [HBV](@id params_hbv)

src/sbm.jl

@@ -194,7 +194,7 @@
     waterlevel_land::Vector{T} | "mm"
     # Water level river [mm]
     waterlevel_river::Vector{T} | "mm"
-    # Total water storage (excluding floodplain volume) [mm]
+    # Total water storage (excluding floodplain volume, lakes and reservoirs) [mm]
     total_storage::Vector{T} | "mm"
 
     function SBM{T,N,M}(args...) where {T,N,M}
@@ -1053,34 +1053,62 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
     end
 end
 
-function update_total_water_storage(sbm::SBM, rnet, rrouting, lrouting)
+@doc """
+Update the total water storage per cell at the end of a timestep.
+
+Takes the following parameters:
+- sbm:
+    The vertical concept (SBM struct)
+- river_network:
+    The indices of the river cells in relation to the active cells, i.e. model.network.index_river
+- cell_xsize:
+    Size in X direction of the cells acquired from model.network.land.xl
+- cell_ysize:
+    Size in Y direction of the cells acquired from model.network.land.yl
+- river_routing:
+    The river routing struct, i.e. model.lateral.river
+- land_routing:
+    The land routing struct, i.e. model.lateral.land
+"""
+function update_total_water_storage(
+    sbm::SBM, 
+    river_network, 
+    cell_xsize, 
+    cell_ysize, 
+    river_routing, 
+    land_routing
+)
+    # Get length active river cells
+    nriv = length(river_network)
 
     # Set the total storage to zero
     fill!(sbm.total_storage, 0)
 
     # Burn the river routing values
-    sbm.total_storage[rnet] = (
-        rrouting.h_av .* sbm.riverfrac[rnet] * 0.001 # convert to mm
+    sbm.total_storage[river_network] = (
+        ( river_routing.h_av[1:nriv] .* river_routing.width[1:nriv] .* 
+        river_routing.dl[1:nriv] ) ./
+        ( cell_xsize[river_network] .* cell_ysize[river_network] ) * 
+        1000 # Convert to mm
     )
 
     # Chunk the data for parallel computing
     threaded_foreach(1:sbm.n, basesize=1000) do i
 
         # Cumulate per vertical type
         # Maybe re-categorize in the future
-        # Unique are those that are not always present
-        unique = sbm.glacierstore[i] * sbm.glacierfrac[i]
         surface = (
+            sbm.glacierstore[i] * sbm.glacierfrac[i] +
             sbm.snow[i] + sbm.snowwater[i] + sbm.canopystorage[i]
         )
-        sub_surfuce = sbm.ustoredepth[i] + sbm.satwaterdepth[i]
+        sub_surface = sbm.ustoredepth[i] + sbm.satwaterdepth[i]
         lateral = (
-            lrouting.h_av[i] * (1 - sbm.riverfrac[i]) * 0.001 # convert to mm
+            land_routing.h_av[i] * (1 - sbm.riverfrac[i]) * 1000 # convert to mm
         )
 
         # Add everything to the total water storage
         sbm.total_storage[i] += (
-            unique + surface + sub_surfuce + lateral
+            surface + sub_surface + lateral
         )
     end
 end

src/sbm_model.jl

@@ -439,6 +439,11 @@ function update_after_subsurfaceflow(
     return model
 end
 
+@doc """
+Update of the total water storage at the end of each timestep per model cell.
+
+This is done here at model level.
+"""
 function update_total_water_storage(
     model::Model{N,L,V,R,W,T},
 ) where {N,L,V,R,W,T<:SbmModel}
@@ -449,6 +454,8 @@ function update_total_water_storage(
     update_total_water_storage(
         vertical,
         network.index_river,
+        network.land.xl,
+        network.land.yl,
         lateral.river,
         lateral.land,
     )

test/run_sbm.jl

@@ -81,7 +81,7 @@ end
     @test sbm.soilevap[50063] == 0.0
     @test sbm.snow[5] ≈ 3.592840840467347f0
     @test sbm.total_storage[50063] ≈ 559.70849973344f0
-    @test sbm.total_storage[429] ≈ 583.9800667536956f0 # river cell
+    @test sbm.total_storage[429] ≈ 597.4578475404879f0 # river cell
 end
 
 # run the second timestep
@@ -95,7 +95,7 @@ model = Wflow.run_timestep(model)
     @test sbm.soilevap[50063] ≈ 0.006358004660566856f0
     @test sbm.snow[5] ≈ 3.667748983774868f0
     @test sbm.total_storage[50063] ≈ 559.7935411649405f0
-    @test sbm.total_storage[429] ≈ 589.7026699582976f0 # river cell
+    @test sbm.total_storage[429] ≈ 617.0062092646873f0 # river cell
 end
 
 @testset "subsurface flow" begin