Richards Equation¶
There are two different forms of Richards equation that differ on how they deal with the nonlinearity in the timestepping term.
The most fundamental form, referred to as the ‘mixed’form of Richards Equation [Celia et al., 1990]
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 timestepping 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:
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]¶ 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 Counterdebug (
Boolean
): Show all messages, a boolean, Default: Falsedo_newton (
Boolean
): Do a Newton iteration vs. a Picard iteration, a boolean, Default: Falsehydraulic_conductivity (
BaseHydraulicConductivity
): hydraulic conductivity function, an instance of BaseHydraulicConductivityinitial_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 BaseMeshmethod (
StringChoice
): Formulation used, See notes in Celia et al., 1990, either “mixed” or “head”, Default: mixedroot_finder_max_iter (
Integer
): Maximum iterations for root_finder iteration, an integer, Default: 30root_finder_tol (
Float
): tolerance of the root_finder, a float, Default: 0.0001sensitivity_path (
String
): path to store the sensitivty, a unicode string, Default: ./sensitivity/None
solver_opts (
Dictionary
): solver options as a kwarg dict, a dictionarysurvey (
BaseSurvey
): a survey object, an instance of BaseSurveyt0 (
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 = [(1e6, 3), 1e5, (1e4, 2)] sim.time_steps = np.r_[1e6,1e6,1e6,1e5,1e4,1e4]
, an array or list of tuples specifying the mesh tensor of <class ‘float’> with shape (*)
 time_steps (
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
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

property
root_finder
¶ Rootfinding 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})

SimPEG.flow.richards.simulation.
SimulationNDCellCentred
[source]¶ alias of
SimPEG.flow.richards.simulation.SimulationNDCellCentered

class
SimPEG.flow.richards.simulation.
RichardsProblem
(*args, **kwargs)[source]¶ 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)[source]¶ Bases:
SimPEG.survey.BaseSurvey
RichardsSurvey
Required Properties:
counter (
Counter
): A SimPEG counter object, an instance of Counterreceiver_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

property
rxList
¶ rxList has been deprecated. See receiver_list for documentation
Richards receivers¶

class
SimPEG.flow.richards.receivers.
Pressure
(*args, **kwargs)[source]¶ 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: CCprojTLoc (
StringChoice
): location on the time mesh where the data are projected from, either “N” or “CC”, Default: NstoreProjections (
Boolean
): Store calls to getP (organized by mesh), a boolean, Default: Truetimes (
Array
): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)

class
SimPEG.flow.richards.receivers.
Saturation
(*args, **kwargs)[source]¶ 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: CCprojTLoc (
StringChoice
): location on the time mesh where the data are projected from, either “N” or “CC”, Default: NstoreProjections (
Boolean
): Store calls to getP (organized by mesh), a boolean, Default: Truetimes (
Array
): times where the recievers measure data, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)
Vadose Zone Empirical Relationships¶

class
SimPEG.flow.richards.empirical.
NonLinearModel
(*args, **kwargs)[source]¶ 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)[source]¶ Bases:
SimPEG.flow.richards.empirical.NonLinearModel
Optional Properties:
model (
Model
): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)

class
SimPEG.flow.richards.empirical.
BaseHydraulicConductivity
(*args, **kwargs)[source]¶ Bases:
SimPEG.flow.richards.empirical.NonLinearModel
Optional Properties:
model (
Model
): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)

class
SimPEG.flow.richards.empirical.
Haverkamp_theta
(*args, **kwargs)[source]¶ Bases:
SimPEG.flow.richards.empirical.BaseWaterRetention
Optional Properties:
alpha (
PhysicalProperty
): , a physical property, Default: 1611000.0alphaMap (
Mapping
): Mapping of to the inversion model., a SimPEG Mapbeta (
PhysicalProperty
): , a physical property, Default: 3.96betaMap (
Mapping
): Mapping of to the inversion model., a SimPEG Mapmodel (
Model
): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)theta_r (
PhysicalProperty
): residual water content [L3L3], a physical property, Default: 0.075theta_rMap (
Mapping
): Mapping of residual water content [L3L3] to the inversion model., a SimPEG Maptheta_s (
PhysicalProperty
): saturated water content [L3L3], a physical property, Default: 0.287theta_sMap (
Mapping
): Mapping of saturated water content [L3L3] 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 [L3L3] wrt the model.theta_sDeriv (
Derivative
): Derivative of saturated water content [L3L3] wrt the model.

property
theta_r
¶ residual water content [L3L3]

property
theta_rMap
¶ Mapping of residual water content [L3L3] to the inversion model.

property
theta_rDeriv
¶ Derivative of residual water content [L3L3] wrt the model.

property
theta_s
¶ saturated water content [L3L3]

property
theta_sMap
¶ Mapping of saturated water content [L3L3] to the inversion model.

property
theta_sDeriv
¶ Derivative of saturated water content [L3L3] 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.

class
SimPEG.flow.richards.empirical.
Haverkamp_k
(*args, **kwargs)[source]¶ Bases:
SimPEG.flow.richards.empirical.BaseHydraulicConductivity
Optional Properties:
A (
PhysicalProperty
): fitting parameter, a physical property, Default: 1175000.0AMap (
Mapping
): Mapping of fitting parameter to the inversion model., a SimPEG MapKs (
PhysicalProperty
): Saturated hydraulic conductivity, a physical property, Default: 0.00944KsMap (
Mapping
): Mapping of Saturated hydraulic conductivity to the inversion model., a SimPEG Mapgamma (
PhysicalProperty
): fitting parameter, a physical property, Default: 4.74gammaMap (
Mapping
): Mapping of fitting parameter to the inversion model., a SimPEG Mapmodel (
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.

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 massconservative numerical solution for the unsaturated flow equation.” Water Resources Research 26.7 (1990): 14831496.

property

class
SimPEG.flow.richards.empirical.
Vangenuchten_theta
(*args, **kwargs)[source]¶ Bases:
SimPEG.flow.richards.empirical.BaseWaterRetention
Optional Properties:
alpha (
PhysicalProperty
): related to the inverse of the air entry suction [L1], >0, a physical property, Default: 0.036alphaMap (
Mapping
): Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion model., a SimPEG Mapmodel (
Model
): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)n (
PhysicalProperty
): measure of the poresize distribution, >1, a physical property, Default: 1.56nMap (
Mapping
): Mapping of measure of the poresize distribution, >1 to the inversion model., a SimPEG Maptheta_r (
PhysicalProperty
): residual water content [L3L3], a physical property, Default: 0.078theta_rMap (
Mapping
): Mapping of residual water content [L3L3] to the inversion model., a SimPEG Maptheta_s (
PhysicalProperty
): saturated water content [L3L3], a physical property, Default: 0.43theta_sMap (
Mapping
): Mapping of saturated water content [L3L3] to the inversion model., a SimPEG Map
Other Properties:
alphaDeriv (
Derivative
): Derivative of related to the inverse of the air entry suction [L1], >0 wrt the model.nDeriv (
Derivative
): Derivative of measure of the poresize distribution, >1 wrt the model.theta_rDeriv (
Derivative
): Derivative of residual water content [L3L3] wrt the model.theta_sDeriv (
Derivative
): Derivative of saturated water content [L3L3] wrt the model.

property
theta_r
¶ residual water content [L3L3]

property
theta_rMap
¶ Mapping of residual water content [L3L3] to the inversion model.

property
theta_rDeriv
¶ Derivative of residual water content [L3L3] wrt the model.

property
theta_s
¶ saturated water content [L3L3]

property
theta_sMap
¶ Mapping of saturated water content [L3L3] to the inversion model.

property
theta_sDeriv
¶ Derivative of saturated water content [L3L3] wrt the model.

property
n
¶ measure of the poresize distribution, >1

property
nMap
¶ Mapping of measure of the poresize distribution, >1 to the inversion model.

property
nDeriv
¶ Derivative of measure of the poresize distribution, >1 wrt the model.

property
alpha
¶ related to the inverse of the air entry suction [L1], >0

property
alphaMap
¶ Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion model.

property
alphaDeriv
¶ Derivative of related to the inverse of the air entry suction [L1], >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 )

class
SimPEG.flow.richards.empirical.
Vangenuchten_k
(*args, **kwargs)[source]¶ Bases:
SimPEG.flow.richards.empirical.BaseHydraulicConductivity
Optional Properties:
I (
PhysicalProperty
): , a physical property, Default: 0.5IMap (
Mapping
): Mapping of to the inversion model., a SimPEG MapKs (
PhysicalProperty
): Saturated hydraulic conductivity, a physical property, Default: 24.96KsMap (
Mapping
): Mapping of Saturated hydraulic conductivity to the inversion model., a SimPEG Mapalpha (
PhysicalProperty
): related to the inverse of the air entry suction [L1], >0, a physical property, Default: 0.036alphaMap (
Mapping
): Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion model., a SimPEG Mapmodel (
Model
): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*, *) or (*)n (
PhysicalProperty
): measure of the poresize distribution, >1, a physical property, Default: 1.56nMap (
Mapping
): Mapping of measure of the poresize 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 [L1], >0 wrt the model.nDeriv (
Derivative
): Derivative of measure of the poresize 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 poresize distribution, >1

property
nMap
¶ Mapping of measure of the poresize distribution, >1 to the inversion model.

property
nDeriv
¶ Derivative of measure of the poresize distribution, >1 wrt the model.

property
alpha
¶ related to the inverse of the air entry suction [L1], >0

property
alphaMap
¶ Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion model.

property
alphaDeriv
¶ Derivative of related to the inverse of the air entry suction [L1], >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 )

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. [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
¶

property