**Objective:**

There are two different forms of Richards equation that differ on how they deal with the non-linearity in the time-stepping term.

The most fundamental form, referred to as the
'mixed'-form of Richards Equation [Celia *et al.*, 1990]

$$\frac{\partial \theta(\psi)}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega$$

where $\theta$ is water content, and $\psi$ is pressure head. This formulation of Richards equation is called the 'mixed'-form because the equation is parameterized in $\psi$ but the time-stepping is in terms of $\theta$.

As noted in [Celia *et al.*, 1990] the 'head'-based form of Richards
equation can be written in the continuous form as:

$$\frac{\partial \theta}{\partial \psi}\frac{\partial \psi}{\partial t} - \nabla \cdot k(\psi) \nabla \psi - \frac{\partial k(\psi)}{\partial z} = 0 \quad \psi \in \Omega$$

However, it can be shown that this does not conserve mass in the discrete formulation.

Here we reproduce the results from Ceilia *et al.* (1990) demonstrating that the head-based formulation is not as good as the mixed-formulation.

*Celia, Michael A., Efthimios T. Bouloutas, and Rebecca L. Zarba. "A general mass-conservative numerical solution for the unsaturated flow equation." Water Resources Research 26.7 (1990): 1483-1496.*

In [1]:

```
from SimPEG import *
import Richards
```

In [2]:

```
M = mesh.TensorMesh([np.ones(40)])
M.setCellGradBC('dirichlet')
Ks = 9.4400e-03
E = Richards.Haverkamp(Ks=np.log(Ks), A=1.1750e+06, gamma=4.74,
alpha=1.6110e+06, theta_s=0.287,
theta_r=0.075, beta=3.96)
bc = np.array([-61.5,-20.7])
h = np.zeros(M.nC) + bc[0]
def getFields(timeStep,method):
prob = Richards.RichardsProblem(M,E, timeStep=timeStep, timeEnd=360,
boundaryConditions=bc, initialConditions=h,
doNewton=False, method=method)
return prob.field(np.log(Ks))
Hs_M10 = getFields(10, 'mixed')
Hs_M30 = getFields(30, 'mixed')
Hs_M120= getFields(120,'mixed')
Hs_H10 = getFields(10, 'head')
Hs_H30 = getFields(30, 'head')
Hs_H120= getFields(120,'head')
```

In [3]:

```
figure(figsize=(13,5))
subplot(121)
plot(40-M.gridCC, Hs_M10[-1],'b-')
plot(40-M.gridCC, Hs_M30[-1],'r-')
plot(40-M.gridCC, Hs_M120[-1],'k-')
title('Mixed Method');xlabel('Depth, cm');ylabel('Pressure Head, cm')
legend(('dt = 10 sec','dt = 30 sec','dt = 120 sec'))
subplot(122)
plot(40-M.gridCC, Hs_H10[-1],'b-')
plot(40-M.gridCC, Hs_H30[-1],'r-')
plot(40-M.gridCC, Hs_H120[-1],'k-')
title('Head-Based Method');xlabel('Depth, cm');ylabel('Pressure Head, cm')
legend(('dt = 10 sec','dt = 30 sec','dt = 120 sec'))
```

Out[3]:

In [4]:

```
def function(var, ax, clim, tlt, i):
ax.cla()
linem, = ax.plot(40-M.gridCC, Hs_M10[i], 'r-', lw=1)
lineh, = ax.plot(40-M.gridCC, Hs_H10[i], 'b-', lw=1)
ax.legend(('mixed','head'))
ax.set_title('Time %4.f seconds'%(i*10))
M.video(Hs_H10,function,colorbar=False,figsize=(6,5))
```

Out[4]: