Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Reinfiltration of surface water and flow threshold for overland flow #470

Draft
wants to merge 28 commits into
base: master
Choose a base branch
from
Draft

Reinfiltration of surface water and flow threshold for overland flow #470

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
df51bea
Update sbm.jl
JoostBuitink Sep 3, 2024
6a984ce
Merge branch 'master' into reinfiltration
JoostBuitink Sep 3, 2024
ebe92bd
add support for h_thres for overland flow
JoostBuitink Sep 9, 2024
30e30fb
Merge branch 'master' into reinfiltration
JoostBuitink Sep 9, 2024
6bd036c
update test and changelog
JoostBuitink Sep 9, 2024
de4c229
clean up overland flow threshold
JoostBuitink Sep 17, 2024
6b70f1b
clean code related to `excesswater`
JoostBuitink Sep 17, 2024
68aed28
improve infilt_ratio
JoostBuitink Sep 17, 2024
127b68b
fix `excesswater`
JoostBuitink Sep 17, 2024
7f2d86c
fix water balance issue; improve readability of excesswater computation
JoostBuitink Sep 25, 2024
6e1fd92
move correction to runoff calculation
JoostBuitink Sep 26, 2024
79ed50b
add flow threshold for overland flow
JoostBuitink Sep 26, 2024
de29798
fix `h_thresh` type in struct definition
JoostBuitink Sep 26, 2024
2cd9577
update changelog
JoostBuitink Sep 26, 2024
f6bfe8a
update docs, fix gwf error
JoostBuitink Sep 27, 2024
ae37c4a
add tests
JoostBuitink Sep 27, 2024
8560072
update docs
JoostBuitink Sep 27, 2024
e5067f7
update WflowServer tests
JoostBuitink Sep 27, 2024
715485e
Merge branch 'master' into reinfiltration
JoostBuitink Sep 27, 2024
db41548
wrap lines
JoostBuitink Sep 27, 2024
b5d4889
Update run_sbm.jl
JoostBuitink Sep 27, 2024
56e7d32
fix computation of pond volume to cell area
JoostBuitink Oct 11, 2024
4a39e4a
fix such that only water above the threshold is allowed to flow
JoostBuitink Oct 11, 2024
01bd103
revert area calculation back to width*dl
JoostBuitink Oct 16, 2024
26f5004
remove dependency of waterfrac; fix flowing_fraction
JoostBuitink Oct 18, 2024
dd921b7
fix open water evaporation
JoostBuitink Oct 24, 2024
d05f673
move pond calculations to fully inside for loop
JoostBuitink Oct 24, 2024
2dfa845
fix flowing_fraction and volume calculation
JoostBuitink Oct 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions src/flow.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -252,8 +252,11 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
# Start kinematic wave if pond volume exceeds threshold
if pond_volume_pot >= threshold_volume[v]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this correct? For example, if pond_volume_current[v] is below the threshold_volume, and the input of water causes pond_volume_pot > threshold_volume, then part of the input of water to the kinematic wave should go first to the pond, and the rest is available for overland flow?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a quick thought: an alternative is to pull the pond out of the kinematic wave solution, and only consider vertical fluxes, similar to the Paddy approach.

# Calculate fraction of pond volume that is above threshold
flowing_fraction = iszero(pot_inflow) ? 1.0 :
(pond_volume_pot - threshold_volume[v]) / pot_inflow
if pot_inflow <= 0.0
flowing_fraction = 1.0
else
flowing_fraction = (pond_volume_pot - threshold_volume[v]) / pot_inflow
end

sf.q[v] = kinematic_wave(
sf.qin[v] * flowing_fraction,
Expand Down
38 changes: 23 additions & 15 deletions src/sbm.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,8 @@
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⁻¹]
Expand Down Expand Up @@ -718,6 +720,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
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),
Expand Down Expand Up @@ -916,26 +919,31 @@ 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],
# Determine remaining fraction of available potential_evaporation
soilevap_fraction = max(
sbm.canopygapfraction[i] - sbm.riverfrac[i] - sbm.glacierfrac[i],
0.0,
)
# Determine amount of surface water
surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would recommend to use the same approach for runoff_land (L. 912) for consistency. Here waterfrac is still used.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This implementation (not using waterfrac) results in a hydrograph that is quite flashy (tested by @JoostBuitink ). Most likely this is caused by a disconnect between the overland flow kinematic wave and soil: water is available in the kinematic wave while the soil (unsaturated zone) has also capacity left (ustorecapacitiy), especially in the case when infiltration of surface water is not allowed.

So, while we model the kinematic wave routing for the complete land fraction, and allow for open water evaporation and surface runoff for that fraction (by not using waterfrac), most likely the soil is not completely saturated. As a consequence, this results most likely in an "overestimation" of surface runoff (runoff_land) and open water evaporation.

Based on this, my recommendation would be to:

  1. Revert back to using waterfrac, the approach is then consistent between surface runoff and open water evaporation, or
  2. Use waterfrac unless the soil in the grid cell is almost completely saturated (ustorecapacitiy $\approx$ 0.0), then use the land fraction.

Another approach (for the longterm) and similar to 2) could be to use subgrid variability, for example by using a pdf of soil storage capacity or of topographic information to derive a dynamic saturated fraction (e.g. see also this paper).


# Update land waterlevel
waterlevel_land = sbm.waterlevel_land[i] - ae_openw_l
# Calculate open water evaporation
ae_openw_l = min(
surface_water,
sbm.potential_evaporation[i] * soilevap_fraction,
)

# 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 += waterlevel_land
avail_forinfilt += contribution_surfacewater
else
contribution_surfacewater = 0.0
end

# 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]
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)
Expand Down Expand Up @@ -1139,7 +1147,7 @@ function update_until_recharge(sbm::SBM, config)
# 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, waterlevel_land * infilt_ratio)
infilt_surfacewater = max(0.0, contribution_surfacewater * infilt_ratio)
else
infilt_surfacewater = 0.0
end
Expand Down Expand Up @@ -1214,6 +1222,7 @@ function update_until_recharge(sbm::SBM, config)
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
Expand Down Expand Up @@ -1274,10 +1283,9 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater, config)
do_surface_water_infiltration =
get(config.model, "surface_water_infiltration", false)::Bool
if do_surface_water_infiltration
contribution_surfacewater = sbm.waterlevel_land[i] - sbm.ae_openw_l[i]
correction_surfacewater =
iszero(sbm.avail_forinfilt[i]) ? 1.0 :
1.0 - (contribution_surfacewater / sbm.avail_forinfilt[i])
1.0 - (sbm.contribution_surfacewater[i] / sbm.avail_forinfilt[i])
else
correction_surfacewater = 1.0
end
Expand Down
Loading