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 Problem¶

class
SimPEG.FLOW.Richards.RichardsProblem.
RichardsProblem
(mesh, **kwargs)¶ Bases:
SimPEG.Problem.BaseTimeProblem
Required Properties:
 Solver (
Property
): Numerical Solver, Default: new instance of type  boundary_conditions (
Array
): boundary conditions, a list or numpy array of <class ‘float’>, <class ‘int’> with shape (*)  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 (*)  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  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 (*, *)

Dz
¶

Jfull
(m=None, f=None)¶

Jtvec
(m, v, f=None)¶

Jvec
(m, v, f=None)¶

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

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  ' ' ' ' ' '

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

fields
(m=None)¶

getBoundaryConditions
(ii, u_ii)¶

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

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

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

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

root_finder
¶ Rootfinding Algorithm

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

solverOpts
= {}¶

water_retention
¶ water_retention (
BaseWaterRetention
): water retention curve, an instance of BaseWaterRetention
 Solver (
Richards Survey¶

class
SimPEG.FLOW.Richards.RichardsSurvey.
BaseRichardsRx
(locs, times)¶ Bases:
SimPEG.Survey.BaseTimeRx
Richards Receiver Object

class
SimPEG.FLOW.Richards.RichardsSurvey.
PressureRx
(locs, times)¶ Bases:
SimPEG.FLOW.Richards.RichardsSurvey.BaseRichardsRx
Richards Receiver Object

deriv
(U, prob, du_dm_v=None, v=None, adjoint=False)¶


class
SimPEG.FLOW.Richards.RichardsSurvey.
RichardsSurvey
(rxList, **kwargs)¶ Bases:
SimPEG.Survey.BaseSurvey

deriv
(U, du_dm_v=None, v=None)¶ The Derivative with respect to the model.
Note
To use survey.deriv(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:
survey.pair(myProblem)

derivAdjoint
(U, v=None)¶ The adjoint derivative with respect to the model.
Note
To use survey.derivAdjoint(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:
survey.pair(myProblem)

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.
Note
To use survey.dpred(), SimPEG requires that a problem be bound to the survey. If a problem has not been bound, an Exception will be raised. To bind a problem to the Data object:
survey.pair(myProblem)

nD
¶

rxList
= None¶


class
SimPEG.FLOW.Richards.RichardsSurvey.
SaturationRx
(locs, times)¶ Bases:
SimPEG.FLOW.Richards.RichardsSurvey.BaseRichardsRx
Richards Receiver Object

deriv
(U, prob, du_dm_v=None, v=None, adjoint=False)¶

Vadose Zone Empirical Relationships¶

class
SimPEG.FLOW.Richards.Empirical.
BaseHydraulicConductivity
(mesh, **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)¶
 model (

class
SimPEG.FLOW.Richards.Empirical.
BaseWaterRetention
(mesh, **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)¶
 model (

class
SimPEG.FLOW.Richards.Empirical.
HaverkampParams
¶ Bases:
object
Holds some default parameterizations for the Haverkamp model.

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.


class
SimPEG.FLOW.Richards.Empirical.
Haverkamp_k
(mesh, **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.

A
¶ fitting parameter

ADeriv
¶ Derivative of fitting parameter wrt the model.

AMap
¶ Mapping of fitting parameter to the inversion model.

Ks
¶ Saturated hydraulic conductivity

KsDeriv
¶ Derivative of Saturated hydraulic conductivity wrt the model.

KsMap
¶ Mapping of Saturated hydraulic conductivity to the inversion model.

derivM
(u)¶

derivU
(u)¶

gamma
¶ fitting parameter

gammaDeriv
¶ Derivative of fitting parameter wrt the model.

gammaMap
¶ Mapping of fitting parameter to the inversion model.
 A (

class
SimPEG.FLOW.Richards.Empirical.
Haverkamp_theta
(mesh, **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 [L3L3], a physical property, Default: 0.075  theta_rMap (
Mapping
): Mapping of residual water content [L3L3] to the inversion model., a SimPEG Map  theta_s (
PhysicalProperty
): saturated water content [L3L3], a physical property, Default: 0.287  theta_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.

alpha
¶

alphaDeriv
¶ Derivative of wrt the model.

alphaMap
¶ Mapping of to the inversion model.

beta
¶

betaDeriv
¶ Derivative of wrt the model.

betaMap
¶ Mapping of to the inversion 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)¶

theta_r
¶ residual water content [L3L3]

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

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

theta_s
¶ saturated water content [L3L3]

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

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

class
SimPEG.FLOW.Richards.Empirical.
NonLinearModel
(mesh, **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¶

mesh
= None¶

nP
¶ Number of parameters in the model.
 model (

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]

clay
¶

clay_loam
¶

loam
¶

loamy_sand
¶

sand
¶

sandy_clay
¶

sandy_clay_loam
¶

sandy_loam
¶

silt_loam
¶

silty_clay
¶

silty_clay_loam
¶


class
SimPEG.FLOW.Richards.Empirical.
Vangenuchten_k
(mesh, **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 [L1], >0, a physical property, Default: 0.036  alphaMap (
Mapping
): Mapping of related to the inverse of the air entry suction [L1], >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 poresize distribution, >1, a physical property, Default: 1.56  nMap (
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.

I
¶

IDeriv
¶ Derivative of wrt the model.

IMap
¶ Mapping of to the inversion model.

Ks
¶ Saturated hydraulic conductivity

KsDeriv
¶ Derivative of Saturated hydraulic conductivity wrt the model.

KsMap
¶ Mapping of Saturated hydraulic conductivity to the inversion model.

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

alphaDeriv
¶ Derivative of related to the inverse of the air entry suction [L1], >0 wrt the model.

alphaMap
¶ Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion 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)¶

n
¶ measure of the poresize distribution, >1

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

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

class
SimPEG.FLOW.Richards.Empirical.
Vangenuchten_theta
(mesh, **kwargs)¶ 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.036  alphaMap (
Mapping
): Mapping of related to the inverse of the air entry suction [L1], >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 poresize distribution, >1, a physical property, Default: 1.56  nMap (
Mapping
): Mapping of measure of the poresize distribution, >1 to the inversion model., a SimPEG Map  theta_r (
PhysicalProperty
): residual water content [L3L3], a physical property, Default: 0.078  theta_rMap (
Mapping
): Mapping of residual water content [L3L3] to the inversion model., a SimPEG Map  theta_s (
PhysicalProperty
): saturated water content [L3L3], a physical property, Default: 0.43  theta_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.

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

alphaDeriv
¶ Derivative of related to the inverse of the air entry suction [L1], >0 wrt the model.

alphaMap
¶ Mapping of related to the inverse of the air entry suction [L1], >0 to the inversion 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)¶

n
¶ measure of the poresize distribution, >1

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

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

theta_r
¶ residual water content [L3L3]

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

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

theta_s
¶ saturated water content [L3L3]

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

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

SimPEG.FLOW.Richards.Empirical.
haverkamp
(mesh, **kwargs)¶

SimPEG.FLOW.Richards.Empirical.
van_genuchten
(mesh, **kwargs)¶