Minor coding error in infiltrationproblem Celia_MPM

[Huite]

Hi @amireson,

I've benchmarking some ODE solvers in Julia, and I found that your paper "A simple, efficient, mass-conservative approach to solving Richards' equation (openRE, v1.0)" provides a nice approach as well as a nice set of benchmarks.

In reproducing some of the results, I noticed that the Celia solution for the infiltration problem is overly diffusive. This is especially apparent for the drainage plot in figure 7 of the paper, but it's something I couldn't reproduce.

The issue appears to be a coding error in the residual function:

openRE/infiltrationproblem/modelbenchmarks/celia/celia_MPM.py

Lines 77 to 78 in bf095e0

	x3=-Kmid[1:]*(psigrid[2:]-psigrid[1:-1])
	x4=Kmid[:-1]*(psigrid[1:-1]-psigrid[:-2])

These terms are missing division by dz; e.g. your implementation of the Miller benchmark does contain them:

openRE/millerproblem/celiasolution/celia_MPM_T1_VG_psiall.py

Lines 74 to 75 in bf095e0

	x3=-Kmid[1:]*(psigrid[2:]-psigrid[1:-1])/dz
	x4=Kmid[:-1]*(psigrid[1:-1]-psigrid[:-2])/dz

Inserting /dz here, re-running the benchmark, shows a solution which is still a bit diffusive due to the large fixed time step, but which is much closer to the ATS and Hydrus solutions (we can now recognize the same individual bumps):

[amireson]

Hello @Huite - thank you very much for finding this and sharing with me. An unfortunate error that slipped in somehow. I have corrected this and the performance of the Celia model is improved. Here's the new figure:

I'm not sure if I can update this in the manuscript or not - I will contact the editor.

Good luck with your work with Julia. I continue to work with the scipy ODE solvers on increasingly complex problems - I'm very satisfied with their performance, especially when combined with numba and intelligent handling of the jacobian pattern. I'd be interested in see your work if and when it's ready to share.

Cheers!

Andrew