.. include:: ../notations/notations.rst

.. _soil-reservoir:

==============
Soil reservoir
==============

The soil reservoir takes as input rainfall |R| (mm) and |PET| |EP| (mm) to
calculate the effective rainfall |RE| (mm), the |AET| |EA| (mm) and the |UPET|
|EU| (mm) of the watershed. |HS| (mm) is the water level reservoir (soil
moisture storage).

In |rameau|, there are actually two methods to compute |HS|: the
Thornthwaite method :cite:p:`1948:thornthwaite_approach` and the GR3
method :cite:p:`1989:edijatno_model`.

.. _soil_reservoir_parameters:

Physical parameters
===================

Soil capacity of the Thornthwaite reservoir
-------------------------------------------

The soil water holding capacity |CT| (mm) is the maximum water level allowed
in the soil reservoir.

Should be optimised. Default value is 70 mm.

.. code-block:: toml

    [watershed.1]
    thornthwaite.capacity = { value = 70.0, lower = 0, upper = 500, opti = true, sameas = 0 } 

Soil capacity of the progressive reservoir
------------------------------------------

The soil water holding capacity |CG| (mm) is the maximum water level allowed
in the soil reservoir.

Should be optimised. Default value is 70 mm.

.. code-block:: toml

    [watershed.1]
    progressive.capacity = { value = 70.0, lower = 0, upper = 500, opti = true, sameas = 0 } 

Exponential decrease of |PET|
-----------------------------

If true, the GR3 formulation to compute |AET| is triggered only if the
reservoir saturation is lower than 50% of its capacity.

Default value is false.

.. code-block:: toml

    [watershed.1]
    progressive.pet_decrease = false

Calculating effective rainfall
==============================

Thornthwaite model
------------------

If |R| > |EP|, then

.. math::

    E_A & = E_P \\
    E_U & = 0 \\
    H_s & = H_s + R - E_P

Then, if |HS| > |CT| then

.. math::

    R_E & = H_s - C_t \\
    H_s & = C_t

else

.. math::
    
    R_E = 0

Conversely, if |EP| > |R| then

.. math::

    E_A & = \min(R + H_s, E_P) \\
    H_s & = \max(0, H_s - R + E_P) \\
    E_U & = \max(0, E_P - E_A) \\
    R_E & = 0.

GR3 model
---------

If |R| > |EP|, the amount of water that feeds into the soil reservoir over a
time step is equal to :math:`R'= R - E_P`. Then, |RE| (mm) is computed as
follows:

.. math::

    R_E = \left ( \frac{H_s}{C_g} \right )^2 R'

The complement to 1 of :math:`(1 - \left ( \frac{H_s}{C_g})^2 \right )R'`
will contribute to varying |HS|. For an elementary amount of rainfall
:math:`dR'` :cite:p:`1995:makhlouf_relation`:

.. math::

    dS = \left ( 1 - \left ( \frac{H_s}{C_g} \right )^2 \right )dR'.

By integrating this equation over the time step:

.. math::

    H_s = \frac{
            H_s + C_g\tanh \left ( \frac{R'}{C_g} \right )
        }{
            1 + \frac{H_s}{C_g}\tanh \left (\frac{R'}{C_g} \right )
        }.

Then,

.. math::

    E_A & = E_P \\
    E_U & = 0.
    
Conversely, if |R| < |EP|, the soil reservoir is subject to the
|PET| :math:`E'= R - E_P` and 

.. math::

    E_A = \frac{H_s}{C_g} \left (2 - \frac{H_s}{C_g} \right )E'.

For an elementary change of |PET| :math:`dE`, the reduction of |HS|
is expressed as follows :cite:p:`1995:makhlouf_relation`:

.. math::

    dH_s = -\frac{H_s}{C_g} \left (2 - \frac{H_s}{C_g} \right )dE.

By integrating this equation over the time step |Deltat| (s)

.. math::
    :label: eq_s2

    H_s(t + \Delta t) = H_s(t)\frac{
        1 - \tanh \left ( \frac{E'}{C_g} \right )
        } {
            1 + \left (1 - \frac{H_s(t)}{C_g} \right )
                \tanh \left (\frac{E'}{C_g} \right )
        }

and then, to make sure that |HS| is always positive:

.. math::

    H_s(t + \Delta t) & = max(0, H_s(t + \Delta t)) \\
    E_A & = R - (H_s(t + \Delta t)- H_s(t)) \\
    E_U & = E_P - E_A.

In |rameau|, it is possible to trigger the |AET| decrease depending on 
the reservoir saturation only when it is below 50% of its capacity.
In this case, we have:

.. math:: 

    C^{50}_g = 0.5C_g

If :math:`H_s(t) > C^{50}_g`, |PET| is fully satisfied similarly to the
Thornthwaite approach:

.. math::

    H_s(t + \Delta t)= H_s(t) - E'

Else, :math:`H_s(t + \Delta t)` is computed using Equation :eq:`eq_s2`
replacing |CG| with :math:`C^{50}_g`.