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): Richards Problem¶

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

Required Properties:

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

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

boundary_conditions

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

debug

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

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

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)

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

deriv(U, prob, du_dm_v=None, v=None, adjoint=False)
class SimPEG.FLOW.Richards.RichardsSurvey.RichardsSurvey(rxList, **kwargs)
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)

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

class SimPEG.FLOW.Richards.Empirical.BaseHydraulicConductivity(mesh, **kwargs)

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘int’> with shape (*) or (*, *)
plot(ax=None)
class SimPEG.FLOW.Richards.Empirical.BaseWaterRetention(mesh, **kwargs)

Optional Properties:

• model (Model): Inversion model., a numpy array of <class ‘float’>, <class ‘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)

Optional Properties:

Other Properties:

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)

Optional Properties:

Other Properties:

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

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)

Optional Properties:

Other Properties:

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)

Optional Properties:

Other Properties:

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)