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)

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 <type ‘float’>, <type ‘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 <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 Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*) or (*, *)
Dz
Jfull(*args, **kwargs)
Jtvec(*args, **kwargs)
Jvec(*args, **kwargs)
Solver

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

boundary_conditions

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

debug

debug (Boolean) – Show all messages, a boolean, Default – False

diagsJacobian(*args, **kwargs)

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(*args, **kwargs)
getBoundaryConditions(ii, u_ii)
getResidual(*args, **kwargs)

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

method

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

root_finder

Root-finding Algorithm

root_finder_max_iter

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

root_finder_tol

root_finder_tol (Float) – tolerance of the root_finder, a float, Default – 0.0001

solverOpts = {}
water_retention

water_retention (BaseWaterRetention) – water retention curve, an instance of BaseWaterRetention

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

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

derivAdjoint(*args, **kwargs)

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

dpred(*args, **kwargs)

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.
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 <type ‘float’>, <type ‘int’> with shape (*) or (*, *)
plot(ax=None)
class SimPEG.FLOW.Richards.Empirical.BaseWaterRetention(mesh, **kwargs)

Bases: SimPEG.FLOW.Richards.Empirical.NonLinearModel

Optional Properties:

  • model (Model): Inversion model., a numpy array of <type ‘float’>, <type ‘int’> with shape (*) or (*, *)
plot(ax=None)
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 mass-conservative numerical solution for the unsaturated flow equation.” Water Resources Research 26.7 (1990): 1483-1496.

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 <type ‘float’>, <type ‘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.

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 <type ‘float’>, <type ‘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.
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 [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.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 <type ‘float’>, <type ‘int’> with shape (*) or (*, *)
counter = None
mesh = None
nP

Number of parameters in the 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 [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 (*) 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.
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.

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

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

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

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