diff --git a/docs/src/changelog.md b/docs/src/changelog.md index de9f4cffb..e9560baab 100644 --- a/docs/src/changelog.md +++ b/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 diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md index 6456500e8..4c58258f0 100644 --- a/docs/src/model_docs/params_lateral.md +++ b/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 diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md index 12e52bcd6..e52431555 100644 --- a/docs/src/model_docs/params_vertical.md +++ b/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}`` | - | diff --git a/docs/src/user_guide/additional_options.md b/docs/src/user_guide/additional_options.md index 6b5baa03b..60e59cace 100644 --- a/docs/src/user_guide/additional_options.md +++ b/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 diff --git a/server/test/client.jl b/server/test/client.jl index 36e6b5dc5..f99c73dbf 100644 --- a/server/test/client.jl +++ b/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", diff --git a/src/flow.jl b/src/flow.jl index c778b9f6c..c30b1697b 100644 --- a/src/flow.jl +++ b/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,23 +234,63 @@ 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 @@ -251,7 +298,7 @@ function update(sf::SurfaceFlowLand, network, frac_toriver) 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, diff --git a/src/sbm.jl b/src/sbm.jl index baab14691..3c32b9500 100644 --- a/src/sbm.jl +++ b/src/sbm.jl @@ -154,6 +154,10 @@ infiltsoilpath::Vector{T} # Infiltration excess water [mm Δt⁻¹] infiltexcess::Vector{T} + # Infiltration from surface water [mm Δt⁻¹] + infilt_surfacewater::Vector{T} + # Contribution from surface water [mm Δt⁻¹] + contribution_surfacewater::Vector{T} # Water that cannot infiltrate due to saturated soil (saturation excess) [mm Δt⁻¹] excesswater::Vector{T} # Water exfiltrating during saturation excess conditions [mm Δt⁻¹] @@ -715,6 +719,8 @@ function initialize_sbm(nc, config, riverfrac, inds) ae_openw_l = fill(mv, n), ae_openw_r = fill(mv, n), avail_forinfilt = fill(mv, n), + infilt_surfacewater = fill(mv, n), + contribution_surfacewater = fill(mv, n), actinfilt = fill(mv, n), actinfiltsoil = fill(mv, n), actinfiltpath = fill(mv, n), @@ -869,6 +875,8 @@ function update_until_recharge(sbm::SBM, config) transfermethod = get(config.model, "transfermethod", false)::Bool ust = get(config.model, "whole_ust_available", false)::Bool # should be removed from optional setting and code? ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String + do_surface_water_infiltration = + get(config.model, "surface_water_infiltration", false)::Bool threaded_foreach(1:sbm.n, basesize = 250) do i if modelsnow @@ -901,7 +909,18 @@ function update_until_recharge(sbm::SBM, config) ustoredepth = sum(@view sbm.ustorelayerdepth[i][1:sbm.nlayers[i]]) runoff_river = min(1.0, sbm.riverfrac[i]) * avail_forinfilt - runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt + + # Calculate the initial capacity of the unsaturated store + ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth + # Add a flag if the soil is saturated + check_saturated = ustorecapacity < 0.1 + + if check_saturated + runoff_land = min(1.0, 1.0 - sbm.riverfrac[i]) * avail_forinfilt + else + runoff_land = min(1.0, sbm.waterfrac[i]) * avail_forinfilt + end + if !isnothing(sbm.paddy) || !isnothing(sbm.nonpaddy) avail_forinfilt = avail_forinfilt + sbm.allocation.irri_alloc[i] end @@ -911,18 +930,38 @@ function update_until_recharge(sbm::SBM, config) sbm.waterlevel_river[i] * sbm.riverfrac[i], sbm.riverfrac[i] * sbm.potential_evaporation[i], ) - ae_openw_l = min( - sbm.waterlevel_land[i] * sbm.waterfrac[i], - sbm.waterfrac[i] * sbm.potential_evaporation[i], - ) - - # evap available for soil evaporation + # Determine remaining fraction of available potential_evaporation soilevap_fraction = max( - sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.waterfrac[i] - - sbm.glacierfrac[i], + sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.glacierfrac[i], 0.0, ) - potsoilevap = soilevap_fraction * sbm.potential_evaporation[i] + # Determine amount of surface water + surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i]) + + # Calculate open water evaporation + if check_saturated + ae_openw_l = min( + surface_water, + sbm.potential_evaporation[i] * soilevap_fraction, + ) + else + ae_openw_l = min( + surface_water, + sbm.potential_evaporation[i] * sbm.waterfrac[i], + ) + end + + # Update land waterlevel and scale to part of cell not covered by rivers + if do_surface_water_infiltration + contribution_surfacewater = surface_water - ae_openw_l + # Add land waterlevel to infiltration + avail_forinfilt += contribution_surfacewater + else + contribution_surfacewater = 0.0 + end + + # evap available for soil evaporation + potsoilevap = (sbm.potential_evaporation[i] * soilevap_fraction) - ae_openw_l if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i] evap_paddy_water = min(sbm.paddy.h[i], potsoilevap) @@ -933,9 +972,6 @@ function update_until_recharge(sbm::SBM, config) evap_paddy_water = 0.0 end - # Calculate the initial capacity of the unsaturated store - ustorecapacity = sbm.soilwatercapacity[i] - sbm.satwaterdepth[i] - ustoredepth - # Calculate the infiltration flux into the soil column infiltsoilpath, infiltsoil, @@ -1121,7 +1157,17 @@ function update_until_recharge(sbm::SBM, config) end actinfilt = infiltsoilpath - du - excesswater = avail_forinfilt - infiltsoilpath - infiltexcess + du + + # Scale infiltration from surface water based on the ratio between actinfil and + # infiltsoilpath (and prevent division by zero) + if do_surface_water_infiltration + infilt_ratio = iszero(avail_forinfilt) ? 0.0 : actinfilt / avail_forinfilt + infilt_surfacewater = max(0.0, contribution_surfacewater * infilt_ratio) + else + infilt_surfacewater = 0.0 + end + + excesswater = avail_forinfilt - actinfilt - infiltexcess # Separation between compacted and non compacted areas (correction with the satflow du) # This is required for D-Emission/Delwaq @@ -1190,6 +1236,8 @@ function update_until_recharge(sbm::SBM, config) sbm.avail_forinfilt[i] = avail_forinfilt sbm.actinfilt[i] = actinfilt sbm.infiltexcess[i] = infiltexcess + sbm.infilt_surfacewater[i] = infilt_surfacewater + sbm.contribution_surfacewater[i] = contribution_surfacewater sbm.recharge[i] = recharge sbm.transpiration[i] = transpiration sbm.soilevap[i] = soilevap @@ -1223,7 +1271,7 @@ function update_until_recharge(sbm::SBM, config) end end -function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater) +function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater, config) threaded_foreach(1:sbm.n, basesize = 1000) do i usl, n_usl = set_layerthickness(zi[i], sbm.sumlayers[i], sbm.act_thickl[i]) @@ -1244,24 +1292,38 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater) ustoredepth = sum(@view usld[1:n_usl]) + # Correct the excess water to only include the water that is not (yet) surface water, + # otherwise this water would be accounted for twice. This is only relevant if + # do_surface_water_infiltration is true + do_surface_water_infiltration = + get(config.model, "surface_water_infiltration", false)::Bool + if do_surface_water_infiltration + correction_surfacewater = + iszero(sbm.avail_forinfilt[i]) ? 1.0 : + 1.0 - (sbm.contribution_surfacewater[i] / sbm.avail_forinfilt[i]) + else + correction_surfacewater = 1.0 + end + if !isnothing(sbm.paddy) && sbm.paddy.irrigation_areas[i] paddy_h_add = exfiltustore + exfiltsatwater[i] + - sbm.excesswater[i] + + sbm.excesswater[i] * correction_surfacewater + sbm.runoff_land[i] + - sbm.infiltexcess[i] + sbm.infiltexcess[i] * correction_surfacewater runoff = max(paddy_h_add - sbm.paddy.h_max[i], 0.0) sbm.paddy.h[i] = paddy_h_add - runoff else runoff = exfiltustore + exfiltsatwater[i] + - sbm.excesswater[i] + + sbm.excesswater[i] * correction_surfacewater + sbm.runoff_land[i] + - sbm.infiltexcess[i] + sbm.infiltexcess[i] * correction_surfacewater end + # volumetric water content per soil layer and root zone vwc = sbm.vwc[i] vwc_perc = sbm.vwc_perc[i] @@ -1305,7 +1367,7 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater) sbm.exfiltsatwater[i] = exfiltsatwater[i] sbm.exfiltustore[i] = exfiltustore sbm.runoff[i] = runoff - sbm.net_runoff[i] = runoff - sbm.ae_openw_l[i] + sbm.net_runoff[i] = runoff - sbm.ae_openw_l[i] - sbm.infilt_surfacewater[i] sbm.vwc[i] = vwc sbm.vwc_perc[i] = vwc_perc sbm.rootstore[i] = rootstore @@ -1429,13 +1491,13 @@ function update_water_demand(sbm::SBM) sbm.nonpaddy.maximum_irrigation_rate[i] if depletion >= raw # start irrigation irri_dem_gross += depletion - # add depletion to irrigation gross demand when the maximum irrigation rate has been + # add depletion to irrigation gross demand when the maximum irrigation rate has been # applied at the previous time step (to get volumetric water content at field capacity) elseif depletion > 0.0 && max_irri_rate_applied # continue irrigation irri_dem_gross += depletion end end - # limit irrigation demand to infiltration capacity + # limit irrigation demand to infiltration capacity infiltration_capacity = sbm.soilinfredu[i] * (1.0 - sbm.pathfrac[i]) * sbm.infiltcapsoil[i] irri_dem_gross = min(irri_dem_gross, infiltration_capacity) @@ -1466,7 +1528,7 @@ function update_water_demand(sbm::SBM) end sbm.paddy.demand_gross[i] = irri_dem_gross end - # update gross water demands + # update gross water demands sbm.allocation.irri_demand_gross[i] = irri_dem_gross sbm.allocation.nonirri_demand_gross[i] = industry_dem + domestic_dem + livestock_dem sbm.allocation.total_gross_demand[i] = diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl index 2dfe4a2ae..f313fc745 100644 --- a/src/sbm_gwf_model.jl +++ b/src/sbm_gwf_model.jl @@ -600,6 +600,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel} vertical, (network.land.altitude .- aquifer.head) .* 1000.0, # zi [mm] in vertical concept SBM exfiltwater .* 1000.0, + config, ) ssf_toriver = zeros(vertical.n) diff --git a/src/sbm_model.jl b/src/sbm_model.jl index ec1032894..c706fa03c 100644 --- a/src/sbm_model.jl +++ b/src/sbm_model.jl @@ -518,6 +518,7 @@ function update_after_subsurfaceflow( vertical, lateral.subsurface.zi * 1000.0, lateral.subsurface.exfiltwater * 1000.0, + config, ) ssf_toriver = lateral.subsurface.to_river ./ tosecond(basetimestep) diff --git a/test/bmi.jl b/test/bmi.jl index ba0f3751b..dc78246a1 100644 --- a/test/bmi.jl +++ b/test/bmi.jl @@ -20,8 +20,8 @@ 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) == 207 - @test BMI.get_output_item_count(model) == 207 + @test BMI.get_input_item_count(model) == 210 + @test BMI.get_output_item_count(model) == 210 to_check = [ "vertical.nlayers", "vertical.theta_r", diff --git a/test/run_sbm.jl b/test/run_sbm.jl index d50ccb752..f97cf9df2 100644 --- a/test/run_sbm.jl +++ b/test/run_sbm.jl @@ -540,3 +540,81 @@ Wflow.close_files(model, delete_output = false) Wflow.close_files(model, delete_output = false) end + + +# test infiltration of surface water +@testset "surface water infiltration" begin + i_surface = 498 # pixel with surface water, infiltration and excesswater + i_nosurface = 2444 # pixel without surface water + + tomlpath = joinpath(@__DIR__, "sbm_config.toml") + config = Wflow.Config(tomlpath) + config.model.surface_water_infiltration = true + + model = Wflow.initialize_sbm_model(config) + + model = Wflow.run_timestep(model) + model = Wflow.run_timestep(model) + + @unpack vertical = model + + # Test is surface water is indeed added to avail_forinfilt + @test vertical.waterlevel_land[i_surface] ≈ 16.191219691562683 + @test vertical.waterlevel_land[i_nosurface] == 0.0 + + @test vertical.avail_forinfilt[i_surface] ≈ 16.24809470532854 + @test vertical.avail_forinfilt[i_nosurface] ≈ 0.08956446991756856 + + # Test if the paritioning of infiltrated water is correct + @test vertical.infilt_surfacewater[i_surface] ≈ 9.1388544728681 + @test vertical.infilt_surfacewater[i_nosurface] == 0.0 + + # Test excesswater and runoff for two locations + @test vertical.excesswater[i_surface] ≈ 7.077138112678853 + @test vertical.runoff[i_surface] ≈ 88.75818647205926 + + # Test is net_runoff is indeed lower than runoff + @test vertical.net_runoff[i_surface] ≈ 79.61933199919116 + @test vertical.net_runoff[i_nosurface] == 0.0 + + Wflow.close_files(model, delete_output = false) +end + + +# test flow threshold for kinematic overland flow +@testset "surface water infiltration" begin + idx = 498 # pixel with surface water, infiltration and excesswater + threshold = 0.02 + + tomlpath = joinpath(@__DIR__, "sbm_config.toml") + config = Wflow.Config(tomlpath) + config.input.lateral.land.h_thresh = Dict("value" => threshold) + + mod_threshold = Wflow.initialize_sbm_model(config) + + + # run two timesteps for both models + mod_threshold = Wflow.run_timestep(mod_threshold) + mod_threshold = Wflow.run_timestep(mod_threshold) + Wflow.close_files(mod_threshold, delete_output = false) + + + config = Wflow.Config(tomlpath) + mod_nothreshold = Wflow.initialize_sbm_model(config) + + mod_nothreshold = Wflow.run_timestep(mod_nothreshold) + mod_nothreshold = Wflow.run_timestep(mod_nothreshold) + Wflow.close_files(mod_nothreshold, delete_output = false) + + + @unpack lateral = mod_threshold + lateral_threshold = lateral + @unpack lateral = mod_nothreshold + lateral_nothreshold = lateral + + @test lateral_threshold.land.q_av[idx] ≈ 1.3364541870336923 + @test lateral_threshold.land.q_av[idx] < lateral_nothreshold.land.q_av[idx] + @test lateral_threshold.land.pond_height[idx] ≈ 0.01591595571103741 + @test lateral_nothreshold.land.pond_height[idx] == 0.0 + +end