Fork me on GitHub

geoscience / scientific computing / visualization / web design / thoughts

notebooks / Richards Equation - Comparison to Ceilia et al. 1990

by Rowan Cockett

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
Warning: mumps solver not available.

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')
NewtonRoot stopped by maxIters. norm: 5.1942e-05
NewtonRoot stopped by maxIters. norm: 1.0531e-06
NewtonRoot stopped by maxIters. norm: 1.1545e-06
NewtonRoot stopped by maxIters. norm: 2.2226e-06
NewtonRoot stopped by maxIters. norm: 2.6211e-06

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]:
<matplotlib.legend.Legend at 0x10fb352d0>
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]: