One interesting fact that is not clear to most of people is that the chain derivation rule is not valid when we go from the continuous domain of calculus to the discrete domain of numerics. This means, for instance that:

$$frac{partial theta(psi(t))}{partial t} neq frac{partial theta(psi(t))}{partial psi} frac{partial psi(t)}{partial t} $$

where, for instance, $theta$ is the volumetric water contente, $psi$ is soil suction, $t$ is time.  The arcane is explained because when you go to discretize the differential we have:

$$ frac{theta^{n+1}-theta^n}{psi^{n+1}-psi^n} $$

where the upper scripts $n+1$ and $n$ indicate respectively the $n+1$ and $n$ time step.  Both on the numerator and denominator we have a dependence on time step $n+1$ which is unknown and the way they are approximated can bring to an inconsistent estimation of the derivative.  As a result,  conservation equations are said written in conservative form when they contain are written as:

$$ frac{partial # }{partial t}= – nabla cdot * $$

where $#$ stands for any state variable to be conserved and $*$ for any flux expression. $nabla cdot$ is the divergence operator. 

This issue was brought to light by the mathematician P.L. Roe in his 1981 paper, and to us by colleagues Vincenzo Casulli and Michael Dumbser and is clearly illustrated in Tubini et al. (2021) below.  

References

Roe, P. L.: Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics, 43, 357–372, 1981

Tubini, N., Gruber, S., and Rigon, R.: A method for solving heat transfer with phase change in ice or soil that allows for large time steps while guaranteeing energy conservation, The Cryosphere Discuss. [preprint], https://doi.org/10.5194/tc-2020-293, in review, 2020.