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

../../_images/sphx_glr_plot_richards_celia1990_001.png

Richards Problem

class SimPEG.FLOW.Richards.RichardsProblem.RichardsProblem(mesh, **kwargs)[source]

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

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 Solver

Solver (Property): Numerical Solver, Default: new instance of type

solverOpts = {}
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]

The field given the model.

Parameters

m (numpy.ndarray) – model

Return type

numpy.ndarray

Returns

u, the fields

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]

Effect of J(m) on a vector v.

Parameters
Return type

numpy.ndarray

Returns

Jv

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

Effect of transpose of J(m) on a vector v.

Parameters
Return type

numpy.ndarray

Returns

JTv

Richards Survey

class SimPEG.FLOW.Richards.RichardsSurvey.BaseRichardsRx(locs, times)[source]

Bases: SimPEG.Survey.BaseTimeRx

Richards Receiver Object

class SimPEG.FLOW.Richards.RichardsSurvey.PressureRx(locs, times)[source]

Bases: SimPEG.FLOW.Richards.RichardsSurvey.BaseRichardsRx

Richards Receiver Object

deriv(U, prob, du_dm_v=None, v=None, adjoint=False)[source]
class SimPEG.FLOW.Richards.RichardsSurvey.SaturationRx(locs, times)[source]

Bases: SimPEG.FLOW.Richards.RichardsSurvey.BaseRichardsRx

Richards Receiver Object

deriv(U, prob, du_dm_v=None, v=None, adjoint=False)[source]
class SimPEG.FLOW.Richards.RichardsSurvey.RichardsSurvey(rxList, **kwargs)[source]

Bases: SimPEG.Survey.BaseSurvey

rxList = None
property nD

Number of data

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.

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)
deriv(U, du_dm_v=None, v=None)[source]

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)[source]

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)

Vadose Zone Empirical Relationships

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

plot(ax=None)[source]
class SimPEG.FLOW.Richards.Empirical.BaseHydraulicConductivity(mesh, **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 (*, *)

plot(ax=None)[source]
class SimPEG.FLOW.Richards.Empirical.Haverkamp_theta(mesh, **kwargs)[source]

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)[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(mesh, **kwargs)[source]

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)[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(mesh, **kwargs)[source]

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)[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(mesh, **kwargs)[source]

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