# 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.

  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()