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)

Bases: SimPEG.simulation.BaseTimeSimulation

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)

property root_finder

Root-finding Algorithm

fields(m=None)

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

dpred(m, f=None)

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)

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)

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)

Jvec(m, v, f=None)

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)

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

alias of SimPEG.flow.richards.simulation.SimulationNDCellCentered

class SimPEG.flow.richards.simulation.RichardsProblem(*args, **kwargs)

Bases: SimPEG.flow.richards.simulation.SimulationNDCellCentered

This class has been deprecated, see SimulationNDCellCentered for documentation

Richards Survey

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

Bases: SimPEG.survey.BaseSurvey

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)

The Derivative with respect to the model.

derivAdjoint(simulation, f, v=None)

The adjoint derivative with respect to the model.

property rxList

rxList has been deprecated. See receiver_list for documentation

Richards receivers

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

Bases: SimPEG.survey.BaseTimeRx

Richards Receiver Object

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)

class SimPEG.flow.richards.receivers.Saturation(*args, **kwargs)

Bases: SimPEG.survey.BaseTimeRx

Richards Receiver Object

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)

Vadose Zone Empirical Relationships

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

Bases: SimPEG.props.HasModel

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)

Bases: SimPEG.flow.richards.empirical.NonLinearModel

Optional Properties:

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

plot(ax=None)

class SimPEG.flow.richards.empirical.BaseHydraulicConductivity(*args, **kwargs)

Bases: SimPEG.flow.richards.empirical.NonLinearModel

Optional Properties:

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

plot(ax=None)

class SimPEG.flow.richards.empirical.Haverkamp_theta(*args, **kwargs)

Bases: SimPEG.flow.richards.empirical.BaseWaterRetention

Optional Properties:

• alpha (PhysicalProperty): , a physical property, Default: 1611000.0

• alphaMap (Mapping): Mapping of to the inversion model., a SimPEG Map

• beta (PhysicalProperty): , a physical property, Default: 3.96

• betaMap (Mapping): Mapping of to the inversion model., a SimPEG Map

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

• theta_r (PhysicalProperty): residual water content [L3L-3], a physical property, Default: 0.075

• theta_rMap (Mapping): Mapping of residual water content [L3L-3] to the inversion model., a SimPEG Map

• theta_s (PhysicalProperty): saturated water content [L3L-3], a physical property, Default: 0.287

• theta_sMap (Mapping): Mapping of saturated water content [L3L-3] to the inversion model., a SimPEG Map

Other Properties:

• alphaDeriv (Derivative): Derivative of wrt the model.

• betaDeriv (Derivative): Derivative of wrt the model.

• theta_rDeriv (Derivative): Derivative of residual water content [L3L-3] wrt the model.

• theta_sDeriv (Derivative): Derivative of saturated water content [L3L-3] wrt the model.

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)

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)

class SimPEG.flow.richards.empirical.Haverkamp_k(*args, **kwargs)

Bases: SimPEG.flow.richards.empirical.BaseHydraulicConductivity

Optional Properties:

• A (PhysicalProperty): fitting parameter, a physical property, Default: 1175000.0

• AMap (Mapping): Mapping of fitting parameter to the inversion model., a SimPEG Map

• Ks (PhysicalProperty): Saturated hydraulic conductivity, a physical property, Default: 0.00944

• KsMap (Mapping): Mapping of Saturated hydraulic conductivity to the inversion model., a SimPEG Map

• gamma (PhysicalProperty): fitting parameter, a physical property, Default: 4.74

• gammaMap (Mapping): Mapping of fitting parameter to the inversion model., a SimPEG Map

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

Other Properties:

• ADeriv (Derivative): Derivative of fitting parameter wrt the model.

• KsDeriv (Derivative): Derivative of Saturated hydraulic conductivity wrt the model.

• gammaDeriv (Derivative): Derivative of fitting parameter wrt the model.

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)

derivM(u)

SimPEG.flow.richards.empirical.haverkamp(mesh, **kwargs)

class SimPEG.flow.richards.empirical.HaverkampParams

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)

Bases: SimPEG.flow.richards.empirical.BaseWaterRetention

Optional Properties:

• alpha (PhysicalProperty): related to the inverse of the air entry suction [L-1], >0, a physical property, Default: 0.036

• alphaMap (Mapping): Mapping of related to the inverse of the air entry suction [L-1], >0 to the inversion model., a SimPEG Map

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

• n (PhysicalProperty): measure of the pore-size distribution, >1, a physical property, Default: 1.56

• nMap (Mapping): Mapping of measure of the pore-size distribution, >1 to the inversion model., a SimPEG Map

• theta_r (PhysicalProperty): residual water content [L3L-3], a physical property, Default: 0.078

• theta_rMap (Mapping): Mapping of residual water content [L3L-3] to the inversion model., a SimPEG Map

• theta_s (PhysicalProperty): saturated water content [L3L-3], a physical property, Default: 0.43

• theta_sMap (Mapping): Mapping of saturated water content [L3L-3] to the inversion model., a SimPEG Map

Other Properties:

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

• nDeriv (Derivative): Derivative of measure of the pore-size distribution, >1 wrt the model.

• theta_rDeriv (Derivative): Derivative of residual water content [L3L-3] wrt the model.

• theta_sDeriv (Derivative): Derivative of saturated water content [L3L-3] wrt the model.

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)

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)

class SimPEG.flow.richards.empirical.Vangenuchten_k(*args, **kwargs)

Bases: SimPEG.flow.richards.empirical.BaseHydraulicConductivity

Optional Properties:

• I (PhysicalProperty): , a physical property, Default: 0.5

• IMap (Mapping): Mapping of to the inversion model., a SimPEG Map

• Ks (PhysicalProperty): Saturated hydraulic conductivity, a physical property, Default: 24.96

• KsMap (Mapping): Mapping of Saturated hydraulic conductivity to the inversion model., a SimPEG Map

• alpha (PhysicalProperty): related to the inverse of the air entry suction [L-1], >0, a physical property, Default: 0.036

• alphaMap (Mapping): Mapping of related to the inverse of the air entry suction [L-1], >0 to the inversion model., a SimPEG Map

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

• n (PhysicalProperty): measure of the pore-size distribution, >1, a physical property, Default: 1.56

• nMap (Mapping): Mapping of measure of the pore-size distribution, >1 to the inversion model., a SimPEG Map

Other Properties:

• IDeriv (Derivative): Derivative of wrt the model.

• KsDeriv (Derivative): Derivative of Saturated hydraulic conductivity wrt the model.

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

• nDeriv (Derivative): Derivative of measure of the pore-size distribution, >1 wrt the model.

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)

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)

SimPEG.flow.richards.empirical.van_genuchten(mesh, **kwargs)

class SimPEG.flow.richards.empirical.VanGenuchtenParams

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. [1982]

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