Green-ampt infiltration not calculated correctly

[lrntct]

Original report by Anonymous.

---

When using the green-ampt infiltration model, the simulation results of infiltration seem to be dependent on the water-depth. Where water-depths are low, the infiltration rates are low. With fixed infiltration this does not happen.

From the publication of itzi the green-ampt Model is implemented as follow:

The infiltration rate is independent on water head in this implementation. Therefore the infiltration simulation output should be constant for the same green-ampt parameters where ponding occurs. Example case: By setting the suction head to 0, the infiltration rate should yield the hydraulic conductivity K according to the formula. The simulation was performed with a homogeneous rain event of 100 mm/h over a 4m x 4m domain with the suction head = 0mm and conductivity K = 10 mm/h. The result output of infiltration after 20min is as follow:

I would appreciate if you could investigate where the problem is coming from. I looked up the code and found maybe an error, not sure that this is the cause. It seems to me that when capping the infiltration rate, the time step dt is used instead of dt_h :

line 423 of flow.pyx : infrate = cap_inf_rate(dt, arr_h[r, c], infrate)

Thanks in advance for the investigation of this issue.

[lrntct]

Original comment by pirmin borer (Bitbucket: PBorer, ).

---

After changing dt to dt_h in the method in line 423 of flow.pyx and compiling with cython, the green-ampt infiltration yields correct results. However I performed a series of tests on a square 10x10m inclined plane with 5% slope and 0.2x0.2m resolution. The rainfall was set to 60 mm/h, conductivity to 5 mm/h, suction head SH to 150mm and effective porosity to 0.4.

In the first simulation I had no infiltration at all on the slope, only at the bottom of the slope of the domain where water accumulated. Infiltration was only computed where the hmin treshold was exceeded. On a flat domain with same parameters, the infiltration was correctly computed also below the treshold. It somehow seems to be linked with the fixed routing velocity. I realized that I had defined a hydrological timestep of 60 seconds and a routing velocity of 0.25 m/s which would yield 15m between each timestep. By reducing the hydrological timestep to 1 second and a velocity of 0.1 m/s, the infiltration was correctly computed compared to the green-ampt formula with same soil parameters (graph below). This indicates that the coupling timestep should be limited to some sort of CFL condition for the hydrology (e.g 0.7 * min(dx,dy)/v_routing). I couldn’t find a plausible explanation why the infiltration is correctly computed for waterheights above the hmin treshold by just looking at the code. Maybe the hydrological losses are applied at the timestep of the flow computation where hmin is exceeded. Maybe you could implement that the timestep for hydrology is calculated from the routing velocity and cellsize.

Thanks in advance for the code improvement and future enhancements.

‌

Pirmin Borer

[lrntct]

Original comment by Laurent Courty (Bitbucket: lrntct, GitHub: lrntct).

---

Thanks a lot for your report! The development of itzi started again after a hiatus, so we will be looking into the issue.

[lrntct]

Original comment by pirmin borer (Bitbucket: PBorer, ).

---

I made several improvements to the code:

• Corrected capping of infiltration and accounting for rain and 2d flow from neighbour cells.

• Added depression storage. Depression storage must be filed before runoff starts. Depression storage is recovered with infiltration.

• Added a parameter in config file to use a swmm raingage as precipitation input.

• Refactored the 1d coupling code to the integral form for stability improvement. Capping of outflow to 1d nodes accounts now for incoming flow from neighbour cells. 2d timestep now accounts for 1d node head. This prevents instability with large outflows from 1d nodes.

• Corrected an error of flow stencil in flow.pyx

• Added swim .out result file generation

• Modified courant condition to be cfl × min(dx,dy) / max(v+sqrt(g×h)). High velocity flows do not cause instability with this condition.

• Modified boundary conditions. 0 means now no flow boundary. Flow is not calculated where boundary type is zero. Therefore prevents calculation on masked areas and speeds up calculation. Boundary type 1 is now the compute domain. Boundary type 2 is free outflow boundary. Can also be used inside domain. Boundary type 5 is only flow, no hydrology applied (for swmm catchments where hydrology is already computed in 1d model).

All these modifications have been thoroughly tested. I can give you the improved code. Please contact me at febex2008@gmail.com

Further planned improvements:

• Porting to python 3
• Adding boundary condition for coupling 2d cells to a conduit (open channel flow)
• Adding a weir boundary condition
• Update the massbalance with boundary type 2.

Cheers.

Pirmin Borer

‌

[lrntct]

I worked on the infiltration module. The solution yields less than 1% error compared to the integral estimated with the Lambert W method. See https://github.com/lrntct/itzi/tree/bmi/itzi Indeed, the default time-step of 1min could induce missing the infiltration on slopes. Considering that the infiltration computation is rather cheap, it could be calculated at every routing time step.