.. 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`. According to
:cite:p:`1995:makhlouf_relation`, an elementary amount of effective rainfall
:math:`dR_E` in the GR3 model will be proportional to an elementary amount of rainfall
:math:`dR'` according to the following equation:

.. math::

    dR_E = \left ( \frac{H_s}{C_g} \right )^2 dR'

Since in this case no evaporation occurs, only rainfall contributes to
changes in the reservoir level. The complement to 1 of :math:`\left (
\frac{H_s}{C_g} \right )^2` will contribute to varying |HS|. Thus, for an
elementary amount of rainfall :math:`dR'`, the elementary variation of water
level :math:`dH_s` will be:

.. math::

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

By integrating this equation over a time step |Deltat|, |HS| will be updated as
follows:

.. math::

    H_s(t  + \Delta t) = \frac{
            H_s(t) + C_g\tanh \left ( \frac{R'}{C_g} \right )
        }{
            1 + \frac{H_s(t)}{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`. The actual evapotranspiration |EA| is
then calculated as follows:

.. math::

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

Since in this case :math:`dR' = 0`, only evaporation contributes to changes in
the reservoir level. Hence, for an elementary change of |PET| :math:`dE'`, the
elementary change of water level :math:`dH_s` is expressed as follows:

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

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

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

For more details about the description of the soil reservoir in the 
GR3 model, see :cite:`1995:makhlouf_relation` and :cite:`1989:edijatno_model`.