.. _api_Richards:


Richards Equation
*****************

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]

.. math::

    \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:

.. math::

    \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):

.. plot::

    from SimPEG import *
    from simpegFLOW import Richards
    import matplotlib.pyplot as plt

    M = Mesh.TensorMesh([np.ones(40)])
    M.setCellGradBC('dirichlet')
    params = Richards.Empirical.HaverkampParams().celia1990
    model = Richards.Empirical.Haverkamp(M, **params)

    bc = np.array([-61.5,-20.7])
    h = np.zeros(M.nC) + bc[0]


    def getFields(timeStep,method):
        prob = Richards.RichardsProblem(M,model, timeStep=timeStep, timeEnd=360,
                                        boundaryConditions=bc, initialConditions=h,
                                        doNewton=False, method=method)
        return prob.fields(params['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')

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


Richards
========

.. automodule:: simpegFLOW.Richards.Empirical
    :show-inheritance:
    :members:
    :undoc-members: