Richards equation

By applying the chain rule on temporal derivative leads to

${\displaystyle {\frac {\partial \theta (h)}{\partial t}}={\frac {{\textrm {d}}\theta }{{\textrm {d}}h}}{\frac {\partial h}{\partial t}}}$,

where ${\displaystyle {\frac {{\textrm {d}}\theta }{{\textrm {d}}h}}}$ is known as the retention water capacity ${\displaystyle C(h)}$. The equation is then stated as

${\displaystyle C(h){\frac {\partial h}{\partial t}}=\nabla \cdot \left(\mathbf {K} (h)\nabla h+\nabla z\right)-S}$.

The head-based Richards equation is prone to the following computational issue: the discretized temporal derivative using the implicit Rothe method yields the following approximation:

${\displaystyle {\frac {\Delta \theta }{\Delta t}}\approx C(h){\frac {\Delta h}{\Delta t}},\quad {\mbox{and so}}\quad {\frac {\Delta \theta }{\Delta t}}-C(h){\frac {\Delta h}{\Delta t}}=\varepsilon .}$

This approximation produces an error ${\displaystyle \varepsilon }$ that affects the mass conservation of the numerical solution, and so special strategies for temporal derivatives treatment are necessary.[7]

By applying the chain rule on the spatial derivative leads to

${\displaystyle \mathbf {K} (h)\nabla h=\mathbf {K} (h){\frac {{\textrm {d}}h}{{\textrm {d}}\theta }}\nabla \theta ,}$

where ${\displaystyle \mathbf {K} (h){\frac {{\textrm {d}}h}{{\textrm {d}}\theta }}}$, which could be further formulated as ${\displaystyle {\frac {\mathbf {K} (\theta )}{C(\theta )}}}$, is known as the soil water diffusivity ${\displaystyle \mathbf {D} (\theta )}$. The equation is then stated as

${\displaystyle {\frac {\partial \theta }{\partial t}}=\nabla \cdot \mathbf {D} (\theta )\nabla \theta -S.}$

The saturation-based Richards equation is prone to the following computational issues. Since the limits ${\displaystyle \lim _{\theta \to \theta _{s}}||\mathbf {D} (\theta )||=\infty }$ and ${\displaystyle \lim _{\theta \to \theta _{r}}||\mathbf {D} (\theta )||=\infty }$, where ${\displaystyle \theta _{s}}$ is the saturated (maximal) water content and ${\displaystyle \theta _{r}}$ is the residual (minimal) water content a successful numerical solution is restricted just for ranges of water content satisfactory below the full saturation (the saturation should be even lower than air entry value) as well as satisfactory above the residual water content.[8]