Description
This issue requires some scientific work first before turning it into code (and that also needs some thought). Copied from Teams:
In the Fortran code, hydraulic conductivity K is defined differently depending on the range of porosities (using a conditional) and Ivan uses equations from the paper by Hsu & Cheng (1991) for the different cases. Fortran has three cases:
a. the porosity is close to one (phi>0.95)
In this case Eq 16 from L'Heureux 2018 is used, which corresponds to Eq 5 in Hsu & Cheng, k(i)=betasV*10.*ph(i)**2/(1-ph(i))
b. the porosity is nearly 0 (phi<eps)
no permeability, only solids move
c. "normal case"
Eq 15 from L'Heureux 2018 is used (Eq 4 Hsu & Cheng, but with empirical factors gathered together into beta),
k(i)=betasV*ph(i)**3/(1-ph(i))**2
k(i)=k(i)*(1-dexp(-10.*(1-ph(i))/ph(i)))
According to this paper, we should be fine using the generalized version with the F factor, as python does currently:
Yet Ivan uses explicitly different expressions for different cases in the new version he sent in January. Perhaps because phi can fall below 0.4?
According to Boudreau's book, Ivan's choice of K is sketchy. About the Carman-Kozeny equation (Eq. 14) he writes:
None of these equations has been well established for very high porosities (ϕ>0.8), as found in muddy surficial sediments.
And specifically cites Hsu & Chang's empirical relation used by Ivan (Eqs 15 and 17, his 4.127) as:
This last equation has not been tested experimentally, and its applicability to natural sediments remains an open question.
The jump between cases a and c might cause instability. I have tested it using the Fortran code (without @jhidding's makefile) but I am not sure how to interpret the results.