FLOW: Richards: 1D: Celia1990ΒΆ

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 Celia1990

\[\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 Celia1990 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 Celia1990 demonstrating the head-based formulation and the mixed-formulation.

(Source code, png, hires.png, pdf)

../../_images/FLOW_Richards_1D_Celia1990-1.png
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
import matplotlib.pyplot as plt
import numpy as np

from SimPEG import Mesh, Maps
from SimPEG.FLOW import Richards


def run(plotIt=True):
    """
        FLOW: Richards: 1D: Celia1990
        =============================

        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 Celia1990_

        .. 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 Celia1990_ 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 Celia1990_ demonstrating the
        head-based formulation and the mixed-formulation.

        .. _Celia1990: http://www.webpages.uidaho.edu/ch/papers/Celia.pdf
    """

    M = Mesh.TensorMesh([np.ones(40)])
    M.setCellGradBC('dirichlet')
    params = Richards.Empirical.HaverkampParams().celia1990
    k_fun, theta_fun = Richards.Empirical.haverkamp(M, **params)
    k_fun.KsMap = Maps.IdentityMap(nP=M.nC)

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

    def getFields(timeStep, method):
        timeSteps = np.ones(int(360/timeStep))*timeStep
        prob = Richards.RichardsProblem(
            M,
            hydraulic_conductivity=k_fun,
            water_retention=theta_fun,
            boundary_conditions=bc, initial_conditions=h,
            do_newton=False, method=method
        )
        prob.timeSteps = timeSteps
        return prob.fields(params['Ks'] * np.ones(M.nC))

    Hs_M010 = getFields(10, 'mixed')
    Hs_M030 = getFields(30, 'mixed')
    Hs_M120 = getFields(120, 'mixed')
    Hs_H010 = getFields(10, 'head')
    Hs_H030 = getFields(30, 'head')
    Hs_H120 = getFields(120, 'head')

    if not plotIt:
        return
    plt.figure(figsize=(13, 5))
    plt.subplot(121)
    plt.plot(40-M.gridCC, Hs_M010[-1], 'b-')
    plt.plot(40-M.gridCC, Hs_M030[-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_H010[-1], 'b-')
    plt.plot(40-M.gridCC, Hs_H030[-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'
    ))

if __name__ == '__main__':
    run()
    plt.show()