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 :cite:`press1992B` when a zero pivot occurs in the
tridiagonal solution.


Solving Richards' Equation
--------------------------

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

.. math::

   \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 :math:`S` has been added to accommodate water extraction by roots
and lateral flow:

.. math::

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

.. list-table::
   :widths: 20 80

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

For a detailed description see :cite:`feddes1978273` and :cite:`belmans1983272`.

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

.. code-block:: c

   /*-
    * 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 :math:`\theta`, :math:`\psi`, and :math:`K_{unsat}`
----------------------------------------------------------------------

VAMPS was first developed using the equations of :cite:`genuchten1980179` (see also
:cite:`mualem1976261`).  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 :cite:`beverly1994274` 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:

.. math::

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

Where:

.. list-table::
   :widths: 20 80

   * - :math:`\theta_s`
     - saturated water content
   * - :math:`\theta_r`
     - residual water content
   * - :math:`\theta`
     - actual soil water content
   * - :math:`\Theta`
     - dimensionless soil water content

A class of :math:`\Theta(\psi)` functions :cite:`genuchten1980179`:

.. math::

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

Where :math:`\psi` is the pressure head and :math:`\alpha`, :math:`n`,
:math:`m` are soil characteristic parameters.

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

.. math::

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

And the hydraulic conductivity:

.. math::

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

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

Values for :math:`\alpha`, :math:`n`, and :math:`L` for most Dutch soil
types are given in :cite:`wosten1987287`.  For tropical soils a nonlinear
regression method :cite:`marquardt1963N` can be used to determine :math:`\alpha`
and :math:`n` from a soil moisture retention curve.


Clapp/Hornberger
~~~~~~~~~~~~~~~~

:cite:`campbell1974` demonstrated that a power curve conveniently describes soil
moisture characteristics.  :cite:`clapp1978263` described these functions and
presented representative values for soil hydraulic parameters:

.. math::

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

.. math::

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

Where:

.. list-table::
   :widths: 20 80

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

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

.. list-table::
   :header-rows: 1
   :widths: 30 15 15 15 15

   * - Soil Texture
     - :math:`\bar{b}`
     - :math:`\bar{\psi_s}` [cm]
     - :math:`\bar{\theta_s}`
     - :math:`\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
~~~~~~~~~~~~~~~~

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

.. math::

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

Where :math:`\theta_s` is the saturated water content, :math:`\theta_r` is
the residual water content, and :math:`\theta` is the actual soil water
content.

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

The unsaturated hydraulic conductivity as a function of :math:`S_e` is:

.. math::

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

Where:

.. list-table::
   :widths: 20 80

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

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

.. math::

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

**Converting from Van Genuchten parameters:**

Starting approximations are :math:`\lambda \approx n_\mathrm{vG} - 1` and
:math:`h_b \approx -1/\alpha`.  The two retention curves agree closely for
:math:`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.

.. list-table::
   :header-rows: 1
   :widths: 15 85

   * - 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 (:doc:`../reference/config_reference`) for details.