Reinfiltration of surface water and flow threshold for overland flow #470
Conversation
|
Nice work @JoostBuitink ! I will submit a review this Friday. One question I did already have: you did check the water balance of the vertical part (runoff), did you also check the water balance of the overland flow? I would recommend this as the code has changed quite a bit (I think mainly the kinematic wave). |
| 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`). |
There was a problem hiding this comment.
| 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`). | |
| overland flow. This parameter can be changed through the `lateral.land.h_thresh` key in the | |
| TOML file. The parameter can be added as a map or as a fixed value (through the use of | |
| `h_thresh.value`). |
| - 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 |
There was a problem hiding this comment.
| - 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 | |
| - Support for reinfiltration of surface overland flow water for both kinematic wave and local inertial routing schemes. Can be set through the `surface_water_infiltration` key in the `[model]` section of the TOML file. This way, the water stored on the |
| @@ -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: | |||
There was a problem hiding this comment.
| To allow for re-infiltration of overland flow, the following model flag needs to be set: | |
| To allow for re-infiltration of overland flow as part of the kinematic wave and local inertial routing schemes, the following model flag needs to be set: |
| surface_water_infiltration = true | ||
| ``` | ||
|
|
||
| When using this option, the average surface water level (of the previous time step) is added to |
There was a problem hiding this comment.
| When using this option, the average surface water level (of the previous time step) is added to | |
| When using this option, the average surface water level (at the previous time step) is added to |
| `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. |
There was a problem hiding this comment.
| `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. | |
| `excesswater` variable, to prevent the surface water from being counted twice as excess water. |
| 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. |
There was a problem hiding this comment.
Maybe remove this? I think it can be a bit confusing to readers (e.g. why is this then not the case for local inertial routing?
| 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: |
There was a problem hiding this comment.
| flow is allowed to occur. This flow threshold can be set as a map from the staticmaps: | |
| flow is allowed to occur. This flow threshold can be set as a static map as part of the input netCDF file: |
src/flow.jl
Outdated
| 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[v] |
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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.
src/flow.jl
Outdated
| end | ||
|
|
||
| # Update values | ||
| sf.pond_height[v] = pond_volume_current[v] / (sf.width[v] * sf.dl[v]) |
There was a problem hiding this comment.
I think first an update is required of pond_volume_current[v] before computing sf.pond_height[v]? When it is exceeding the threshold_volume[v].
|
|
||
|
|
||
| # test flow threshold for kinematic overland flow | ||
| @testset "surface water infiltration" begin |
There was a problem hiding this comment.
| @testset "surface water infiltration" begin | |
| @testset "flow threshold for kinematic overland flow" begin |
src/flow.jl
Outdated
| threshold_volume = sf.h_thresh .* sf.width .* sf.dl | ||
|
|
||
| # Get pond volume at the start of the time step | ||
| pond_volume_current = sf.pond_height .* sf.width .* sf.dl |
There was a problem hiding this comment.
This will result in more allocation and could have an impact on performance. Would be nice if this can be avoided (in the for loop?).
src/flow.jl
Outdated
| @@ -207,10 +207,10 @@ function update(sf::SurfaceFlowLand, network, frac_toriver) | |||
| @. sf.qlat = sf.inwater / sf.dl | |||
|
|
|||
| # Convert threshold to volume | |||
| threshold_volume = sf.h_thresh .* sf.width .* sf.dl | |||
| threshold_volume = sf.h_thresh .* area | |||
There was a problem hiding this comment.
Why the use of area here? It's confusing as it deviates from the overland flow dimensions width and dl. For river cells area is higher than sf.width .* sf.dl, for land cells it has the same value. I would suggest to use the overland flow dimensions for the computation of the threshold_volume, and the pond_volume_current.
What most probably is not correct yet when using the overland flow dimensions, is the amount of surface water that can infiltrate from river cells. I think this amount should be scaled by the land fraction of these cells. We have a similar approach for evaporation and direct runoff by using waterfrac.
src/sbm.jl
Outdated
|
|
||
| if do_surface_water_infiltration | ||
| # Add land waterlevel to infiltration | ||
| avail_forinfilt += waterlevel_land |
There was a problem hiding this comment.
The amount of surface water that can infiltrate needs to be scaled by the land fraction. We use a similar approach for runoff and evaporation from open water (using waterfrac).
In the wflow python code (version with support for re-infiltration of surface water) a similar approach is used (in this case for rivers): https://github.com/openstreams/wflow/blob/bdc9b7966a0caa3c826b3a1f25150186b37d683f/wflow/wflow_sbm_old.py#L1849.
| ) | ||
| # Determine amount of surface water | ||
| surface_water = sbm.waterlevel_land[i] * (1.0 - sbm.riverfrac[i]) |
There was a problem hiding this comment.
Would recommend to use the same approach for runoff_land (L. 912) for consistency. Here waterfrac is still used.
There was a problem hiding this comment.
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:
- Revert back to using
waterfrac, the approach is then consistent between surface runoff and open water evaporation, or - Use
waterfracunless 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).
Issue addressed
Fixes #237
Explanation
Explain how you addressed the bug/feature request, what choices you made and why.
Checklist
masterAdditional Notes (optional)
Add any additional notes or information that may be helpful.