Merge branch 'master' into cyclic_times

Project.toml

@@ -1,7 +1,7 @@
 name = "Wflow"
 uuid = "d48b7d99-76e7-47ae-b1d5-ff0c1cf9a818"
 authors = ["Deltares and contributors"]
-version = "0.7.1"
+version = "0.7.2"
 
 [deps]
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"

build/create_binaries/create_app.jl

@@ -18,7 +18,7 @@ create_app(
     executables=["wflow_cli" => "julia_main"],
     precompile_execution_file="precompile.jl",
     filter_stdlibs=false,
-    force=true
+    force=true,
 )
 
 include("add_metadata.jl")

build/wflow_cli/Setup/Installer/Application.Setup.wixproj

@@ -6,7 +6,7 @@
     <ProductVersion>3.10</ProductVersion>
     <ProjectGuid>{39da8083-e405-42e6-850c-d25685e91f81}</ProjectGuid>
     <SchemaVersion>2.0</SchemaVersion>
-    <ReleaseVersion>0.7.1</ReleaseVersion>
+    <ReleaseVersion>0.7.2</ReleaseVersion>
     <OutputType>Package</OutputType>
     <DefineSolutionProperties>false</DefineSolutionProperties>
   </PropertyGroup>

docs/src/changelog.md

@@ -5,12 +5,27 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+
 ## [unreleased]
 
 ### Changed
 - For cyclic parameters different cyclic time inputs are supported (only one common cyclic
   time (for example daily or monthly) was allowed).
 
+## v0.7.2 - 2023-09-27
+
+### Fixed
+- Water balance of modified Rutter interception model. The sum of the stemflow partitioning
+  coefficient `pt` and free throughfall coefficient `p` could get larger than 1, resulting
+  in an overestimation of stemflow and throughfall amounts and negative net interception
+  amounts. And the first drainage amount `dd`, controlled by a change over time in canopy
+  storage capacity `cmax`, should not be subtracted from precipitation to compute net
+  interception.
+- The `netinterception` (net interception) computed by the modified Rutter interception
+  model was stored as `interception` in `SBM`, while this should be the `interception`
+  (interception loss by evaporation) output of the modified Rutter interception model. The
+  `interception` of `SBM` is used to compute the total actual evapotranspiration `actevap`.
+
 ## v0.7.1 - 2023-06-30
 
 ### Fixed

docs/src/model_docs/params_vertical.md

@@ -76,7 +76,7 @@ specific_leaf = "Sl"
 | `pottrans_soil` | interception subtracted from potential evaporation) | mm Δt``^{-1}`` | - |
 | `transpiration` | transpiration | mm Δt``^{-1}`` | - |
 | `ae_ustore` | actual evaporation from unsaturated store | mm Δt``^{-1}`` | - |
-| `interception` | interception | mm Δt``^{-1}`` | - |
+| `interception` | interception loss by evaporation | mm Δt``^{-1}`` | - |
 | `soilevap` | total soil evaporation from unsaturated and saturated store | mm Δt``^{-1}`` | - |
 | `soilevapsat` | soil evaporation from saturated store | mm Δt``^{-1}`` | - |
 | `actcapflux` | actual capillary rise | mm Δt``^{-1}`` | - |

src/sbm.jl

@@ -84,7 +84,7 @@
     transpiration::Vector{T}
     # Actual evaporation from unsaturated store [mm Δt⁻¹]
     ae_ustore::Vector{T}
-    # Interception [mm Δt⁻¹]
+    # Interception loss by evaporation [mm Δt⁻¹]
     interception::Vector{T}
     # Soil evaporation from unsaturated and saturated store [mm Δt⁻¹]
     soilevap::Vector{T}
@@ -650,7 +650,6 @@ function update_until_snow(sbm::SBM, config)
                     cmax,
                 )
             pottrans_soil = max(0.0, leftover)  # now in mm
-            interception = netinterception
         end
 
         if modelsnow

src/vertical_process.jl

@@ -32,9 +32,9 @@ function rainfall_interception_gash(
 )
     # TODO: add other rainfall interception method (lui)
     # TODO: include subdaily Gash model
-    # Hack for stemflow (pt)
-    pt = 0.1 * canopygapfraction
-    pfrac = max((1.0 - canopygapfraction - pt), 0.0)
+    # TODO: improve computation of stemflow partitioning coefficient pt (0.1 * canopygapfraction)
+    pt = min(0.1 * canopygapfraction, 1.0 - canopygapfraction)
+    pfrac = 1.0 - canopygapfraction - pt
     p_sat = (-cmax / e_r) * log(1.0 - min((e_r / pfrac), 1.0))
     p_sat = isinf(p_sat) ? 0.0 : p_sat
 
@@ -58,11 +58,11 @@ function rainfall_interception_gash(
     stemflow = cmaxzero ? 0.0 : stemflow
 
     # Now corect for maximum potential evap
-    overestimate = interception > maxevap ? interception - maxevap : 0.0
+    canopy_drainage = interception > maxevap ? interception - maxevap : 0.0
     interception = min(interception, maxevap)
 
-    # Add surpluss to the thoughdfall
-    throughfall = throughfall + overestimate
+    # Add surpluss to the throughfall
+    throughfall = throughfall + canopy_drainage
 
     return throughfall, interception, stemflow, canopystorage
 
@@ -82,37 +82,38 @@ function rainfall_interception_modrut(
     cmax,
 )
 
-    p = canopygapfraction
-    pt = 0.1 * p
+    # TODO: improve computation of stemflow partitioning coefficient pt (0.1 * canopygapfraction)
+    pt = min(0.1 * canopygapfraction, 1.0 - canopygapfraction)
 
     # Amount of p that falls on the canopy
-    pfrac = max((1.0 - p - pt), 0.0) * precipitation
+    precip_canopy = (1.0 - canopygapfraction - pt) * precipitation
 
-    # S cannot be larger than Cmax, no gravity drainage below that
-    dd = canopystorage > cmax ? canopystorage - cmax : 0.0
-    canopystorage = canopystorage - dd
+    # Canopystorage cannot be larger than cmax, no gravity drainage below that. This check
+    # is required because cmax can change over time
+    canopy_drainage1 = canopystorage > cmax ? canopystorage - cmax : 0.0
+    canopystorage = canopystorage - canopy_drainage1
 
     # Add the precipitation that falls on the canopy to the store
-    canopystorage = canopystorage + pfrac
+    canopystorage = canopystorage + precip_canopy
 
     # Now do the Evap, make sure the store does not get negative
-    dc = -1.0 * min(canopystorage, potential_evaporation)
-    canopystorage = canopystorage + dc
+    canopy_evap = min(canopystorage, potential_evaporation)
+    canopystorage = canopystorage - canopy_evap
 
-    leftover = potential_evaporation + dc
     # Amount of evap not used
+    leftover = potential_evaporation - canopy_evap
 
-    # Now drain the canopy storage again if needed...
-    d = canopystorage > cmax ? canopystorage - cmax : 0.0
-    canopystorage = canopystorage - d
+    # Now drain the canopystorage again if needed...
+    canopy_drainage2 = canopystorage > cmax ? canopystorage - cmax : 0.0
+    canopystorage = canopystorage - canopy_drainage2
 
-    # Calculate throughfall
-    throughfall = dd + d + p * precipitation
+    # Calculate throughfall and stemflow
+    throughfall = canopy_drainage1 + canopy_drainage2 + canopygapfraction * precipitation
     stemflow = precipitation * pt
 
     # Calculate interception, this is NET Interception
-    netinterception = precipitation - throughfall - stemflow
-    interception = -dc
+    netinterception = precipitation + canopy_drainage1 - throughfall - stemflow
+    interception = canopy_evap
 
     return netinterception, throughfall, stemflow, leftover, interception, canopystorage
 
@@ -153,7 +154,7 @@ function acttransp_unsat_sbm(
     θₛ,
     θᵣ,
     hb,
-    ust::Bool = false,
+    ust::Bool=false,
 )
 
     # AvailCap is fraction of unsat zone containing roots
@@ -362,9 +363,9 @@ function snowpack_hbv(
     ttm,
     cfmax,
     whc;
-    rfcf = 1.0,
-    sfcf = 1.0,
-    cfr = 0.05,
+    rfcf=1.0,
+    sfcf=1.0,
+    cfr=0.05
 )
 
     # fraction of precipitation which falls as rain