Make it easier to replace the functional form of the hydraulic conductivity

[HannoSpreeuw]

Request from @EmiliaJarochowska

Make it straightforward to change its functional form, now it is completely hidden.

[EmiliaJarochowska]

This was dictated by the idea that the expression used for $K$ currently may need to be substituted by a different one. Currently it's split between the expressions presum and rhorat. Based on Hsu & Cheng (1991), L'Heureux used two different expressions for the cases on either side of $\phi = 0.95$. But that seems to be a misunderstunding of the paper.
As Jourabchi (2007) explains in chapter 5:

So Eq 15 includes already the limit case and is supposed to cover the entire interval of interest $K \in [0.5, 1]$.

Why therefore the different formulations in Fortran, which introduce additional instability? Something to ask Ivan about. He was a co-author of this chapter by Jourabchi (2007) so he should have known.

Does it mean this issue is obsolete? It is needed for #24 and it's not difficult to change. But it means that the equation we wanted to substitute current formula for K with is not needed.

[HannoSpreeuw]

I guess there are two options for introducing the function K as an argument to the rhs calculations in both evolution_rate and _make_pde_rhs_numba:

1. Something similar to the call to calculate_sigma from within evolution_rate (for _make_pde_rhs_numba this is calculated within, so no equivalent call there) . So one could define the function K in the LHeureux_model.py module.
2. Make the function defining K an argument of the LMAHeureuxPorosityDiff class in which case the definition of K would be added to a module different from LHeureux_model.py. Perhaps to parameters.py.

In both cases the implementations would differ between evolution_rate and _make_pde_rhs_numba.

[HannoSpreeuw]

I guess the second option is the most elegant and flexible.
So that implies defining the function K within parameters.py. I think that makes sense since the user may decide how to define K. So one could think of it as a parameter chosen by the user.

Now how do we propagate K into evolution_rate and _make_pde_rhs_numba?
Well , the same way as all the other parameters.
But evolution_rate needs a version of K that returns K for all depths, while _make_pde_rhs_numba will need a version of K that returns K for a specific depth, since there is a loop over all depths in _make_pde_rhs_numba.

This can likely be accommodated for by first defining a function of K that returns K for a specific depth. The jitted version of that function K should be propagated when the Numba backend is used, into _make_pde_rhs_numba.
When the Numpy backend is used the vectorized version of that function K should be propagated into evolution_rate.

[HannoSpreeuw]

K only enters the differential equations through the velocities, i.e. through U and W, right?
I.e. only through equations 46 and 47 from L'Heureux?

[HannoSpreeuw]

@EmiliaJarochowska I removed compressibility from the title of this issue, since our discussion moved towards hydraulic conductivity. You explained that compressibility and conductivity are not directly connected, so please open a separate issue Make it easier to replace compressibility if you feel the need.

With that settled, do you agree that fixing this issue also fixes #21?

[HannoSpreeuw]

This issue turns out to be more involved than it seemed at first sight.
The problem is equation 43 from L'Heureux, which involves a spatial derivative of W.
Since W depends on K, the user must not only define K, but also its spatial derivative, as some expression times the spatial derivative of the porosity.

More accurately, the user has to supply the spatial derivative of K * (1 - Phi)**2 / Phi. See equation 47.

[EmiliaJarochowska]

Yes, I agree it would fix it. Apologies for making such redundant issues 🤦‍♀️

[HannoSpreeuw]

No problem.
Please check the Define_K_in_parameters.py branch before I submit a PR.

[EmiliaJarochowska]

This is very useful, although our investigations showed now that the model is not reliable for reasons other than the choice of hydraulic conductivity. I would still recommend including it in main.