From f02a8eb728887ad91b0cdcc1e4e313db703099d3 Mon Sep 17 00:00:00 2001 From: verseve Date: Fri, 1 Sep 2023 09:00:25 +0200 Subject: [PATCH 1/4] fix waterbalance rutter interception --- src/vertical_process.jl | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/vertical_process.jl b/src/vertical_process.jl index b9f073d23..0f246e935 100644 --- a/src/vertical_process.jl +++ b/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 @@ -82,13 +82,14 @@ 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 + pfrac = (1.0 - canopygapfraction - pt) * precipitation - # S cannot be larger than Cmax, no gravity drainage below that + # Canopystorage cannot be larger than cmax, no gravity drainage below that. This check + # is required because cmax can change over time dd = canopystorage > cmax ? canopystorage - cmax : 0.0 canopystorage = canopystorage - dd @@ -99,19 +100,19 @@ function rainfall_interception_modrut( dc = -1.0 * min(canopystorage, potential_evaporation) canopystorage = canopystorage + dc - leftover = potential_evaporation + dc # Amount of evap not used - - # Now drain the canopy storage again if needed... + leftover = potential_evaporation + dc + + # Now drain the canopystorage again if needed... d = canopystorage > cmax ? canopystorage - cmax : 0.0 canopystorage = canopystorage - d - # Calculate throughfall - throughfall = dd + d + p * precipitation + # Calculate throughfall and stemflow + throughfall = dd + d + canopygapfraction * precipitation stemflow = precipitation * pt # Calculate interception, this is NET Interception - netinterception = precipitation - throughfall - stemflow + netinterception = precipitation + dd - throughfall - stemflow interception = -dc return netinterception, throughfall, stemflow, leftover, interception, canopystorage From d02f9007ac5853517ca64ddbcd057c1d41bafe7d Mon Sep 17 00:00:00 2001 From: verseve Date: Fri, 1 Sep 2023 10:00:30 +0200 Subject: [PATCH 2/4] Update changelog --- docs/src/changelog.md | 10 ++++++++++ src/vertical_process.jl | 4 ++-- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index e7cfc5d17..659ce68ff 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -5,6 +5,16 @@ 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] + +### 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. + ## v0.7.1 - 2023-06-30 ### Fixed diff --git a/src/vertical_process.jl b/src/vertical_process.jl index 0f246e935..f18ec2d2a 100644 --- a/src/vertical_process.jl +++ b/src/vertical_process.jl @@ -61,7 +61,7 @@ function rainfall_interception_gash( overestimate = interception > maxevap ? interception - maxevap : 0.0 interception = min(interception, maxevap) - # Add surpluss to the thoughdfall + # Add surpluss to the throughfall throughfall = throughfall + overestimate return throughfall, interception, stemflow, canopystorage @@ -82,7 +82,7 @@ function rainfall_interception_modrut( cmax, ) - # TODO: improve computation of stemflow partitioning coefficient pt ((0.1 * canopygapfraction) + # 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 From c7c615ef001e643dc1f3d5055d85e276d57796a1 Mon Sep 17 00:00:00 2001 From: verseve Date: Fri, 1 Sep 2023 12:25:16 +0200 Subject: [PATCH 3/4] Interception from modified Rutter model - The `interception` output (interception loss by evaporation) from the modified Rutter interception model should be used in `SBM`, also as input for the computation of total actual evapotranspiration `actevap`. - Update of docs and changelog. --- docs/src/changelog.md | 4 ++++ docs/src/model_docs/params_vertical.md | 2 +- src/sbm.jl | 3 +-- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/src/changelog.md b/docs/src/changelog.md index 659ce68ff..5313f7194 100644 --- a/docs/src/changelog.md +++ b/docs/src/changelog.md @@ -14,6 +14,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 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 diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md index 56bac933c..de39ff770 100644 --- a/docs/src/model_docs/params_vertical.md +++ b/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}`` | - | diff --git a/src/sbm.jl b/src/sbm.jl index bd19aa5cc..00631753a 100644 --- a/src/sbm.jl +++ b/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 From 00361d240a3706320171faf0929c6fb4f0d40868 Mon Sep 17 00:00:00 2001 From: JoostBuitink <44062204+JoostBuitink@users.noreply.github.com> Date: Tue, 5 Sep 2023 10:57:33 +0200 Subject: [PATCH 4/4] improved variable names in interception functions --- src/vertical_process.jl | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/vertical_process.jl b/src/vertical_process.jl index f18ec2d2a..f3fe22c1d 100644 --- a/src/vertical_process.jl +++ b/src/vertical_process.jl @@ -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 throughfall - throughfall = throughfall + overestimate + throughfall = throughfall + canopy_drainage return throughfall, interception, stemflow, canopystorage @@ -86,34 +86,34 @@ function rainfall_interception_modrut( pt = min(0.1 * canopygapfraction, 1.0 - canopygapfraction) # Amount of p that falls on the canopy - pfrac = (1.0 - canopygapfraction - pt) * precipitation + precip_canopy = (1.0 - canopygapfraction - pt) * precipitation # Canopystorage cannot be larger than cmax, no gravity drainage below that. This check # is required because cmax can change over time - dd = canopystorage > cmax ? canopystorage - cmax : 0.0 - canopystorage = canopystorage - dd + 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 # Amount of evap not used - leftover = potential_evaporation + dc - + leftover = potential_evaporation - canopy_evap + # Now drain the canopystorage again if needed... - d = canopystorage > cmax ? canopystorage - cmax : 0.0 - canopystorage = canopystorage - d + canopy_drainage2 = canopystorage > cmax ? canopystorage - cmax : 0.0 + canopystorage = canopystorage - canopy_drainage2 # Calculate throughfall and stemflow - throughfall = dd + d + canopygapfraction * precipitation + throughfall = canopy_drainage1 + canopy_drainage2 + canopygapfraction * precipitation stemflow = precipitation * pt # Calculate interception, this is NET Interception - netinterception = precipitation + dd - throughfall - stemflow - interception = -dc + netinterception = precipitation + canopy_drainage1 - throughfall - stemflow + interception = canopy_evap return netinterception, throughfall, stemflow, leftover, interception, canopystorage @@ -154,7 +154,7 @@ function acttransp_unsat_sbm( θₛ, θᵣ, hb, - ust::Bool = false, + ust::Bool=false, ) # AvailCap is fraction of unsat zone containing roots @@ -363,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