The Soil

Introduction

The soil module can be regarded as the core of VAMPS. The model was developed using the SWAP94 model as a framework. About 80 % of the soil module consists of hand-translated Fortran code from the original SWAP94 programme. The head-based solution of Richards’ equation remains more or less the same in the present version, extended with a band-diagonal solver and an optional step to regain full machine precision [15] when a zero pivot occurs in the tridiagonal solution.

Solving Richards’ Equation

The soil module solves Richards’ equation using a \(\psi\)-based form. A differential moisture capacity \(C(h)\) has been introduced as the single independent variable. This differential soil moisture capacity equals the slope of the soil-moisture retention curve, \(d\theta/d\psi\).

\[\frac{\partial \psi}{\partial t} = \frac{1}{C(\psi)} \frac{\partial}{\partial z} \left[ K(\psi) \left( \frac{\partial \psi}{\partial z} - 1 \right) \right]\]

A sink term \(S\) has been added to accommodate water extraction by roots and lateral flow:

\[\frac{\partial \psi}{\partial t} = \frac{1}{C(\psi)} \frac{\partial}{\partial z} \left[ K(\psi) \left( \frac{\partial \psi}{\partial z} - 1 \right) \right] - \frac{S(\psi)}{C(\psi)}\]

Where:

\(\psi\)	suction head [cm]
\(t\)	time [days]
\(C\)	differential moisture capacity (\(d\theta/d\psi\))
\(S\)	sink term
\(z\)	height [cm]

For a detailed description see [16] and [17].

VAMPS estimates the largest time-step possible from the previous dt and the thetol parameter (maximum allowed change in \(\theta\) per time-step) across all soil layers — the layer with the largest rate of change determines the maximum dt:

/*-
 * void timestep(int step)
 * Calculation of timestep (dt) depending on theta changes and new time (t)
 */
double
timestep (int step)
{
  int i;

  /* Store previous values */
  tm1 = t;
  dtm1 = dt;

  /* start with max dt */
  dt = dtmax;

  mdtheta = fabs (thetm1[0] - theta[0]);
  for (i = 0; i < layers; i++)
    {
      dtheta = fabs (thetm1[i] - theta[i]);
      mdtheta = mdtheta < dtheta ? dtheta : mdtheta;
    }

  if (mdtheta < 1.0E-6)
    dt = dtmax;
  else
      dt = 10.0 * thetol * dtm1 / mdtheta;

  /* restrict to dtmin and dtmax */
  dt = dt > dtmin ? dt : dtmin;
  dt = dt < dtmax ? dt : dtmax;

return dt;
}

If the result is not accurate enough (controlled by thetol) VAMPS performs iterations. If these fail dt is halved and the procedure repeats until the required accuracy is reached or dtmin is hit. The maximum number of iterations is set with maxitr; the minimum allowed dt with dtmin.

Relation between \(\theta\), \(\psi\), and \(K_{unsat}\)

VAMPS was first developed using the equations of [18] (see also [19]). This remains the default method. For extra speed and flexibility VAMPS can also use lookup tables that describe these relations, allowing any user-defined model to be plugged in.

At present VAMPS can generate tables for the Van Genuchten method, and can read tables generated by the topog_soil program [20] or user-made tables. Only the Broadbridge–White (option 1) and Van Genuchten (option 5) table types from topog_soil have been tested with VAMPS.

Van Genuchten

The default and most thoroughly tested method. The dimensionless water content is:

\[\Theta = \frac{\theta - \theta_r}{\theta_s - \theta_r}\]

Where:

\(\theta_s\)	saturated water content
\(\theta_r\)	residual water content
\(\theta\)	actual soil water content
\(\Theta\)	dimensionless soil water content

A class of \(\Theta(\psi)\) functions [18]:

\[\Theta = \left[ \frac{1}{1 + (\alpha \psi)^n} \right]^m\]

Where \(\psi\) is the pressure head and \(\alpha\), \(n\), \(m\) are soil characteristic parameters.

With \(m = 1 - 1/n\) the following closed-form expression for \(\Theta(\psi)\) is obtained:

\[\begin{split}\Theta = \left\{ \begin{array}{ll} \dfrac{1}{[1 + (\alpha \psi)^n]^{1-1/n}} & \psi < 0 \\[8pt] 1 & \psi \geq 0 \end{array} \right.\end{split}\]

And the hydraulic conductivity:

\[K = K_s \, \Theta^{1/2} \left[ 1 - \left( 1 - \Theta^{n/(n-1)} \right)^{1-1/n} \right]^2\]

The \(K_s \Theta^{1/2}\) term can be replaced by \(K_s \Theta^L\) where \(L\) is a fitting parameter describing the \(\Theta\)–\(K_{unsat}\) relation if unsaturated conductivity measurements are available.

Values for \(\alpha\), \(n\), and \(L\) for most Dutch soil types are given in [21]. For tropical soils a nonlinear regression method [22] can be used to determine \(\alpha\) and \(n\) from a soil moisture retention curve.

Clapp/Hornberger

[23] demonstrated that a power curve conveniently describes soil moisture characteristics. [24] described these functions and presented representative values for soil hydraulic parameters:

\[\psi = \psi_s \left( \frac{\theta}{\theta_s} \right)^{-b}\]

\[\frac{K}{K_s} = \left( \frac{\theta}{\theta_s} \right)^{2b + 2}\]

Where:

\(\psi\)	suction head
\(\psi_s\)	suction head at saturation (air entry value)
\(\theta\)	volumetric water content
\(\theta_s\)	water content at saturation (porosity)
\(b\)	soil characteristic parameter
\(K\)	hydraulic conductivity
\(K_s\)	saturated hydraulic conductivity

Representative values for the Clapp/Hornberger parameters (after Clapp 1978):

Soil Texture	\(\bar{b}\)	\(\bar{\psi_s}\) [cm]	\(\bar{\theta_s}\)	\(\bar{K_s}\) [cm/day]
Sand	4.05	12.2	0.395	1520.64
Loamy sand	4.38	9.0	0.410	1350.72
Sandy loam	4.90	21.8	0.435	299.52
Silt loam	5.3	78.6	0.485	62.21
Loam	5.39	47.8	0.451	60.05
Sandy clay loam	7.12	29.9	0.420	54.43
Silty clay loam	7.75	35.6	0.477	14.69
Clay loam	8.52	63.0	0.476	21.17
Sandy clay	10.4	15.3	0.426	18.72
Silty clay	10.4	49.0	0.492	8.93
Clay	11.4	40.5	0.482	11.09

An input file containing all these soil types is soillib.inp in the src/sllib directory.

Brooks and Corey

[25] proposed a power-law description of the soil moisture characteristic with an explicit air-entry (bubbling) pressure head \(h_b\) below which the soil begins to desaturate. The effective saturation \(S_e\) is defined as:

\[S_e = \frac{\theta - \theta_r}{\theta_s - \theta_r}\]

Where \(\theta_s\) is the saturated water content, \(\theta_r\) is the residual water content, and \(\theta\) is the actual soil water content.

For \(h < h_b\) (unsaturated) \(S_e = (h_b/h)^\lambda\), otherwise \(S_e = 1\).

The unsaturated hydraulic conductivity as a function of \(S_e\) is:

\[K = K_s \, S_e^{\,2/\lambda + 3}\]

Where:

\(h\)	pressure head (negative for unsaturated conditions)
\(h_b\)	air-entry (bubbling) pressure head (\(h_b < 0\))
\(\lambda\)	pore-size distribution index
\(K_s\)	saturated hydraulic conductivity

The conductivity can equivalently be written as a function of pressure head:

\[K(h) = K_s \left(\frac{h_b}{h}\right)^{n}, \quad n = 2 + 3\lambda\]

Converting from Van Genuchten parameters:

Starting approximations are \(\lambda \approx n_\mathrm{vG} - 1\) and \(h_b \approx -1/\alpha\). The two retention curves agree closely for \(S_e \leq 0.80\); the main difference is near saturation where Brooks–Corey has a discrete air-entry while Van Genuchten desaturates smoothly.

See examples/fiji/fiji_bandC.inp for a complete example using Brooks–Corey parameters equivalent to the Van Genuchten Fiji example.

Bottom Boundary Conditions

VAMPS can handle several boundary conditions at the bottom of the profile. Not all conditions have been fully tested in the present version.

Code	Description
0	Daily groundwater table depth (cm) as input
1	Fixed flux at the bottom of the profile
2	Seepage or infiltration from/to deep groundwater
3	Flux calculated as a function of head
4	Fixed head at the bottom of the profile (interpolation between daily values)
5	No-flow bottom boundary
6	Free drainage at the bottom of the profile

The choice of bottom boundary condition is of great importance to the modelling results. In most cases the choice is determined by the available data.

See the bottom variable in the [soil] section of the configuration reference (Configuration File Reference — vamps(5)) for details.