Gr rainfall intensity

doc/source/_static/bib/references.bib

@@ -1295,5 +1295,22 @@ @book{monnier2021coursevariational
   url = {https://www.math.univ-toulouse.fr/~jmonnie/Enseignement/CourseVDA.pdf}
 }
 
+@article{Astagneau_2022,
+author = {Paul C. Astagneau and François Bourgin and Vazken Andréassian and Charles Perrin},
+title = {Catchment response to intense rainfall: Evaluating modelling hypotheses},
+journal = {Hydrological Processes},
+volume = {26},
+number = {e14676},
+year  = {2022},
+doi = {10.1080/02626667.2021.1923720},
+
+URL = { 
+        https://doi.org/10.1002/hyp.14676
+
+},
+eprint = { 
+        https://doi.org/10.1002/hyp.14676
+}
+}
 
 @Comment{jabref-meta: databaseType:bibtex;}

doc/source/math_num_documentation/forward_structure.rst

@@ -351,8 +351,120 @@ Hydrological processes can be described at pixel scale in `smash` with one of th
     Same as ``gr4`` transfer, see :ref:`GR4 Transfer <math_num_documentation.forward_structure.hydrological_module.gr4>`
 
 
-.. _math_num_documentation.forward_structure.hydrological_module.gr6:
+.. _math_num_documentation.forward_structure.hydrological_module.gr5_ri:
+
+.. dropdown:: gr5_ri (Génie Rural 5 with rainfall intensity terms)
+    :animate: fade-in-slide-down
+
+    This hydrological module is derived from the model introduced in :cite:p:`Astagneau_2022`.
+
+
+    **Production**
+
+
+    In the classical gr production reservoir formulation, the instantaneous production rate is the ratio between the state and the capacity of the reservoir :
+    :math:`\eta = \left( \frac{h_p}{c_p} \right)^2`. The integration of the complementary of :math:`\eta` the rainfall infiltration :math:`p_s`.
+
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+            &p_s = \int_{t-\Delta t}^{t} (1 - \eta) dt \\
+
+        \end{eqnarray}
+
+    We assume that the production reservoir receives a rainfall of :math:`p_n` at time :math:`t`, then
+
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+            &p_s = & c_p \tanh\left(\frac{p_n}{c_p}\right) \frac{1 - \left( \frac{h_p}{c_p} \right)^2}{1 + \frac{h_p}{c_p} \tanh\left( \frac{p_n}{c_p} \right)} \\
+
+        \end{eqnarray}
+
+    To improve runoff production by a gr reservoir, 
+    even with low production level in dry condition, 
+    in the case of high rainfall intensity, :cite:p:`Astagneau_2022` suggests a modification 
+    of the infiltration rate :math:p_s depending on rainfall intensity :math:p_n:
+    :math:`\eta = \left( 1 - \gamma \right) \left( \frac{h_p}{c_p} \right)^2 + \gamma` with :math:`\gamma = 1 - \exp(-p_n \times \alpha_1)`
+    and :math:`\alpha_1` in :math:`mm` per time unit.
+
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+        &p_s& &=& &\int_{t-\Delta t}^{t} (1 - \eta) dt\\
+
+        && &=& &\int_{t-\Delta t}^{t} \left(1 - (1-\gamma) \left(\frac{h_p}{c_p} \right)^2 \right) dt - \int_{t-\Delta t}^{t} \gamma dt\\
+
+        && &=& &\left[ \frac{ c_p }{ \sqrt{1-\gamma} } \tanh \left( \frac{\sqrt{1-\gamma} \  h_p}{c_p} \right) \right]_{t-\Delta t}^t - \gamma \Delta t
+
+        \end{eqnarray}
+
+
+    We denote :math:`\lambda := \sqrt{1 - \gamma}`, then
 
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+        \tanh \left( \lambda \frac{h_p + p_n}{c_p} \right) - \tanh\left( \lambda \frac{h_p}{c_p} \right) &=& 
+        \tanh \left( \lambda \frac{p_n}{c_p} \right) \left(1 - \tanh \left( \lambda \frac{h_p + p_n}{c_p} \right) \tanh \left( \lambda \frac{h_p}{c_p} \right) \right) \\
+        &=& \tanh \left( \lambda \frac{p_n}{c_p} \right) \left(1 - \frac{ \tanh \left( \lambda \frac{h_p}{c_p} \right) + \tanh \left( \lambda \frac{p_n}{c_p} \right) } { 1 + \tanh \left( \lambda \frac{h_p}{c_p} \right) \tanh \left( \lambda \frac{p_n}{c_p} \right) } \tanh \left( \lambda \frac{h_p}{c_p} \right) \right) \\
+        &\sim& \tanh \left( \lambda \frac{p_n}{c_p} \right) \left(1 - \frac{ \lambda \frac{h_p}{c_p} + \tanh \left( \lambda \frac{p_n}{c_p} \right) } { 1 + \lambda \frac{h_p}{c_p} \tanh \left( \lambda \frac{p_n}{c_p} \right) }  \lambda \frac{h_p}{c_p} \right) \\
+        &=& \tanh \left( \lambda \frac{p_n}{c_p} \right) \frac{1 - \left( \lambda \frac{h_p}{c_p} \right)^2}{1 + \lambda \frac{h_p}{c_p} \tanh \left( \lambda \frac{p_n}{c_p} \right)}
+        \end{eqnarray}
+
+    Thus
+
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+        p_s &=& \frac{c_p}{\lambda} \tanh \left( \lambda \frac{p_n}{c_p} \right) \frac{1 - \left( \lambda \frac{h_p}{c_p} \right)^2}{1 + \lambda \frac{h_p}{c_p} \tanh \left( \lambda \frac{p_n}{c_p} \right)} - \gamma \Delta t
+        \end{eqnarray}
+
+
+    .. note::
+
+        Note that if :math:`\alpha_1 = 0`, we return to the general writting of the instantaneous production rate.
+
+
+    **Transfer**
+
+    In context of high rainfall intensities triggering flash flood responses, it is crucial to account for fast dynamics related to surface/hypodermic runoff 
+    and slower responses due to delayed/deeper flows (e.g. https://doi.org/10.5194/hess-26-5185-2022). 
+    Following :cite:p:`Astagneau_2022` for a lumped GR model, we introduce at pixel scale in smash a function to modify the partitionning between fast 
+    and slower transfert branches depending on rainfall intensity of the current time step only (small pixel size):
+
+    .. math::
+        :nowrap:
+
+        \begin{eqnarray}
+
+            &p_{rr}& =& (1 - Q_9)(p_r + p_{erc}) + l_{exc}\\
+            &p_{rd}& =& Q_9(p_r + p_{erc}) \\
+            &Q_9& =& 0.9 \tanh(\alpha_2 p_n)^2 + 0.1
+
+        \end{eqnarray}
+
+    with :math:`\alpha_2` in :math:`mm` per time unit.
+
+
+    .. note::
+
+        If :math:`\alpha_2 = 0`, we return to the ``gr-4/gr-5`` writting of the transfer.
+        If :math:`\alpha_2 = \alpha_1 = 0`, it is equivalent to ``gr-5`` structure.
+
+
+
+.. _math_num_documentation.forward_structure.hydrological_module.gr6:
 
 .. dropdown:: gr6 (Génie Rural 6)
     :animate: fade-in-slide-down

smash/_constant.py

@@ -64,10 +64,12 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
 
 HYDROLOGICAL_MODULE = [
     "gr4",
+    "gr4_ri",
     "gr4_mlp",
     "gr4_ode",
     "gr4_ode_mlp",
     "gr5",
+    "gr5_ri",
     "gr6",
     "grc",
     "grd",
@@ -97,10 +99,12 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
         HYDROLOGICAL_MODULE,
         [
             ["ci", "cp", "ct", "kexc"],  # % gr4
+            ["ci", "cp", "ct", "alpha1", "alpha2", "kexc"],  # % gr4_ri
             ["ci", "cp", "ct", "kexc"],  # % gr4_mlp
             ["ci", "cp", "ct", "kexc"],  # % gr4_ode
             ["ci", "cp", "ct", "kexc"],  # % gr4_ode_mlp
             ["ci", "cp", "ct", "kexc", "aexc"],  # % gr5
+            ["ci", "cp", "ct", "alpha1", "alpha2", "kexc", "aexc"],  # % gr5_ri
             ["ci", "cp", "ct", "be", "kexc", "aexc"],  # % gr6
             ["ci", "cp", "ct", "cl", "kexc"],  # % grc
             ["cp", "ct"],  # % grd
@@ -139,10 +143,12 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
         HYDROLOGICAL_MODULE,
         [
             ["hi", "hp", "ht"],  # % gr4
+            ["hi", "hp", "ht"],  # % gr4_ri
             ["hi", "hp", "ht"],  # % gr4_mlp
             ["hi", "hp", "ht"],  # % gr4_ode
             ["hi", "hp", "ht"],  # % gr4_ode_mlp
             ["hi", "hp", "ht"],  # % gr5
+            ["hi", "hp", "ht"],  # % gr5_ri
             ["hi", "hp", "ht", "he"],  # % gr6
             ["hi", "hp", "ht", "hl"],  # % grc
             ["hp", "ht"],  # % grd
@@ -183,10 +189,12 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
         HYDROLOGICAL_MODULE,
         [
             ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "qr", "qd", "qt"],  # % gr4
+            ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "qr", "qd", "qt"],  # % gr4-ri
             ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "qr", "qd", "qt"],  # % gr4_mlp
             ["pn", "en", "lexc", "qt"],  # % gr4_ode
             ["pn", "en", "lexc", "qt"],  # % gr4_ode_mlp
             ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "qr", "qd", "qt"],  # % gr5
+            ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "qr", "qd", "qt"],  # % gr5-ri
             ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "pre", "qr", "qd", "qe", "qt"],  # % gr6
             ["pn", "en", "pr", "perc", "lexc", "prr", "prd", "prl", "qr", "qd", "ql", "qt"],  # % grc
             ["ei", "pn", "en", "pr", "perc", "prr", "qr", "qt"],  # % grd
@@ -250,13 +258,15 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
 
 RR_PARAMETERS = [
     "kmlt",  # % ssn
-    "ci",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc)
-    "cp",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc, grd)
-    "ct",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc, grd)
+    "ci",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6)
+    "cp",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6, grd)
+    "ct",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6, grd)
+    "alpha1",  # % (gr4_ri, gr5_ri)
+    "alpha2",  # % (gr4_ri, gr5_ri)
     "cl",  # % grc
-    "be",  # % gr6
-    "kexc",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc)
-    "aexc",  # % (gr5, gr6)
+    "be",  # % (gr6)
+    "kexc",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6)
+    "aexc",  # % (gr5, gr5_ri, gr6)
     "ca",  # % loieau
     "cc",  # % loieau
     "kb",  # % loieau
@@ -276,9 +286,9 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
 
 RR_STATES = [
     "hs",  # % ssn
-    "hi",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc)
-    "hp",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc, grd)
-    "ht",  # % (gr4, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr6, grc, grd)
+    "hi",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6)
+    "hp",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6, grd)
+    "ht",  # % (gr4, gr4_ri, gr4_mlp, gr4_ode, gr4_ode_mlp, gr5, gr5_ri, grc, gr6, grd)
     "hl",  # % grc
     "he",  # % gr6
     "ha",  # % loieau
@@ -302,6 +312,8 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
             (0, np.inf),  # % ci
             (0, np.inf),  # % cp
             (0, np.inf),  # % ct
+            (0, np.inf),  # % alpha1
+            (0, np.inf),  # % alpha2
             (0, np.inf),  # % cl
             (0, np.inf),  # % be
             (-np.inf, np.inf),  # % kexc
@@ -361,6 +373,8 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
             1e-6,  # % ci
             200,  # % cp
             500,  # % ct
+            3e-4,  # % alpha1
+            1e-3,  # % alpha2
             500,  # % cl
             10,  # % be
             0,  # % kexc
@@ -418,6 +432,8 @@ def get_neurons_from_hydrological_module(hydrological_module: str, hidden_neuron
             (1e-6, 1e2),  # % ci
             (1e-6, 1e3),  # % cp
             (1e-6, 1e3),  # % ct
+            (1e-6, 5e-4),  # % alpha1
+            (1e-5, 1.0),  # % alpha2
             (1e-6, 1e3),  # % cl
             (1e-3, 20),  # % be
             (-50, 50),  # % kexc