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

$\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 Celia et al. (1990): ## Richards Simulation¶

class SimPEG.flow.richards.simulation.SimulationNDCellCentered(*args, **kwargs)[source]

Richards Simulation

Required Properties:

• boundary_conditions (Array): boundary conditions, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

• counter (Counter): A SimPEG.utils.Counter object, an instance of Counter

• debug (Boolean): Show all messages, a boolean, Default: False

• do_newton (Boolean): Do a Newton iteration vs. a Picard iteration, a boolean, Default: False

• hydraulic_conductivity (BaseHydraulicConductivity): hydraulic conductivity function, an instance of BaseHydraulicConductivity

• initial_conditions (Array): initial conditions, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

• mesh (BaseMesh): a discretize mesh instance, an instance of BaseMesh

• method (StringChoice): Formulation used, See notes in Celia et al., 1990, either “mixed” or “head”, Default: mixed

• root_finder_max_iter (Integer): Maximum iterations for root_finder iteration, an integer, Default: 30

• root_finder_tol (Float): tolerance of the root_finder, a float, Default: 0.0001

• sensitivity_path (String): path to store the sensitivty, a unicode string, Default: ./sensitivity/

• None

• solver_opts (Dictionary): solver options as a kwarg dict, a dictionary

• survey (BaseSurvey): a survey object, an instance of BaseSurvey

• t0 (Float): Origin of the time discretization, a float, Default: 0.0

• time_steps (TimeStepArray):

Sets/gets the time steps for the time domain simulation. You can set as an array of dt’s or as a list of tuples/floats. Tuples must be length two with […, (dt, repeat), …] For example, the following setters are the same:

sim.time_steps = [(1e-6, 3), 1e-5, (1e-4, 2)]
sim.time_steps = np.r_[1e-6,1e-6,1e-6,1e-5,1e-4,1e-4]


, an array or list of tuples specifying the mesh tensor of <class ‘float’> with shape (*)

• water_retention (BaseWaterRetention): water retention curve, an instance of BaseWaterRetention

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*) or (*, *)

property hydraulic_conductivity

hydraulic_conductivity (BaseHydraulicConductivity): hydraulic conductivity function, an instance of BaseHydraulicConductivity

property water_retention

water_retention (BaseWaterRetention): water retention curve, an instance of BaseWaterRetention

property boundary_conditions

boundary_conditions (Array): boundary conditions, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property initial_conditions

initial_conditions (Array): initial conditions, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

property debug

debug (Boolean): Show all messages, a boolean, Default: False

property method

method (StringChoice): Formulation used, See notes in Celia et al., 1990, either “mixed” or “head”, Default: mixed

property do_newton

do_newton (Boolean): Do a Newton iteration vs. a Picard iteration, a boolean, Default: False

property root_finder_max_iter

root_finder_max_iter (Integer): Maximum iterations for root_finder iteration, an integer, Default: 30

property root_finder_tol

root_finder_tol (Float): tolerance of the root_finder, a float, Default: 0.0001

getBoundaryConditions(ii, u_ii)[source]
property root_finder

Root-finding Algorithm

fields(m=None)[source]

u = fields(m) The field given the model. :param numpy.ndarray m: model :rtype: numpy.ndarray :return: u, the fields

dpred(m, f=None)[source]

Create the projected data from a model. The field, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!).

$d_\text{pred} = P(f(m), m)$

Where P is a projection of the fields onto the data space.

property Dz
diagsJacobian(m, hn, hn1, dt, bc)[source]

Diagonals and rhs of the jacobian system

The matrix that we are computing has the form:

.-                                      -. .-  -.   .-  -.
|  Adiag                                 | | h1 |   | b1 |
|   Asub    Adiag                        | | h2 |   | b2 |
|            Asub    Adiag               | | h3 | = | b3 |
|                 ...     ...            | | .. |   | .. |
|                         Asub    Adiag  | | hn |   | bn |
'-                                      -' '-  -'   '-  -'

getResidual(m, hn, h, dt, bc, return_g=True)[source]

Used by the root finder when going between timesteps

Where h is the proposed value for the next time iterate (h_{n+1})

Jfull(m=None, f=None)[source]
Jvec(m, v, f=None)[source]

Jv = Jvec(m, v, f=None) Effect of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: Jv

Jtvec(m, v, f=None)[source]

Jtv = Jtvec(m, v, f=None) Effect of transpose of J(m) on a vector v. :param numpy.ndarray m: model :param numpy.ndarray v: vector to multiply :param Fields f: fields :rtype: numpy.ndarray :return: JTv

SimPEG.flow.richards.simulation.SimulationNDCellCentred[source]
class SimPEG.flow.richards.simulation.RichardsProblem(*args, **kwargs)[source]

This class has been deprecated, see SimulationNDCellCentered for documentation

## Richards Survey¶

class SimPEG.flow.richards.survey.Survey(*args, **kwargs)[source]

RichardsSurvey

Required Properties:

• counter (Counter): A SimPEG counter object, an instance of Counter

• receiver_list (a list of BaseRx): list of receivers for flow simulations, a list (each item is an instance of BaseRx)

• source_list (a list of BaseSrc): A list of sources for the survey, a list (each item is an instance of BaseSrc)

property receiver_list

receiver_list (a list of BaseRx): list of receivers for flow simulations, a list (each item is an instance of BaseRx)

property nD

Number of data

deriv(simulation, f, du_dm_v=None, v=None)[source]

The Derivative with respect to the model.

derivAdjoint(simulation, f, v=None)[source]

The adjoint derivative with respect to the model.

property rxList

rxList has been deprecated. See receiver_list for documentation

class SimPEG.flow.richards.receivers.Pressure(*args, **kwargs)[source]

Required Properties:

• locations (RxLocationArray): Locations of the receivers (nRx x nDim), an array of receiver locations of <class ‘float’>, <class ‘int’> with shape (*, *)

• projGLoc (StringChoice): Projection grid location, default is CC, any of “CC”, “Fx”, “Fy”, “Fz”, “Ex”, “Ey”, “Ez”, “N”, Default: CC

• projTLoc (StringChoice): location on the time mesh where the data are projected from, either “N” or “CC”, Default: N

• storeProjections (Boolean): Store calls to getP (organized by mesh), a boolean, Default: True

• times (Array): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

deriv(U, simulation, du_dm_v=None, v=None, adjoint=False)[source]
class SimPEG.flow.richards.receivers.Saturation(*args, **kwargs)[source]

Required Properties:

• locations (RxLocationArray): Locations of the receivers (nRx x nDim), an array of receiver locations of <class ‘float’>, <class ‘int’> with shape (*, *)

• projGLoc (StringChoice): Projection grid location, default is CC, any of “CC”, “Fx”, “Fy”, “Fz”, “Ex”, “Ey”, “Ez”, “N”, Default: CC

• projTLoc (StringChoice): location on the time mesh where the data are projected from, either “N” or “CC”, Default: N

• storeProjections (Boolean): Store calls to getP (organized by mesh), a boolean, Default: True

• times (Array): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

deriv(U, simulation, du_dm_v=None, v=None, adjoint=False)[source]

class SimPEG.flow.richards.empirical.NonLinearModel(*args, **kwargs)[source]

A non linear model that has dependence on the fields and a model

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*) or (*, *)

counter = None

A SimPEG.utils.Counter object

mesh = None

A SimPEG Mesh

property nP

Number of parameters in the model.

class SimPEG.flow.richards.empirical.BaseWaterRetention(*args, **kwargs)[source]

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*) or (*, *)

plot(ax=None)[source]
class SimPEG.flow.richards.empirical.BaseHydraulicConductivity(*args, **kwargs)[source]

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*) or (*, *)

plot(ax=None)[source]
class SimPEG.flow.richards.empirical.Haverkamp_theta(*args, **kwargs)[source]

Optional Properties:

Other Properties:

property theta_r

residual water content [L3L-3]

property theta_rMap

Mapping of residual water content [L3L-3] to the inversion model.

property theta_rDeriv

Derivative of residual water content [L3L-3] wrt the model.

property theta_s

saturated water content [L3L-3]

property theta_sMap

Mapping of saturated water content [L3L-3] to the inversion model.

property theta_sDeriv

Derivative of saturated water content [L3L-3] wrt the model.

property alpha
property alphaMap

Mapping of to the inversion model.

property alphaDeriv

Derivative of wrt the model.

property beta
property betaMap

Mapping of to the inversion model.

property betaDeriv

Derivative of wrt the model.

derivM(u)[source]

derivative with respect to m

import sympy as sy

alpha, u, beta, theta_r, theta_s = sy.symbols(
'alpha u beta theta_r theta_s', real=True
)

f_n = (
alpha *
(theta_s - theta_r) /
(alpha + abs(u)**beta) +
theta_r
)

derivU(u)[source]
class SimPEG.flow.richards.empirical.Haverkamp_k(*args, **kwargs)[source]

Optional Properties:

Other Properties:

property Ks

Saturated hydraulic conductivity

property KsMap

Mapping of Saturated hydraulic conductivity to the inversion model.

property KsDeriv

Derivative of Saturated hydraulic conductivity wrt the model.

property A

fitting parameter

property AMap

Mapping of fitting parameter to the inversion model.

property ADeriv

Derivative of fitting parameter wrt the model.

property gamma

fitting parameter

property gammaMap

Mapping of fitting parameter to the inversion model.

property gammaDeriv

Derivative of fitting parameter wrt the model.

derivU(u)[source]
derivM(u)[source]
SimPEG.flow.richards.empirical.haverkamp(mesh, **kwargs)[source]
class SimPEG.flow.richards.empirical.HaverkampParams[source]

Bases: object

Holds some default parameterizations for the Haverkamp model.

property celia1990

Parameters used in:

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.

class SimPEG.flow.richards.empirical.Vangenuchten_theta(*args, **kwargs)[source]

Optional Properties:

Other Properties:

property theta_r

residual water content [L3L-3]

property theta_rMap

Mapping of residual water content [L3L-3] to the inversion model.

property theta_rDeriv

Derivative of residual water content [L3L-3] wrt the model.

property theta_s

saturated water content [L3L-3]

property theta_sMap

Mapping of saturated water content [L3L-3] to the inversion model.

property theta_sDeriv

Derivative of saturated water content [L3L-3] wrt the model.

property n

measure of the pore-size distribution, >1

property nMap

Mapping of measure of the pore-size distribution, >1 to the inversion model.

property nDeriv

Derivative of measure of the pore-size distribution, >1 wrt the model.

property alpha

related to the inverse of the air entry suction [L-1], >0

property alphaMap

Mapping of related to the inverse of the air entry suction [L-1], >0 to the inversion model.

property alphaDeriv

Derivative of related to the inverse of the air entry suction [L-1], >0 wrt the model.

derivM(u)[source]

derivative with respect to m

import sympy as sy

alpha, u, n, I, Ks, theta_r, theta_s = sy.symbols(
'alpha u n I Ks theta_r theta_s', real=True
)

m = 1.0 - 1.0/n
theta_e = 1.0 / ((1.0 + sy.functions.Abs(alpha * u) ** n) ** m)

f_n = (
(
theta_s - theta_r
) /
(
(1.0 + abs(alpha * u)**n) ** (1.0 - 1.0 / n)
) +
theta_r
)

derivU(u)[source]
class SimPEG.flow.richards.empirical.Vangenuchten_k(*args, **kwargs)[source]

Optional Properties:

Other Properties:

property Ks

Saturated hydraulic conductivity

property KsMap

Mapping of Saturated hydraulic conductivity to the inversion model.

property KsDeriv

Derivative of Saturated hydraulic conductivity wrt the model.

property I
property IMap

Mapping of to the inversion model.

property IDeriv

Derivative of wrt the model.

property n

measure of the pore-size distribution, >1

property nMap

Mapping of measure of the pore-size distribution, >1 to the inversion model.

property nDeriv

Derivative of measure of the pore-size distribution, >1 wrt the model.

property alpha

related to the inverse of the air entry suction [L-1], >0

property alphaMap

Mapping of related to the inverse of the air entry suction [L-1], >0 to the inversion model.

property alphaDeriv

Derivative of related to the inverse of the air entry suction [L-1], >0 wrt the model.

derivM(u)[source]

derivative with respect to m

import sympy as sy

alpha, u, n, I, Ks, theta_r, theta_s = sy.symbols(
'alpha u n I Ks theta_r theta_s', real=True
)

m = 1.0 - 1.0/n
theta_e = 1.0 / ((1.0 + sy.functions.Abs(alpha * u) ** n) ** m)

f_n = Ks * theta_e ** I * (
(1.0 - (1.0 - theta_e ** (1.0 / m)) ** m) ** 2
)

f_n = (
(
theta_s - theta_r
) /
(
(1.0 + abs(alpha * u)**n) ** (1.0 - 1.0 / n)
) +
theta_r
)

derivU(u)[source]
SimPEG.flow.richards.empirical.van_genuchten(mesh, **kwargs)[source]
class SimPEG.flow.richards.empirical.VanGenuchtenParams[source]

Bases: object

The RETC code for quantifying the hydraulic functions of unsaturated soils, Van Genuchten, M Th, Leij, F J, Yates, S R

Table 3: Average values for selected soil water retention and hydraulic conductivity parameters for 11 major soil textural groups according to Rawls et al. 

property sand
property loamy_sand
property sandy_loam
property loam
property silt_loam
property sandy_clay_loam
property clay_loam
property silty_clay_loam
property sandy_clay
property silty_clay
property clay