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

(Source code, png, hires.png, pdf)

../../_images/index-11.png

Richards Problem

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

Bases: SimPEG.Problem.BaseTimeProblem

Required

Parameters:
  • Solver (Property) – Numerical Solver, Default: new instance of type
  • boundary_conditions (Array) – boundary conditions, a list or numpy array of <type ‘float’>, <type ‘int’> with shape (*)
  • debug (Bool) – Show all messages, a boolean, Default: False
  • do_newton (Bool) – 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 <type ‘float’>, <type ‘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

Parameters:model (Model) – Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
hydraulic_conductivity

hydraulic conductivity function

water_retention

water retention curve

boundary_conditions

boundary conditions

initial_conditions

initial conditions

debug

Show all messages

Solver

Numerical Solver

solverOpts = {}
method

Formulation used, See notes in Celia et al., 1990

do_newton

Do a Newton iteration vs. a Picard iteration

root_finder_max_iter

Maximum iterations for root_finder iteration

root_finder_tol

tolerance of the root_finder

getBoundaryConditions(ii, u_ii)[source]
root_finder

Root-finding Algorithm

fields(*args, **kwargs)[source]
Dz
diagsJacobian(*args, **kwargs)[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(*args, **kwargs)[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(*args, **kwargs)[source]
Jvec(*args, **kwargs)[source]
Jtvec(*args, **kwargs)[source]

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
nD
dpred(*args, **kwargs)[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.

To use dpred method, SimPEG requires that the prob be specified.
deriv(*args, **kwargs)[source]

The Derivative with respect to the model. To use deriv method, SimPEG requires that the prob be specified.

derivAdjoint(*args, **kwargs)[source]

The adjoint derivative with respect to the model. To use derivAdjoint method, SimPEG requires that the prob be specified.

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

Parameters:model (Model) – Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
counter = None

A SimPEG.Utils.Counter object

mesh = None

A SimPEG Mesh

nP

Number of parameters in the model.

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

Bases: SimPEG.FLOW.Richards.Empirical.NonLinearModel

Optional

Parameters:model (Model) – Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
plot(ax=None)[source]
class SimPEG.FLOW.Richards.Empirical.BaseHydraulicConductivity(mesh, **kwargs)[source]

Bases: SimPEG.FLOW.Richards.Empirical.NonLinearModel

Optional

Parameters:model (Model) – Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*)
plot(ax=None)[source]
class SimPEG.FLOW.Richards.Empirical.Haverkamp_theta(mesh, **kwargs)[source]

Bases: SimPEG.FLOW.Richards.Empirical.BaseWaterRetention

Optional

Parameters:
  • 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 <type ‘float’>, <type ‘int’> with shape (*)
  • 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

Immutable

Attribute alphaDeriv:
 (Derivative) - Derivative of wrt the model.
Attribute betaDeriv:
 (Derivative) - Derivative of wrt the model.
Attribute theta_rDeriv:
 (Derivative) - Derivative of residual water content [L3L-3] wrt the model.
Attribute theta_sDeriv:
 (Derivative) - Derivative of saturated water content [L3L-3] 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]
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.

theta_r

residual water content [L3L-3]

theta_rDeriv

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

theta_rMap

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

theta_s

saturated water content [L3L-3]

theta_sDeriv

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

theta_sMap

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

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

Bases: SimPEG.FLOW.Richards.Empirical.BaseHydraulicConductivity

Optional

Parameters:
  • 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 <type ‘float’>, <type ‘int’> with shape (*)

Immutable

Attribute ADeriv:
 (Derivative) - Derivative of fitting parameter wrt the model.
Attribute KsDeriv:
 (Derivative) - Derivative of Saturated hydraulic conductivity wrt the model.
Attribute gammaDeriv:
 (Derivative) - Derivative of fitting parameter wrt the model.
derivU(u)[source]
derivM(u)[source]
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.

gamma

fitting parameter

gammaDeriv

Derivative of fitting parameter wrt the model.

gammaMap

Mapping of fitting parameter to the inversion model.

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.

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

Parameters:
  • 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 <type ‘float’>, <type ‘int’> with shape (*)
  • 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

Immutable

Attribute alphaDeriv:
 (Derivative) - Derivative of related to the inverse of the air entry suction [L-1], >0 wrt the model.
Attribute nDeriv:
 (Derivative) - Derivative of measure of the pore-size distribution, >1 wrt the model.
Attribute theta_rDeriv:
 (Derivative) - Derivative of residual water content [L3L-3] wrt the model.
Attribute theta_sDeriv:
 (Derivative) - Derivative of saturated water content [L3L-3] 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]
alpha

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

alphaDeriv

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

alphaMap

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

n

measure of the pore-size distribution, >1

nDeriv

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

nMap

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

theta_r

residual water content [L3L-3]

theta_rDeriv

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

theta_rMap

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

theta_s

saturated water content [L3L-3]

theta_sDeriv

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

theta_sMap

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

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

Bases: SimPEG.FLOW.Richards.Empirical.BaseHydraulicConductivity

Optional

Parameters:
  • 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 <type ‘float’>, <type ‘int’> with shape (*)
  • 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

Immutable

Attribute IDeriv:
 (Derivative) - Derivative of wrt the model.
Attribute KsDeriv:
 (Derivative) - Derivative of Saturated hydraulic conductivity wrt the model.
Attribute alphaDeriv:
 (Derivative) - Derivative of related to the inverse of the air entry suction [L-1], >0 wrt the model.
Attribute nDeriv:
 (Derivative) - Derivative of measure of the pore-size distribution, >1 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]
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 [L-1], >0

alphaDeriv

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

alphaMap

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

n

measure of the pore-size distribution, >1

nDeriv

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

nMap

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

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]

sand
loamy_sand
sandy_loam
loam
silt_loam
sandy_clay_loam
clay_loam
silty_clay_loam
sandy_clay
silty_clay
clay