SBM evaporation and use of fractions (#381)

* Use of fractions for potential evaporation
Always use the fraction (`riverfrac`, `waterfrac`, `canopygapfraction` and bare soil) to divide potential evaporation over different processes (interception, transpiration, soil evaporation and open water evaporation).

* Replace `et_reftopot` by crop coefficient `kc`

* Update docs and fix `ewet` and `e_r` computation
For computation of `ewet` the crop coefficient `kc` should be used, and for `e_r` the canopy fraction should be used for the precipitation amount.

* Update changelog

* Soil evaporation and `glacierfrac`
Besides `riverfrac` and `waterfrac` also include `glacierfrac` to correct `canopygapfraction` for soil evaporation.

docs/src/changelog.md

@@ -22,6 +22,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   `to_river` variable from overland flow and lateral subsurface flow was not added to the
   inflow of these locations.
 - Close netCDF `NCDataset` with state variables in extended BMI function `save_state`.
+- For the computation of Gash interception model parameter `e_r` multiply the precipitation
+  input with the canopy fraction (this was only done for the potential evapotranspiration
+  input).
+
 
 ### Changed
 - Stop exposing scalar variables through BMI. The `BMI.get_value_ptr` function was not
@@ -46,6 +50,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   When in the Python version of wflow\_sbm the kinematic wave for surface flow was split
   into a river and land component, the lower limit was removed for river runoff
   (`net_runoff_river`), but not for land runoff.
+- Always use fractions for the computation of potential evapotranspiration (interception and
+  transpiration) and potential evaporation (open water and soil). Replaced variable
+  `et_reftopot` by crop coefficient `kc` (use of `et_reftopot` has been deprecated).
 
 ### Added
 - Total water storage as an export variable for `SBM` concept. This is the total water stored

docs/src/model_docs/params_vertical.md

@@ -64,7 +64,7 @@ profile `kv` is used and `z_layered` is required as input.
 | **`rootdistpar`** | controls how roots are linked to water table | - | -500.0  |
 | **`cap_hmax`** | water depth beyond which capillary flux ceases | mm | 2000.0  |
 | **`cap_n`** | coefficient controlling capillary rise | - | 2.0  |
-| **`et_reftopot`** | multiplication factor to correct reference evaporation | - | 1.0  |
+| **`kc`** | crop coefficient Kc | - | 1.0  |
 | **`sl`** (`specific_leaf`) | specific leaf storage  | mm | - |
 | **`swood`** (`storage_wood`) | storage woody part of vegetation | mm | - |
 | **`kext`** | extinction coefficient (to calculate canopy gap fraction) | - | - |
@@ -89,8 +89,8 @@ profile `kv` is used and `z_layered` is required as input.
 | `canopystorage`|  canopy storage | mm | - |
 |**`precipitation`** | precipitation | mm Δt``^{-1}``|  - |
 | **`temperature`** | temperature | ᵒC | - |
-| **`potential_evaporation`** | potential evaporation | mm Δt``^{-1}`` | - |
-| `pottrans_soil` | interception subtracted from potential evaporation) | mm Δt``^{-1}`` | - |
+| **`potential_evaporation`** | potential reference evapotranspiration | mm Δt``^{-1}`` | - |
+| `pottrans` | interception subtracted from potential evapotranspiration | mm Δt``^{-1}`` | - |
 | `transpiration` | transpiration | mm Δt``^{-1}`` | - |
 | `ae_ustore` | actual evaporation from unsaturated store | mm Δt``^{-1}`` | - |
 | `interception` | interception loss by evaporation | mm Δt``^{-1}`` | - |

docs/src/model_docs/vertical/sbm.md

@@ -134,22 +134,32 @@ The extinction coefficient `kext` can be related to land cover.
 
 ## [Evaporation](@id evap)
 
-The wflow\_sbm model assumes the input to be potential evaporation. A multiplication factor
-(`et_reftopot`, set to 1 by default) is present to correct the input evaporation if
-required.
-
-The potential evaporation left over after interception and open water evaporation (rivers
-and water bodies) is split in potential soil evaporation and potential transpiration based
-on the canopy gap fraction (assumed to be identical to the amount of bare soil).
+The wflow\_sbm model assumes the input to be potential reference evapotranspiration. A crop
+coefficient (`kc`, set to 1 by default) is used to convert the potential evapotranspiration
+rate of a reference crop fully covering the soil to the potential evapotranspiration rate of
+vegetation (natural and agricultural) fully covering the soil. The crop coefficient `kc` of
+wflow\_sbm is used for a surface completely covered by vegetation, and does not include the
+effect of growing stages of vegetation and soil cover. These effects are handled separately
+through the use of the canopy gap fraction. 
+
+It is assumed that the potential evaporation rate of intercepted water by vegetation is
+equal to the potential evapotranspiration rate of vegetation (fully covering the soil)
+multiplied by the canopy fraction. The potential evapotranspiration rate left over after
+interception is available for transpiration. For potential open water evaporation (river and
+water bodies) the potential reference evapotranspiration rate is used (multipled by the
+river fraction `riverfrac`, and open water fraction `waterfrac`). Also for potential soil
+evaporation the potential reference evapotranspiration rate is used, multiplied by the
+canopy gap fraction corrected by the sum of total water fraction (`riverfrac` and
+`waterfrac`) and the fraction covered by a glacier (`glacierfrac`).
 
 ### Bare soil evaporation
 
 If there is only one soil layer present in the wflow\_sbm model, the bare soil evaporation
 is scaled according to the wetness of the soil layer. The fraction of bare soil is assumed
-to be equal to the fraction not covered by the canopy (`canopygapfraction`). When the soil
-is fully saturated, evaporation is set to equal the potential evaporation. When the soil is
-not fully saturated, actual evaporation decreases linearly with decreasing soil moisture
-values, as indicated by the figure below.
+to be equal to the fraction not covered by the canopy (`canopygapfraction`) corrected by the
+total water fraction. When the soil is fully saturated, evaporation is set to equal the
+potential reference evaporation. When the soil is not fully saturated, actual evaporation
+decreases linearly with decreasing soil moisture values, as indicated by the figure below.
 
 ![soil_evap](../../images/soil_evap.png)
 
@@ -182,10 +192,10 @@ index of the vector that contains all active cells within the spatial model doma
 
 ```julia
     # transpiration from saturated store
-    wetroots = scurve(sbm.zi[i], rootingdepth, 1.0, sbm.rootdistpar[i])
-    actevapsat = min(pottrans * wetroots, satwaterdepth)
+    wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
+    actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
     satwaterdepth = satwaterdepth - actevapsat
-    restpottrans = pottrans - actevapsat
+    restpottrans = sbm.pottrans[i] - actevapsat
 ```
 
 ![soil_wetroots](../../images/soil_wetroots.png)

server/test/client.jl

@@ -15,8 +15,7 @@ function request(message)
 end
 
 @testset "initialization and time functions" begin
-    msg =
-        (fn = "initialize", config_file = joinpath(@__DIR__, "sbm_config.toml"))
+    msg = (fn = "initialize", config_file = joinpath(@__DIR__, "sbm_config.toml"))
     @test request(msg) == Dict("status" => "OK")
     @test request((fn = "get_end_time",)) == Dict("end_time" => 2678400)
     @test request((fn = "get_start_time",)) == Dict("start_time" => 0)
@@ -67,7 +66,7 @@ vwc_1_size = 0
     vwc_1_size = Int(vwc_1_nbytes / vwc_1_itemsize)
     @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
-    @test mean(request(msg)["value"]) ≈ 278.1510965581235
+    @test mean(request(msg)["value"]) ≈ 277.8362161218515
     msg = (fn = "get_value_ptr", name = "vertical.θₛ")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
@@ -77,7 +76,7 @@ vwc_1_size = 0
         inds = [1, 5, 10],
     )
     @test request(msg)["value_at_indices"] ≈
-          [2.1901434445889123, 2.6778265820894545, 3.470059871798756]
+          [2.196562540141252f0, 2.6853997213000866f0, 3.4814142040324985f0]
     msg = (fn = "set_value", name = "vertical.zi", src = fill(300.0, zi_size))
     @test request(msg) == Dict("status" => "OK")
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
@@ -97,15 +96,15 @@ vwc_1_size = 0
     )
     @test request(msg)["value_at_indices"] == [250.0, 350.0, 300.0]
     msg = (fn = "get_value", name = "vertical.vwc[1]", dest = fill(0.0, vwc_1_size))
-    @test mean(request(msg)["value"]) ≈ 0.1845159140308566f0
+    @test mean(request(msg)["value"]) ≈ 0.18600013563085036f0
     msg = (
         fn = "get_value_at_indices",
         name = "vertical.vwc[1]",
         dest = [0.0, 0.0, 0.0],
         inds = [1, 2, 3],
     )
     @test request(msg)["value_at_indices"] ≈
-          [0.12089607119560242f0, 0.11967185322648062f0, 0.14503555864288548f0]
+          [0.12089607119560242f0, 0.11968416924304527f0, 0.14602328618707333f0]
     msg = (fn = "set_value", name = "vertical.vwc[1]", src = fill(0.3, vwc_1_size))
     @test request(msg) == Dict("status" => "OK")
     msg = (fn = "get_value", name = "vertical.vwc[1]", dest = fill(0.0, vwc_1_size))

src/sbm.jl

@@ -49,8 +49,8 @@
     cap_hmax::Vector{T} | "mm"
     # Coefficient [-] controlling capillary rise
     cap_n::Vector{T} | "-"
-    # Multiplication factor [-] to correct
-    et_reftopot::Vector{T} | "-"
+    # Crop coefficient Kc [-]
+    kc::Vector{T} | "-"
     # Brooks-Corey power coefficient [-] for each soil layer
     c::Vector{SVector{N,T}} | "-"
     # Stemflow [mm Δt⁻¹]
@@ -84,10 +84,10 @@
     precipitation::Vector{T}
     # Temperature [ᵒC]
     temperature::Vector{T} | "°C"
-    # Potential evapotranspiration [mm Δt⁻¹]
+    # Potential reference evapotranspiration [mm Δt⁻¹]
     potential_evaporation::Vector{T}
-    # Potential transpiration, open water, river and soil evaporation (after subtracting interception from potential_evaporation)
-    pottrans_soil::Vector{T}
+    # Potential transpiration (after subtracting interception from potential_evaporation)
+    pottrans::Vector{T}
     # Transpiration [mm Δt⁻¹]
     transpiration::Vector{T}
     # Actual evaporation from unsaturated store [mm Δt⁻¹]
@@ -469,8 +469,21 @@ function initialize_sbm(nc, config, riverfrac, inds)
     cap_hmax =
         ncread(nc, config, "vertical.cap_hmax"; sel = inds, defaults = 2000.0, type = Float)
     cap_n = ncread(nc, config, "vertical.cap_n"; sel = inds, defaults = 2.0, type = Float)
-    et_reftopot =
-        ncread(nc, config, "vertical.et_reftopot"; sel = inds, defaults = 1.0, type = Float)
+    kc = ncread(
+        nc,
+        config,
+        "vertical.kc";
+        alias = "vertical.et_reftopot",
+        sel = inds,
+        defaults = 1.0,
+        type = Float,
+    )
+    if haskey(config.input.vertical, "et_reftopot")
+        @warn string(
+            "The `et_reftopot` key in `[input.vertical]` is now called ",
+            "`kc`. Please update your TOML file.",
+        )
+    end
 
     cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
 
@@ -591,7 +604,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         rootdistpar = rootdistpar,
         cap_hmax = cap_hmax,
         cap_n = cap_n,
-        et_reftopot = et_reftopot,
+        kc = kc,
         c = svectorscopy(c, Val{maxlayers}()),
         stemflow = fill(mv, n),
         throughfall = fill(mv, n),
@@ -609,7 +622,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
         precipitation = fill(mv, n),
         temperature = fill(mv, n),
         potential_evaporation = fill(mv, n),
-        pottrans_soil = fill(mv, n),
+        pottrans = fill(mv, n),
         transpiration = fill(mv, n),
         ae_ustore = fill(mv, n),
         interception = fill(mv, n),
@@ -688,42 +701,39 @@ function update_until_snow(sbm::SBM, config)
         if do_lai
             cmax = sbm.sl[i] * sbm.leaf_area_index[i] + sbm.swood[i]
             canopygapfraction = exp(-sbm.kext[i] * sbm.leaf_area_index[i])
-            ewet =
-                (1.0 - exp(-sbm.kext[i] * sbm.leaf_area_index[i])) *
-                sbm.potential_evaporation[i]
+            canopyfraction = 1.0 - canopygapfraction
+            ewet = canopyfraction * sbm.potential_evaporation[i] * sbm.kc[i]
             e_r =
                 sbm.precipitation[i] > 0.0 ?
-                min(0.25, ewet / max(0.0001, sbm.precipitation[i])) : 0.0
+                min(0.25, ewet / max(0.0001, canopyfraction * sbm.precipitation[i])) : 0.0
         else
             cmax = sbm.cmax[i]
             canopygapfraction = sbm.canopygapfraction[i]
             e_r = sbm.e_r[i]
         end
 
-        potential_evaporation = sbm.potential_evaporation[i] * sbm.et_reftopot[i]
-        # should we include tempcor in SBM?
-        # potential_evaporation = PotenEvap #??
-
+        canopy_potevap =
+            sbm.kc[i] * sbm.potential_evaporation[i] * (1.0 - canopygapfraction)
         if Second(sbm.Δt) >= Hour(23)
             throughfall, interception, stemflow, canopystorage = rainfall_interception_gash(
                 cmax,
                 e_r,
                 canopygapfraction,
                 sbm.precipitation[i],
                 sbm.canopystorage[i],
-                potential_evaporation,
+                canopy_potevap,
             )
-            pottrans_soil = max(0.0, potential_evaporation - interception) # now in mm
+            pottrans = max(0.0, canopy_potevap - interception) # now in mm
         else
             netinterception, throughfall, stemflow, leftover, interception, canopystorage =
                 rainfall_interception_modrut(
                     sbm.precipitation[i],
-                    potential_evaporation,
+                    canopy_potevap,
                     sbm.canopystorage[i],
                     canopygapfraction,
                     cmax,
                 )
-            pottrans_soil = max(0.0, leftover)  # now in mm
+            pottrans = max(0.0, leftover)  # now in mm
         end
 
         if modelsnow
@@ -748,7 +758,7 @@ function update_until_snow(sbm::SBM, config)
         sbm.interception[i] = interception
         sbm.stemflow[i] = stemflow
         sbm.throughfall[i] = throughfall
-        sbm.pottrans_soil[i] = pottrans_soil
+        sbm.pottrans[i] = pottrans
         if modelsnow
             sbm.snow[i] = snow
             sbm.snowwater[i] = snowwater
@@ -806,18 +816,20 @@ function update_until_recharge(sbm::SBM, config)
 
         ae_openw_r = min(
             sbm.waterlevel_river[i] * sbm.riverfrac[i],
-            sbm.riverfrac[i] * sbm.pottrans_soil[i],
+            sbm.riverfrac[i] * sbm.potential_evaporation[i],
         )
         ae_openw_l = min(
             sbm.waterlevel_land[i] * sbm.waterfrac[i],
-            sbm.waterfrac[i] * sbm.pottrans_soil[i],
+            sbm.waterfrac[i] * sbm.potential_evaporation[i],
         )
 
-        restevap = sbm.pottrans_soil[i] - ae_openw_r - ae_openw_l
-
-        # evap available for soil evaporation and transpiration
-        potsoilevap = restevap * sbm.canopygapfraction[i]
-        pottrans = restevap * (1.0 - sbm.canopygapfraction[i])
+        # evap available for soil evaporation
+        soilevap_fraction = max(
+            sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.waterfrac[i] -
+            sbm.glacierfrac[i],
+            0.0,
+        )
+        potsoilevap = soilevap_fraction * sbm.potential_evaporation[i]
 
         # Calculate the initial capacity of the unsaturated store
         ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth
@@ -926,9 +938,9 @@ function update_until_recharge(sbm::SBM, config)
 
         # transpiration from saturated store
         wetroots = scurve(sbm.zi[i], rootingdepth, Float(1.0), sbm.rootdistpar[i])
-        actevapsat = min(pottrans * wetroots, satwaterdepth)
+        actevapsat = min(sbm.pottrans[i] * wetroots, satwaterdepth)
         satwaterdepth = satwaterdepth - actevapsat
-        restpottrans = pottrans - actevapsat
+        restpottrans = sbm.pottrans[i] - actevapsat
 
         # actual transpiration from ustore
         actevapustore = 0.0

test/bmi.jl

@@ -57,7 +57,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test_throws ErrorException BMI.get_value_ptr(model, "vertical.")
             dest = zeros(Float, size(model.vertical.zi))
             BMI.get_value(model, "vertical.zi", dest)
-            @test mean(dest) ≈ 276.3767651555451
+            @test mean(dest) ≈ 276.16325589542333
             @test BMI.get_value_at_indices(
                 model,
                 "vertical.vwc[1]",
@@ -75,7 +75,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
                 "lateral.river.q",
                 zeros(Float, 3),
                 [1, 100, 5617],
-            ) ≈ [0.6211503865184697, 5.219305686635002, 0.026163746306482282]
+            ) ≈ [0.623325399343309, 5.227139951657074, 0.02794287432778194]
             BMI.set_value(model, "vertical.zi", fill(300.0, length(model.vertical.zi)))
             @test mean(
                 BMI.get_value(model, "vertical.zi", zeros(Float, size(model.vertical.zi))),
@@ -150,10 +150,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
 
         @testset "recharge part of SBM" begin
             sbm = model.vertical
-            @test sbm.interception[1] ≈ 0.6299999952316284f0
+            @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.1912317735997524f0
-            @test sbm.recharge[5] ≈ -0.0727941579914808f0
+            @test sbm.snow[1] ≈ 3.484789961176288f0
+            @test sbm.recharge[5] ≈ -0.0f0
             @test sbm.zi[5] ≈ 300.0f0
         end
 
@@ -174,10 +174,10 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         @testset "SBM after subsurface flow" begin
             sbm = model.vertical
             sub = model.lateral.subsurface
-            @test sbm.interception[1] ≈ 0.6299999952316284
+            @test sbm.interception[1] ≈ 0.32734913737568716f0
             @test sbm.ustorelayerdepth[1][1] ≈ 0.0f0
-            @test sbm.snow[1] ≈ 3.1912317735997524f0
-            @test sbm.recharge[5] ≈ -0.0727941579914808f0
+            @test sbm.snow[1] ≈ 3.484789961176288f0
+            @test sbm.recharge[5] ≈ 0.0f0
             @test sbm.zi[5] ≈ 250.0f0
             @test sub.zi[5] ≈ 0.25f0
             @test sub.exfiltwater[1] ≈ 1.0f-5