The RICHARDS mode applies to single phase, variably saturated, isothermal systems. The governing mass conservation equation is given by

(1)\[\frac{{{\partial}}}{{{\partial}}t}\left(\varphi s\eta\right) + {\boldsymbol{\nabla}}\cdot\left(\eta{\boldsymbol{q}}\right) = Q_w,\]

with Darcy flux \({\boldsymbol{q}}\) defined as

(2)\[{\boldsymbol{q}}= -\frac{kk_r(s)}{\mu}\left({\boldsymbol{\nabla}}P-W_w\rho {\boldsymbol{g}}\right).\]

Here, \(\varphi\) denotes porosity [-], \(s\) saturation [m\(^3\) m\(^{-3}\)], \(\eta\) molar water density [kmol m\(^{-3}\)], \(\rho\) mass water density [kg m\(^{-3}\)], \({\boldsymbol{q}}\) Darcy flux [m s\(^{-1}\)], \(k\) intrinsic permeability [m\(^2\)], \(k_r\) relative permeability [-], \(\mu\) viscosity [Pa s], \(P\) pressure [Pa], \(W_w\) formula weight of water [kg kmol\(^{-1}\)], and \({\boldsymbol{g}}\) gravity [m s\(^{-2}\)]. Supported relative permeability functions \(k_r\) for Richards’ equation include van Genuchten, Books-Corey and Thomeer-Corey, while the saturation functions include Burdine and Mualem. Water density and viscosity are computed as a function of temperature and pressure through an equation of state for water. The source/sink term \(Q_w\) [kmol m\(^{-3}\) s\(^{-1}\)] has the form

(3)\[Q_w = \frac{q_M}{W_w} \delta({\boldsymbol{r}}-{\boldsymbol{r}}_{ss}),\]

where \(q_M\) denotes a mass rate in kg/m\(^{3}\)/s, and \({\boldsymbol{r}}_{ss}\) denotes the location of the source/sink.

Capillary Pressure - Saturation Functions

van Genuchten Saturation Function

Capillary pressure is related to saturation by various phenomenological relations, one of which is the van Genuchten (1980) relation

(4)\[s_e = \left[1+\left( \frac{p_c}{p_c^0} \right)^n\right]^{-m},\]

where \(p_c\) represents the capillary pressure [Pa], and the effective saturation \(s_e\) is defined by

(5)\[s_e = \frac{s - s_r}{s_0 - s_r},\]

where \(s_r\) denotes the residual saturation, and \(s_0\) denotes the maximum saturation. The inverse relation is given by

(6)\[p_c = p_c^0 \left(s_e^{-1/m}-1\right)^{1/n}.\]

The quantities \(m\), \(n\) and \(p_c^0\) are impirical constants determined by fitting to experimental data.

Brooks-Corey Saturation Function

The Brooks-Corey saturation function is a limiting form of the van Genuchten relation for \(p_c/p_c^0 \gg 1\), with the form

(7)\[s_e = \left(\frac{p_c}{p_c^0}\right)^{-\lambda},\]

with \(\lambda=mn\) and inverse relation

(8)\[p_c = p_c^0 s_e^{-1/\lambda}.\]

Relative Permeability Functions

Two forms of the relative permeability function are implemented based on the Mualem and Burdine formulations. The quantity \(n\) is related to \(m\) by the expression

(9)\[m = 1-\frac{1}{n}, \ \ \ \ \ n = \frac{1}{1-m},\]

for the Mualem formulation and by

(10)\[m = 1-\frac{2}{n}, \ \ \ \ \ n = \frac{2}{1-m},\]

for the Burdine formulation.

Mualem Relative Permeability

For the Mualem relative permeability function based on the van Genuchten saturation function is given by the expression

(11)\[k_{r} = \sqrt{s_e} \left\{1 - \left[1- \left( s_e \right)^{1/m} \right]^m \right\}^2.\]

The Mualem relative permeability function based on the Brooks-Corey saturation function is defined by

(12)\[\begin{split}k_r &= \big(s_e\big)^{5/2+2/\lambda} \\ &=\big(p_c/p_c^0\big)^{-(5\lambda/2+2)}.\end{split}\]

Burdine Relative Permeability

For the Burdine relative permeability function based on the van Genuchten saturation function is given by the expression

(13)\[k_{r} = s_e^2 \left\{1 - \left[1- \left( s_e \right)^{1/m} \right]^m \right\}.\]

The Burdine relative permeability function based on the Brooks-Corey saturation function has the form

(14)\[\begin{split}k_r &= \big(s_e\big)^{3+2/\lambda} \\ &= \left(\frac{p_c}{p_c^0}\right)^{-(3+2\lambda)}.\end{split}\]


At the end points of the saturation and relative permeability functions it is sometimes necessary to smooth the functions in order for the Newton-Raphson equations to converge. This is accomplished using a third order polynomial interpolation by matching the values of the function to be fit (capillary pressure or relative permeability), and imposing zero slope at the fully saturated end point and matching the derivative at a chosen variably saturated point that is close to fully saturated. The resulting equations for coefficients \(a_i\), \(i=0-3\), are given by

(15)\[\begin{split}a_0 + a_1 x_1 + a_2 x_1^2 + a_3 x_1^3 &= f_1,\\ a_0 + a_1 x_2 + a_2 x_2^2 + a_3 x_2^3 &= f_2,\\ a_1 x_1 + 2a_2 x_1 + 3a_3 x_1^2 &= f_1',\\ a_1 x_2 + 2a_2 x_2 + 3a_3 x_2^2 &= f_2',\end{split}\]

for chosen points \(x_1\) and \(x_2\). In matrix form these equations become

(16)\[\begin{split}\begin{bmatrix} 1 & x_1 & x_1^2 & x_1^3\\ 1 & x_2 & x_2^2 & x_2^3\\ 0 & 1 & 2x_1 & 3x_1^2\\ 0 & 1 & 2x_2 & 3x_2^2 \end{bmatrix} \begin{bmatrix} a_0\\ a_1\\ a_2\\ a_3 \end{bmatrix} = \begin{bmatrix} f_1\\ f_2\\ f_1'\\ f_2' \end{bmatrix}.\end{split}\]

The conditions imposed on the smoothing equations for capillary pressure \(f=s_e(p_c)\) are \(x_1=2 p_c^0\), \(x_2=p_c^0/2\), \(f_1 = (s_e)_1\), \(f_2 = 1\), \(f_1' = (s_e')_1\), \(f_2' = 0\). For relative permeability \(f=k_r(s_e)\), \(x_1 = 1\), \(x_2 = 0.99\), \(f_1 = 1\), \(f_2 = (k_r)_2\), \(f_1' = 0\), \(f_2' = (k_r')_2\).