Merge pull request #129 from Deltares/kinwave_vars

Kinwave vars

Author: vers-w

docs/src/changelog.md

@@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `restrapefficiency` in the sediment model.
 - New variables added to the `LandSediment` and `RiverSediment` structs in order to save 
   more output from the sediment model.
+- Added variables `volume` and `inwater` to `SurfaceFlow` struct, this is convenient for the
+  coupling with the water quality model Delwaq. 
 
 ### Added
 - Modify model parameters and forcing through the TOML file (see [Modify parameters](@ref)).
@@ -22,6 +24,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 - Corrected a bug in sediment deposition in the river (case when incoming sediment load is
   more than the river transport capacity).
+- Fixed update of `snow` and `glacierstore` fields of HBV and SBM concepts by the
+  `glacier_hbv` module. 
 
 ## v0.2.0 - 2021-03-26
 

src/flow.jl

@@ -6,7 +6,9 @@
     q::Vector{T} | "m3 s-1"                 # Discharge [m³ s⁻¹]
     qin::Vector{T} | "m3 s-1"               # Inflow from upstream cells [m³ s⁻¹]
     q_av::Vector{T} | "m3 s-1"              # Average discharge [m³ s⁻¹]
-    qlat::Vector{T} | "m3 s-1"              # Lateral discharge [m³ s⁻¹]
+    qlat::Vector{T} | "m2 s-1"              # Lateral inflow per unit length [m² s⁻¹]
+    inwater::Vector{T} | "m3 s-1"           # Lateral inflow [m³ s⁻¹]
+    volume::Vector{T} | "m3"                # Kinematic wave volume [m³] (based on water level h)
     h::Vector{T} | "m"                      # Water level [m]
     h_av::Vector{T} | "m"                   # Average water level [m]
     Δt::T | "s"                             # Model time step [s]
@@ -152,6 +154,7 @@ function update(
     sf.q_av .= q_sum ./ its
     sf.h_av .= h_sum ./ its
     sf.to_river .= sf.to_river ./ its
+    sf.volume .= sf.dl .* sf.width .* sf.h
 end
 
 @get_units @with_kw struct LateralSSF{T}

src/hbv.jl

@@ -136,10 +136,10 @@ function update_after_snow(hbv::HBV, config)
             # Estimate the fraction of snow turned into ice (HBV-light).
             # Estimate glacier melt.
 
-            drysnow, snow2glacier, glacierstore, glaciermelt = glacier_hbv(
+            hbv.snow[i], snow2glacier, hbv.glacierstore[i], glaciermelt = glacier_hbv(
                 hbv.glacierfrac[i],
                 hbv.glacierstore[i],
-                hbv.drysnow[i],
+                hbv.snow[i],
                 hbv.temperature[i],
                 hbv.g_tt[i],
                 hbv.g_cfmax[i],

src/hbv_model.jl

@@ -414,6 +414,8 @@ function initialize_hbv_model(config::Config)
         qin = zeros(Float, n),
         q_av = zeros(Float, n),
         qlat = zeros(Float, n),
+        inwater = zeros(Float, n),
+        volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
         Δt = Float(tosecond(Δt)),
@@ -435,12 +437,8 @@ function initialize_hbv_model(config::Config)
     pcr_dir = dims_xy ? permute_indices(pcrdir) : pcrdir
     graph = flowgraph(ldd, inds, pcr_dir)
 
-    riverslope = ncread(
-        nc,
-        param(config, "input.lateral.river.slope");
-        sel = inds_riv,
-        type = Float,
-    )
+    riverslope =
+        ncread(nc, param(config, "input.lateral.river.slope"); sel = inds_riv, type = Float)
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
     riverwidth = riverwidth_2d[inds_riv]
@@ -470,6 +468,8 @@ function initialize_hbv_model(config::Config)
         qin = zeros(Float, nriv),
         q_av = zeros(Float, nriv),
         qlat = zeros(Float, nriv),
+        inwater = zeros(Float, nriv),
+        volume = zeros(Float, nriv),
         h = zeros(Float, nriv),
         h_av = zeros(Float, nriv),
         Δt = Float(tosecond(Δt)),
@@ -534,6 +534,10 @@ function initialize_hbv_model(config::Config)
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
         set_states(instate_path, model, state_ncnames; type = Float)
+        # update kinematic wave volume for river and land domain
+        @unpack lateral = model
+        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
+        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
     end
 
     # make sure the forcing is already loaded
@@ -575,9 +579,9 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:HBV,R,W}
 
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # hbv concept
-    lateral.land.qlat .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.Δt ./ lateral.land.dl
+    lateral.land.inwater .=
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
+    lateral.land.qlat .= lateral.land.inwater ./ lateral.land.dl
 
     # run kinematic wave for overland flow
     update(
@@ -588,7 +592,8 @@ function update(model::Model{N,L,V,R,W}) where {N,L,V<:HBV,R,W}
     )
 
     # determine lateral inflow (from overland flow) for river flow
-    lateral.river.qlat .= (lateral.land.to_river[inds_riv]) ./ lateral.river.dl
+    lateral.river.inwater .= copy(lateral.land.to_river[inds_riv])
+    lateral.river.qlat .= lateral.river.inwater ./ lateral.river.dl
 
     # run kinematic wave for river flow
     # check if reservoirs or lakes are defined, the inflow from overland flow is required

src/reservoir_lake.jl

@@ -6,7 +6,7 @@
     targetminfrac::Vector{T} | "-"                      # target minimum full fraction (of max storage) [-]
     targetfullfrac::Vector{T} | "-"                     # target fraction full (of max storage) [-]
     volume::Vector{T} | "m3"                            # reservoir volume [m³]
-    inflow::Vector{T} | "m3"                        # inflow into reservoir [m³]
+    inflow::Vector{T} | "m3"                            # inflow into reservoir [m³]
     outflow::Vector{T} | "m3 s-1"                       # outflow from reservoir [m³ s⁻¹]
     totaloutflow::Vector{T} | "m3"                      # total outflow from reservoir [m³]
     percfull::Vector{T} | "-"                           # fraction full (of max storage) [-]

src/sbm.jl

@@ -743,7 +743,7 @@ function update_until_recharge(sbm::SBM, config)
                 # Estimate the fraction of snow turned into ice (HBV-light).
                 # Estimate glacier melt.
 
-                snow, snow2glacier, glacierstore, glaciermelt = glacier_hbv(
+                sbm.snow[i], snow2glacier, sbm.glacierstore[i], glaciermelt = glacier_hbv(
                     sbm.glacierfrac[i],
                     sbm.glacierstore[i],
                     sbm.snow[i],

src/sbm_gwf_model.jl

@@ -59,8 +59,7 @@ function initialize_sbm_gwf_model(config::Config)
         ncread(nc, param(config, "input.lateral.river.length"); type = Float, fill = 0)
     riverlength = riverlength_2d[inds]
 
-    altitude =
-        ncread(nc, param(config, "input.altitude"); sel = inds, type = Float)
+    altitude = ncread(nc, param(config, "input.altitude"); sel = inds, type = Float)
     # read x, y coordinates and calculate cell length [m]
     y_nc = "y" in keys(nc.dim) ? ncread(nc, "y") : ncread(nc, "lat")
     x_nc = "x" in keys(nc.dim) ? ncread(nc, "x") : ncread(nc, "lon")
@@ -139,6 +138,8 @@ function initialize_sbm_gwf_model(config::Config)
         qin = zeros(Float, n),
         q_av = zeros(Float, n),
         qlat = zeros(Float, n),
+        inwater = zeros(Float, n),
+        volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
         Δt = Float(tosecond(Δt)),
@@ -161,12 +162,8 @@ function initialize_sbm_gwf_model(config::Config)
     graph = flowgraph(ldd, inds, pcr_dir)
 
     # river flow (kinematic wave)
-    riverslope = ncread(
-        nc,
-        param(config, "input.lateral.river.slope");
-        sel = inds_riv,
-        type = Float,
-    )
+    riverslope =
+        ncread(nc, param(config, "input.lateral.river.slope"); sel = inds_riv, type = Float)
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
     riverwidth = riverwidth_2d[inds_riv]
@@ -193,6 +190,8 @@ function initialize_sbm_gwf_model(config::Config)
         q_av = zeros(Float, nriv),
         qin = zeros(Float, nriv),
         qlat = zeros(Float, nriv),
+        inwater = zeros(Float, nriv),
+        volume = zeros(Float, nriv),
         h = zeros(Float, nriv),
         h_av = zeros(Float, nriv),
         Δt = Float(tosecond(Δt)),
@@ -408,6 +407,10 @@ function initialize_sbm_gwf_model(config::Config)
         instate_path = joinpath(tomldir, config.state.path_input)
         state_ncnames = ncnames(config.state)
         set_states(instate_path, model, state_ncnames, type = Float)
+        # update kinematic wave volume for river and land domain
+        @unpack lateral = model
+        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
+        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
     end
 
     # make sure the forcing is already loaded
@@ -489,9 +492,9 @@ function update_sbm_gwf(model)
 
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # sbm concept
-    lateral.land.qlat .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.Δt ./ lateral.land.dl
+    lateral.land.inwater .=
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
+    lateral.land.qlat .= lateral.land.inwater ./ lateral.land.dl
     # run kinematic wave for overland flow
     update(
         lateral.land,
@@ -516,11 +519,10 @@ function update_sbm_gwf(model)
     end
 
     # determine lateral inflow to river from groundwater, drains and overland flow
-    lateral.river.qlat .=
-        (
-            flux_gw[inds_riv] ./ lateral.river.Δt .+ lateral.land.to_river[inds_riv] .+
-            net_runoff_river
-        ) ./ lateral.river.dl
+    lateral.river.inwater .=
+        flux_gw[inds_riv] ./ lateral.river.Δt .+ lateral.land.to_river[inds_riv] .+
+        net_runoff_river
+    lateral.river.qlat .= lateral.river.inwater ./ lateral.river.dl
 
     # run kinematic wave for river flow 
     # check if reservoirs or lakes are defined, then the inflow from overland flow is

src/sbm_model.jl

@@ -177,6 +177,8 @@ function initialize_sbm_model(config::Config)
         qin = zeros(Float, n),
         q_av = zeros(Float, n),
         qlat = zeros(Float, n),
+        inwater = zeros(Float, n),
+        volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
         Δt = Float(tosecond(Δt)),
@@ -199,12 +201,8 @@ function initialize_sbm_model(config::Config)
     graph = flowgraph(ldd, inds, pcr_dir)
 
 
-    riverslope = ncread(
-        nc,
-        param(config, "input.lateral.river.slope");
-        sel = inds_riv,
-        type = Float,
-    )
+    riverslope =
+        ncread(nc, param(config, "input.lateral.river.slope"); sel = inds_riv, type = Float)
     clamp!(riverslope, 0.00001, Inf)
     riverlength = riverlength_2d[inds_riv]
     riverwidth = riverwidth_2d[inds_riv]
@@ -234,6 +232,8 @@ function initialize_sbm_model(config::Config)
         qin = zeros(Float, nriv),
         q_av = zeros(Float, nriv),
         qlat = zeros(Float, nriv),
+        inwater = zeros(Float, nriv),
+        volume = zeros(Float, nriv),
         h = zeros(Float, nriv),
         h_av = zeros(Float, nriv),
         Δt = Float(tosecond(Δt)),
@@ -310,14 +310,16 @@ function initialize_sbm_model(config::Config)
         state_ncnames = ncnames(config.state)
         set_states(instate_path, model, state_ncnames; type = Float)
         @unpack lateral, vertical = model
-        # update zi for vertical sbm and α for river and overland flow
+        # update zi for vertical sbm and kinematic wave volume for river and land domain
         zi =
             max.(
                 0.0,
                 vertical.soilthickness .-
                 vertical.satwaterdepth ./ (vertical.θₛ .- vertical.θᵣ),
             )
         vertical.zi .= zi
+        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
+        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
     end
 
     # make sure the forcing is already loaded
@@ -405,9 +407,10 @@ function update_after_subsurfaceflow(model::Model{N,L,V,R,W}) where {N,L,V<:SBM,
 
     # determine lateral inflow for overland flow based on vertical runoff [mm] from vertical
     # sbm concept
-    lateral.land.qlat .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.Δt ./ lateral.land.dl
+    lateral.land.inwater .=
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
+    lateral.land.qlat .= lateral.land.inwater ./ lateral.land.dl
+
     # run kinematic wave for overland flow
     update(
         lateral.land,
@@ -423,11 +426,11 @@ function update_after_subsurfaceflow(model::Model{N,L,V,R,W}) where {N,L,V<:SBM,
             vertical.net_runoff_river[inds_riv] .* network.land.xl[inds_riv] .*
             network.land.yl[inds_riv] .* 0.001
         ) ./ vertical.Δt
-    lateral.river.qlat .=
-        (
-            lateral.subsurface.to_river[inds_riv] ./ 1.0e9 ./ lateral.river.Δt .+
-            lateral.land.to_river[inds_riv] .+ net_runoff_river
-        ) ./ lateral.river.dl
+    lateral.river.inwater .= (
+        lateral.subsurface.to_river[inds_riv] ./ 1.0e9 ./ lateral.river.Δt .+
+        lateral.land.to_river[inds_riv] .+ net_runoff_river
+    )
+    lateral.river.qlat .= lateral.river.inwater ./ lateral.river.dl
 
     # run kinematic wave for river flow
     # check if reservoirs or lakes are defined, the inflow from lateral subsurface and

test/bmi.jl

@@ -18,9 +18,9 @@ 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) == 175
-            @test BMI.get_output_item_count(model) == 175
-            @test BMI.get_input_var_names(model)[[1, 5, 120, 175]] == [
+            @test BMI.get_input_item_count(model) == 179
+            @test BMI.get_output_item_count(model) == 179
+            @test BMI.get_input_var_names(model)[[1, 5, 120, 179]] == [
                 "vertical.Δt",
                 "vertical.n_unsatlayers",
                 "lateral.land.q_av",
@@ -35,7 +35,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test_throws ErrorException BMI.get_var_grid(model, "lateral.river.lake.volume")
             # Vector{Float64} printing on Julia 1.6+
             @test BMI.get_var_type(model, "lateral.river.reservoir.inflow") in
-                ["Array{$Float,1}", "Vector{$Float}"]
+                  ["Array{$Float,1}", "Vector{$Float}"]
             @test BMI.get_var_type(model, "vertical.n") == "Int64"
             @test BMI.get_var_units(model, "vertical.θₛ") == "mm mm-1"
             @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == sizeof(Float)

test/io.jl

@@ -184,11 +184,14 @@ end
     @test Wflow.param(model, "vertical.ustorelayerdepth")[1][1] ≈ 16.73013687133789
     @test Wflow.param(model, "vertical.snowwater")[41120] ≈ 0.42188167572021484
     @test Wflow.param(model, "vertical.canopystorage")[1] ≈ 0.0
+    @test Wflow.param(model, "vertical.zi")[1] ≈ 380.67793082060416
     @test Wflow.param(model, "lateral.subsurface.ssf")[39308] ≈ 1.1614302208e11
     @test Wflow.param(model, "lateral.river.q")[5501] ≈ 111.46229553222656
     @test Wflow.param(model, "lateral.river.h")[5501] ≈ 8.555977821350098
+    @test Wflow.param(model, "lateral.river.volume")[5501] ≈ 249163.33754433296
     @test Wflow.param(model, "lateral.land.q")[39626] ≈ 1.0575231313705444
     @test Wflow.param(model, "lateral.land.h")[39626] ≈ 0.03456519544124603
+    @test Wflow.param(model, "lateral.land.volume")[39626] ≈ 19379.23274404319
 end
 
 @testset "reducer" begin

test/run_hbv.jl

@@ -45,10 +45,13 @@ model = Wflow.update(model)
     @test hbv.snow[1] == 0.0
 end
 
-@testset "overland flow" begin
+@testset "overland domain" begin
     q = model.lateral.land.q_av
+    land = model.lateral.land
     @test sum(q) ≈ 3278.9423039643552f0
     @test q[500] ≈ 0.2348646912841589f0
+    @test land.volume[500] ≈ 2057.7314802432425f0
+    @test land.inwater[500] ≈ 0.027351853491789667f0
     @test q[1566] ≈ 0.030463805623037462f0
     @test q[network.land.order[end]] ≈ 0.296794455575309f0
 end

test/run_sbm_gwf.jl

@@ -43,10 +43,13 @@ end
     @test sum(q) ≈ 0.0
 end
 
-@testset "river flow" begin
+@testset "river domain" begin
     q = model.lateral.river.q_av
+    river = model.lateral.river
     @test sum(q) ≈ 0.03515257228971982f0
     @test q[6] ≈ 0.007952670041402076f0
+    @test river.volume[6] ≈ 3.8923529548071025f0
+    @test river.inwater[6] ≈ 0.0003950369148075657f0
     @test q[13] ≈ 0.0005980549969548235f0
     @test q[network.river.order[end]] ≈ 0.008482648782941154f0
 end