Reinfiltration of surface water and flow threshold for overland flow

docs/src/changelog.md

@@ -5,6 +5,17 @@ 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
+
+### Added
+- Support for reinfiltration of surface overland flow water. Can be set through the
+  `surface_water_infiltration` key in the `[model]` section. This way, the water stored on the
+  land surface is added to the water that is allowed to infiltrate into the soil.
+- Support for adding a water level flow-threshold for both kinematic wave and local inertial
+  overland flow. This setting can be changed by changing the `lateral.land.h_thresh` key in the
+  settings file. This can be added as a map or as a fixed value (through the use of
+  `h_thresh.value`).
+
 ## v0.8.1 - 2024-08-27
 
 ### Fixed
@@ -35,7 +46,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   is fixed by using `divrem` for the computation of the number of `steps` in this function.
   An error is thrown when the absolute remainder of `divrem` is larger than `eps()`, or when
   the number of `steps` is negative.
- - Fixed internal and external broken links in docs. 
+ - Fixed internal and external broken links in docs.
  - The internal time step of the local inertial model (`stable_timestep` function) can get
   zero when `LoopVectorization` is applied (`@tturbo`) to the for loop of these functions.
   This issue occured on a virtual machine, Windows 10 Enterprise, with Intel(R) Xeon(R) Gold

docs/src/model_docs/params_lateral.md

@@ -74,6 +74,8 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `cel`           | celerity of kinematic wave | m s``^{-1}`` | - |
 | `to_river`      | part of overland flow that flows to the river | m``^3`` s``^{-1}`` | - |
 | `kinwave_it`    | boolean for kinematic wave iterations | - | false |
+| **`h_thresh`**  | threshold for water level before flow occurs | m | 0 |
+| `pond_height`   | waterlevel of the pond (water stored on land surface) | m | - |
 
 ### [Reservoirs](@id reservoir_params)
 The Table below shows the parameters (fields) of struct `SimpleReservoir`, including a

docs/src/model_docs/params_vertical.md

@@ -117,6 +117,7 @@ profile `kv` is used and `z_layered` is required as input.
 | `actinfiltpath` | actual infiltration into compacted fraction  | mm Δt``^{-1}`` | - |
 | `infiltsoilpath` | infiltration into the unsaturated zone | mm Δt``^{-1}`` | - |
 | `infiltexcess` | infiltration excess water | mm Δt``^{-1}`` | - |
+| `infilt_surfacewater` | part of infiltration that originates from surface water | mm Δt``^{-1}`` | - |
 | `excesswater` | water that cannot infiltrate due to saturated soil (saturation excess) | mm Δt``^{-1}`` | - |
 | `exfiltsatwater` | water exfiltrating during saturation excess conditions | mm Δt``^{-1}`` | - |
 | `exfiltustore` | water exfiltrating from unsaturated store because of change in water table | mm Δt``^{-1}`` | - |

docs/src/user_guide/additional_options.md

@@ -183,6 +183,38 @@ irrigation_trigger = "irrigation_trigger"
 h = "h_paddy"
 ```
 
+## Enabling re-infiltration of overland flow
+To allow for re-infiltration of overland flow, the following model flag needs to be set:
+
+```toml
+[model]
+surface_water_infiltration = true
+```
+
+When using this option, the average surface water level (of the previous time step) is added to
+the volume of water that is allowed to infiltrate (`avail_forinfilt`). To keep track of the
+infiltrated water that originates from the surface water, the variable `infilt_surfacewater` is
+used. This variable is then subtracted from the net runoff, such that the overland flow is
+corrected for the loss of this water. Lastly, some corrections are applied on the
+`excesswater`, to prevent the surface water from being counted twice as excess water. This
+option works for both kinematic wave and local inertial overland flow.
+
+## Enabling a water level threshold for overland flow
+A water level threshold can be set for the overland flow. As long as this threshold is not
+exceeded, flow is not allowed to occur, and the water will be considered a "pond" (the variable
+`lateral.land.pond_height` can be used to keep track of its water level for kinematic wave
+overland flow, and the variable `lateral.land.h` can be used). When this threshold is exceeded,
+flow is allowed to occur. This flow threshold can be set as a map from the staticmaps:
+
+```toml
+[input.lateral.land]
+h_thresh = "overland_flow_threshold"
+```
+
+It should be noted that (for kinematic wave overland flow) flow is allowed to occur when the
+`pond_height` plus the incoming fluxes would be greater than the `h_thresh`. As a result, flow
+can occur at pond heights that are slightly lower than the set threshold.
+
 ## [Using multithreading] (@id multi_threading)
 
 ### Using wflow in Julia

server/test/client.jl

@@ -32,8 +32,8 @@ end
 
 @testset "model information functions" begin
     @test request((fn = "get_component_name",)) == Dict("component_name" => "sbm")
-    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 207)
-    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 207)
+    @test request((fn = "get_input_item_count",)) == Dict("input_item_count" => 210)
+    @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 210)
     to_check = [
         "vertical.nlayers",
         "vertical.theta_r",

src/flow.jl

@@ -63,6 +63,8 @@ end
     cel::Vector{T} | "m s-1"                        # Celerity of the kinematic wave
     to_river::Vector{T} | "m3 s-1"                  # Part of overland flow [m³ s⁻¹] that flows to the river
     kinwave_it::Bool | "-" | 0 | "none" | "none"    # Boolean for iterations kinematic wave
+    h_thresh::Vector{T} | "m"                       # Threshold for water level before flow occurs
+    pond_height::Vector{T} | "m"                    # Waterlevel of the pond (water stored on land surface)
 end
 
 function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, tstep, dt)
@@ -75,6 +77,9 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
         ncread(nc, config, "lateral.land.n"; sel = inds, defaults = 0.072, type = Float)
     n = length(inds)
 
+    h_thresh =
+        ncread(nc, config, "lateral.land.h_thresh"; sel = inds, defaults = 0, type = Float)
+
     sf_land = SurfaceFlowLand(
         beta = Float(0.6),
         sl = sl,
@@ -97,6 +102,8 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
         cel = zeros(Float, n),
         to_river = zeros(Float, n),
         kinwave_it = iterate,
+        h_thresh = h_thresh,
+        pond_height = zeros(Float, n),
     )
 
     return sf_land
@@ -227,31 +234,71 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         )
                     end
 
-                    sf.q[v] = kinematic_wave(
-                        sf.qin[v],
-                        sf.q[v],
-                        sf.qlat[v],
-                        sf.alpha[v],
-                        sf.beta,
-                        dt,
-                        sf.dl[v],
+                    # Convert threshold to volume
+                    threshold_volume = sf.h_thresh[v] * sf.width[v] * sf.dl[v]
+                    pond_volume_current = sf.pond_height[v] * sf.width[v] * sf.dl[v]
+
+                    # Calculate potential inflow
+                    pot_inflow = (sf.qlat[v] * sf.dl[v] * dt) + (sf.qin[v] * dt)
+
+                    # Check potential pond volume, to see if flow is allowed to occur
+                    pond_volume_pot = max(
+                        pond_volume_current + pot_inflow,
+                        0.0,
                     )
 
-                    # update h, only if surface width > 0.0
-                    if sf.width[v] > 0.0
-                        crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
-                        sf.h[v] = crossarea / sf.width[v]
+                    # Start kinematic wave if pond volume exceeds threshold
+                    if pond_volume_pot >= threshold_volume
+                        # Calculate fraction of pond volume that is above threshold
+                        if pot_inflow <= 0.0
+                            flowing_fraction = 1.0
+                        else
+                            flowing_fraction = (pond_volume_pot - threshold_volume) / pot_inflow
+                            flowing_fraction = min(flowing_fraction, 1.0)
+                        end
+
+                        sf.q[v] = kinematic_wave(
+                            sf.qin[v] * flowing_fraction,
+                            sf.q[v],
+                            sf.qlat[v] * flowing_fraction,
+                            sf.alpha[v],
+                            sf.beta,
+                            dt,
+                            sf.dl[v],
+                        )
+
+                        # update h, only if surface width > 0.0
+                        if sf.width[v] > 0.0
+                            crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
+                            sf.h[v] = crossarea / sf.width[v]
+                        end
+
+                        # Update pond volume
+                        pond_volume_current = threshold_volume
+
+                    else
+                        # No flow if pond volume is below threshold
+                        sf.q[v] = 0.0
+                        sf.h[v] = 0.0
+                        # Update pond volume
+                        pond_volume_current = pond_volume_pot
                     end
+
+                    # Update values
+                    sf.pond_height[v] = pond_volume_current / (sf.width[v] * sf.dl[v])
+
+                    # Update average values
                     sf.q_av[v] += sf.q[v]
-                    sf.h_av[v] += sf.h[v]
+                    sf.h_av[v] += (sf.h[v] + sf.pond_height[v])
+
                 end
             end
         end
     end
     sf.q_av ./= its
     sf.h_av ./= its
     sf.to_river ./= its
-    sf.volume .= sf.dl .* sf.width .* sf.h
+    sf.volume .= sf.dl .* sf.width .* (sf.h .+ sf.pond_height)
 end
 
 function update(sf::SurfaceFlowRiver, network, doy)
@@ -1015,7 +1062,7 @@ const dirs = (:yd, :xd, :xu, :yu)
     g::T | "m2 s-1" | 0 | "scalar"                          # acceleration due to gravity
     theta::T | "-" | 0 | "scalar"                           # weighting factor (de Almeida et al., 2012)
     alpha::T | "-" | 0 | "scalar"                           # stability coefficient (de Almeida et al., 2012)
-    h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
+    h_thresh::Vector{T} | "m" | 0 | "scalar"                # depth threshold for calculating flow
     dt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     qy0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in y direction at previous time step
     qx0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in x direction at previous time step
@@ -1054,9 +1101,17 @@ function initialize_shallowwater_land(
     froude_limit = get(config.model, "froude_limit", true)::Bool # limit flow to subcritical according to Froude number
     alpha = get(config.model, "inertial_flow_alpha", 0.7)::Float64 # stability coefficient for model time step (0.2-0.7)
     theta = get(config.model, "inertial_flow_theta", 0.8)::Float64 # weighting factor
-    h_thresh = get(config.model, "h_thresh", 1.0e-03)::Float64 # depth threshold for flow at link
+    # depth threshold for flow at link
+    h_thresh = ncread(
+        nc,
+        config,
+        "lateral.land.h_thresh";
+        sel = inds,
+        defaults = 1.0e-03,
+        type = Float,
+    )
 
-    @info "Local inertial approach is used for overlandflow." alpha theta h_thresh froude_limit
+    @info "Local inertial approach is used for overlandflow." alpha theta froude_limit
 
     n_land =
         ncread(nc, config, "lateral.land.n"; sel = inds, defaults = 0.072, type = Float)
@@ -1257,7 +1312,7 @@ function shallowwater_update(
             zs_max = max(zs_x, zs_xu)
             hf = (zs_max - sw.zx_max[i])
 
-            if hf > sw.h_thresh
+            if hf > sw.h_thresh[i]
                 length = T(0.5) * (sw.xl[i] + sw.xl[xu]) # can be precalculated
                 sw.qx[i] = local_inertial_flow(
                     sw.theta,
@@ -1297,7 +1352,7 @@ function shallowwater_update(
             zs_max = max(zs_y, zs_yu)
             hf = (zs_max - sw.zy_max[i])
 
-            if hf > sw.h_thresh
+            if hf > sw.h_thresh[i]
                 length = T(0.5) * (sw.yl[i] + sw.yl[yu]) # can be precalculated
                 sw.qy[i] = local_inertial_flow(
                     sw.theta,