diff --git a/doc/source/user_guide/others/input_data_convention.rst b/doc/source/user_guide/others/input_data_convention.rst
index 7f0b2e22..c648a41c 100644
--- a/doc/source/user_guide/others/input_data_convention.rst
+++ b/doc/source/user_guide/others/input_data_convention.rst
@@ -130,7 +130,7 @@ Below is the most basic directory structure you can have, with one subdirectory
     │   ├── V3504010.csv
     │   └── ...
     └── descriptor
-        ├── slope.ti
+        ├── slope.tif
         └── dd.tif
 
 This results in the following ``setup``:
diff --git a/doc/source/user_guide/quickstart/cance_first_simulation.rst b/doc/source/user_guide/quickstart/cance_first_simulation.rst
index 160ba3db..2a0fab59 100644
--- a/doc/source/user_guide/quickstart/cance_first_simulation.rst
+++ b/doc/source/user_guide/quickstart/cance_first_simulation.rst
@@ -670,17 +670,19 @@ First, several information were displayed on the screen during optimization
 
 .. code-block:: text
 
-    At iterate      0    nfg =     1    J =      0.643190    ddx = 0.64
-    At iterate      1    nfg =    30    J =      0.097397    ddx = 0.64
-    At iterate      2    nfg =    59    J =      0.052158    ddx = 0.32
-    At iterate      3    nfg =    88    J =      0.043086    ddx = 0.08
-    At iterate      4    nfg =   118    J =      0.040684    ddx = 0.02
-    At iterate      5    nfg =   152    J =      0.040604    ddx = 0.01
+    At iterate      0    nfg =     1    J =      0.695010    ddx = 0.64
+    At iterate      1    nfg =    30    J =      0.098411    ddx = 0.64
+    At iterate      2    nfg =    59    J =      0.045409    ddx = 0.32
+    At iterate      3    nfg =    88    J =      0.038182    ddx = 0.16
+    At iterate      4    nfg =   117    J =      0.037362    ddx = 0.08
+    At iterate      5    nfg =   150    J =      0.037087    ddx = 0.02
+    At iterate      6    nfg =   183    J =      0.036800    ddx = 0.02
+    At iterate      7    nfg =   216    J =      0.036763    ddx = 0.01
     CONVERGENCE: DDX < 0.01
 
 These lines show the different iterations of the optimization with information on the number of iterations, the number of cumulative evaluations ``nfg`` 
 (number of foward runs performed within each iteration of the optimization algorithm), the value of the cost function to minimize ``J`` and the value of the adaptive descent step ``ddx`` of this heuristic search algorihtm. 
-So, to summarize, the optimization algorithm has converged after 5 iterations by reaching the descent step tolerance criterion of 0.01. This optimization required to perform 152 forward run evaluations and leads to a final cost function value on the order of 0.04.
+So, to summarize, the optimization algorithm has converged after 7 iterations by reaching the descent step tolerance criterion of 0.01. This optimization required to perform 216 forward run evaluations and leads to a final cost function value on the order of 0.04.
 
 Then, we can ask which cost function ``J`` has been minimized and which parameters have been optimized. So, by default, the cost function to be minimized is one minus the Nash-Sutcliffe efficiency ``nse`` (:math:`1 - \text{NSE}`)
 and the optimized parameters are the set of rainfall-runoff parameters (``cp``, ``ct``, ``kexc`` and ``llr``). In the current configuration spatially
diff --git a/smash/core/model/model.py b/smash/core/model/model.py
index 4a8f5d1b..9017baea 100644
--- a/smash/core/model/model.py
+++ b/smash/core/model/model.py
@@ -1134,8 +1134,8 @@ def response(self) -> ResponseDT:
         >>> ind
         0
         >>> model.response.q[ind, :]
-        array([1.9826430e-03, 1.3466669e-07, 6.7617895e-12, ..., 3.2273201e+01,
-               3.2118713e+01, 3.1965160e+01], dtype=float32)
+        array([1.9826430e-03, 1.3466669e-07, 6.7617895e-12, ..., 2.2796249e+01,
+               2.2655941e+01, 2.2517307e+01], dtype=float32)
 
         If you are using IPython, tab completion allows you to visualize all the attributes and methods
 
@@ -1190,11 +1190,13 @@ def rr_final_states(self) -> RR_StatesDT:
         `get_rr_final_states <Model.get_rr_final_states>`
 
         >>> model.get_rr_final_states("hp")
-        array([[0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.5892611 , 0.60285664, 0.01      ,
+        array([[-99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                  0.8682228 ,   0.88014543, -99.        , -99.        ,
                 ...
-                0.01      , 0.01      , 0.01      ]], dtype=float32)
+                -99.        , -99.        , -99.        , -99.        ]],
+              dtype=float32)
 
         .. note::
             Unlike rainfall-runoff initial states, there is no setter for rainfall-runoff final states.
@@ -1999,11 +2001,13 @@ def get_rr_final_states(self, key: str) -> NDArray[np.float32]:
         Access to a specific rainfall-runoff final state grid
 
         >>> model.get_rr_final_states("hp")
-        array([[0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.5892611 , 0.60285664, 0.01      ,
+        array([[-99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                  0.8682228 ,   0.88014543, -99.        , -99.        ,
                 ...
-                0.01      , 0.01      , 0.01      ]], dtype=float32)
+                -99.        , -99.        , -99.        , -99.        ]],
+              dtype=float32)
 
         .. note::
             This method is equivalent to directly slicing the ``rr_final_states.values`` array (as shown
@@ -2023,11 +2027,13 @@ def get_rr_final_states(self, key: str) -> NDArray[np.float32]:
         Slice the ``rr_final_states.values`` array on the last axis
 
         >>> model.rr_final_states.values[..., ind]
-        array([[0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.01      , 0.01      , 0.01      ,
-                0.01      , 0.01      , 0.5892611 , 0.60285664, 0.01      ,
+        array([[-99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                -99.        , -99.        , -99.        , -99.        ,
+                  0.8682228 ,   0.88014543, -99.        , -99.        ,
                 ...
-                0.01      , 0.01      , 0.01      ]], dtype=float32)
+                -99.        , -99.        , -99.        , -99.        ]],
+              dtype=float32)
         """
 
         key = _standardize_get_rr_final_states_args(self, key)
@@ -2205,18 +2211,18 @@ def get_serr_mu(self) -> NDArray[np.float32]:
         Set arbitrary values to structural error mu parameters
 
         >>> model.set_serr_mu_parameters("mg0", 1)
-        >>> model.set_serr_mu_paramters("mg1", 2)
+        >>> model.set_serr_mu_parameters("mg1", 2)
 
         Retrieve the mu value with the `get_serr_mu <Model.get_serr_mu>` method
 
         >>> mu = model.get_serr_mu()
         >>> mu
-        array([[ 1.0039653,  1.0000002,  1.       , ..., 65.5464   , 65.23743  ,
-                64.93032  ],
-               [ 1.0004755,  1.       ,  1.       , ..., 16.804424 , 16.740883 ,
-                16.677677 ],
-               [ 1.0000595,  1.       ,  1.       , ...,  5.186602 ,  5.1694865,
-                5.1524224]], dtype=float32)
+        array([[ 1.0039653,  1.0000002,  1.       , ..., 46.5925   , 46.311882 ,
+                46.034615 ],
+               [ 1.0004755,  1.       ,  1.       , ..., 10.65963  , 10.61587  ,
+                10.572574 ],
+               [ 1.0000595,  1.       ,  1.       , ...,  3.563775 ,  3.5520396,
+                 3.5404253]], dtype=float32)
 
         This is equivalent to
 
@@ -2276,12 +2282,12 @@ def get_serr_sigma(self) -> NDArray[np.float32]:
 
         >>> sigma = model.get_serr_sigma()
         >>> sigma
-        array([[1.0003965, 1.       , 1.       , ..., 7.4546404, 7.423743 ,
-                7.393032 ],
-               [1.0000476, 1.       , 1.       , ..., 2.5804424, 2.5740883,
-                2.5677679],
-               [1.000006 , 1.       , 1.       , ..., 1.4186602, 1.4169487,
-                1.4152422]], dtype=float32)
+        array([[1.0003965, 1.       , 1.       , ..., 5.55925  , 5.5311885,
+                5.5034614],
+               [1.0000476, 1.       , 1.       , ..., 1.965963 , 1.9615871,
+                1.9572574],
+               [1.000006 , 1.       , 1.       , ..., 1.2563775, 1.255204 ,
+                1.2540425]], dtype=float32)
 
         This is equivalent to
 
diff --git a/smash/core/signal_analysis/metrics/metrics.py b/smash/core/signal_analysis/metrics/metrics.py
index 431ff9a3..b9770c03 100644
--- a/smash/core/signal_analysis/metrics/metrics.py
+++ b/smash/core/signal_analysis/metrics/metrics.py
@@ -62,23 +62,25 @@ def metrics(model: Model, metric: str = "nse", end_warmup: str | Timestamp | Non
 
     >>> model.optimize()
     </> Optimize
-        At iterate      0    nfg =     1    J =      0.643190    ddx = 0.64
-        At iterate      1    nfg =    30    J =      0.097397    ddx = 0.64
-        At iterate      2    nfg =    59    J =      0.052158    ddx = 0.32
-        At iterate      3    nfg =    88    J =      0.043086    ddx = 0.08
-        At iterate      4    nfg =   118    J =      0.040684    ddx = 0.02
-        At iterate      5    nfg =   152    J =      0.040604    ddx = 0.01
+        At iterate      0    nfg =     1    J =      0.695010    ddx = 0.64
+        At iterate      1    nfg =    30    J =      0.098411    ddx = 0.64
+        At iterate      2    nfg =    59    J =      0.045409    ddx = 0.32
+        At iterate      3    nfg =    88    J =      0.038182    ddx = 0.16
+        At iterate      4    nfg =   117    J =      0.037362    ddx = 0.08
+        At iterate      5    nfg =   150    J =      0.037087    ddx = 0.02
+        At iterate      6    nfg =   183    J =      0.036800    ddx = 0.02
+        At iterate      7    nfg =   216    J =      0.036763    ddx = 0.01
         CONVERGENCE: DDX < 0.01
 
     Compute the Mean Squared Error for all catchments
 
     >>> smash.metrics(model, metric="mse")
-    array([49.50750351,  4.31146955,  0.57386672])
+    array([44.78328323,  4.38410997,  0.50611502])
 
     The Kling-Gupta Efficiency
 
     >>> smash.metrics(model, metric="kge")
-    array([0.96001273, 0.86749017, 0.81862521])
+    array([0.94752783, 0.84582865, 0.8045246 ])
 
     Add end warm-up date (1 month)
 
@@ -89,9 +91,9 @@ def metrics(model: Model, metric: str = "nse", end_warmup: str | Timestamp | Non
     Compute the Nash-Sutcliffe Efficiency with and without warm-up
 
     >>> smash.metrics(model, metric="nse")
-    array([0.95939589, 0.90611488, 0.82942516])
+    array([0.96327233, 0.90453297, 0.84956211])
     >>> smash.metrics(model, metric="nse", end_warmup=end_warmup)
-    array([0.9774977 , 0.92978942, 0.87382156])
+    array([0.97294462, 0.91321117, 0.86573941])
     """
     metric, end_warmup = _standardize_metrics_args(metric, end_warmup, model.setup)
 
diff --git a/smash/core/signal_analysis/signatures/signatures.py b/smash/core/signal_analysis/signatures/signatures.py
index b40f82d0..891a626f 100644
--- a/smash/core/signal_analysis/signatures/signatures.py
+++ b/smash/core/signal_analysis/signatures/signatures.py
@@ -154,46 +154,46 @@ def signatures(
     Access simulated flood event signatures as a `pandas.DataFrame`
 
     >>> sign.event
-           code  season               start  ...    Erch2r  Elt        Epf
-    0  V3524010  autumn 2014-11-03 03:00:00  ...  0.214772  3.0  85.736832
-    1  V3515010  autumn 2014-11-03 10:00:00  ...  0.202139  0.0  17.256138
-    2  V3517010  autumn 2014-11-03 08:00:00  ...  0.187440  1.0   4.770674
+           code  season               start  ...    Erch2r  Elt         Epf
+    0  V3524010  autumn 2014-11-03 03:00:00  ...  0.324836  3.0  100.448631
+    1  V3515010  autumn 2014-11-03 10:00:00  ...  0.339471  0.0   20.168104
+    2  V3517010  autumn 2014-11-03 08:00:00  ...  0.319605  1.0    5.356519
     [3 rows x 12 columns]
 
     Access simulated continuous signatures as a `pandas.DataFrame`
 
     >>> sign.cont
-           code       Crc     Crchf  ...     Cfp10     Cfp50      Cfp90
-    0  V3524010  0.316813  0.058713  ...  0.013092  5.895140  39.439846
-    1  V3515010  0.251831  0.048244  ...  0.001120  1.072775   9.016475
-    2  V3517010  0.251551  0.048258  ...  0.000020  0.242764   2.433439
+           code       Crc     Crchf  ...         Cfp10     Cfp50      Cfp90
+    0  V3524010  0.208410  0.064490  ...  6.270661e-04  1.767450  27.806826
+    1  V3515010  0.139829  0.049087  ...  4.919724e-05  0.150414   5.161743
+    2  V3517010  0.139431  0.047900  ...  6.545733e-07  0.031947   1.355473
     [3 rows x 9 columns]
 
     Access specific simulated continuous signature for each gauge
 
     >>> sign.cont["Crc"]
-    0    0.516206
-    1    0.509181s
-    2    0.514302
+    0    0.208410
+    1    0.139829
+    2    0.139431
     Name: Crc, dtype: float64
 
     Access all simulated continuous signatures for a single gauge
 
     >>> sign.cont[sign.cont["code"] == "V3524010"]
-           code       Crc     Crchf  ...   Cfp10   Cfp50      Cfp90
-    0  V3524010  0.516206  0.191349  ...  1.1709  3.3225  42.631798
+           code      Crc    Crchf  ...     Cfp10    Cfp50      Cfp90
+    0  V3524010  0.20841  0.06449  ...  0.000627  1.76745  27.806826
     [1 rows x 9 columns]
 
     Access specific simulated continuous signatures for a single gauge
 
     >>> sign.cont["Crc"][sign.cont["code"] == "V3524010"]
-    0    0.516206
+    0    0.20841
     Name: Crc, dtype: float64
 
     Use the `item <pandas.Series.item>` method to return the underlying data as a Python scalar
 
     >>> sign.cont["Crc"][sign.cont["code"] == "V3524010"].item()
-    0.5162064433097839
+    0.20841039717197418
     """
 
     cs, es, domain, event_seg = _standardize_signatures_args(sign, domain, event_seg)
diff --git a/smash/core/simulation/_doc.py b/smash/core/simulation/_doc.py
index 39285459..df44f3bd 100644
--- a/smash/core/simulation/_doc.py
+++ b/smash/core/simulation/_doc.py
@@ -700,12 +700,12 @@ def _gen_docstring_from_base_doc(
 Get the simulated discharges
 
 >>> %(model_example_response)s.response.q
-array([[1.9826430e-03, 1.3466669e-07, 6.7617895e-12, ..., 3.2273201e+01,
-        3.2118713e+01, 3.1965160e+01],
-       [2.3777038e-04, 7.3761623e-09, 1.7551447e-13, ..., 7.9022121e+00,
-        7.8704414e+00, 7.8388391e+00],
-       [2.9721676e-05, 5.4272520e-10, 8.4623445e-15, ..., 2.0933011e+00,
-        2.0847433e+00, 2.0762112e+00]], dtype=float32)
+array([[1.9826430e-03, 1.3466669e-07, 6.7617895e-12, ..., 2.2796249e+01,
+        2.2655941e+01, 2.2517307e+01],
+       [2.3777038e-04, 7.3761623e-09, 1.7551447e-13, ..., 4.8298149e+00,
+        4.8079352e+00, 4.7862868e+00],
+       [2.9721676e-05, 5.4272520e-10, 8.4623445e-15, ..., 1.2818875e+00,
+        1.2760198e+00, 1.2702127e+00]], dtype=float32)
 """
 )
 
@@ -819,23 +819,25 @@ def _gen_docstring_from_base_doc(
 
 >>> %(model_example_func)s
 </> Optimize
-    At iterate      0    nfg =     1    J =      0.643190    ddx = 0.64
-    At iterate      1    nfg =    30    J =      0.097397    ddx = 0.64
-    At iterate      2    nfg =    59    J =      0.052158    ddx = 0.32
-    At iterate      3    nfg =    88    J =      0.043086    ddx = 0.08
-    At iterate      4    nfg =   118    J =      0.040684    ddx = 0.02
-    At iterate      5    nfg =   152    J =      0.040604    ddx = 0.01
+    At iterate      0    nfg =     1    J =      0.695010    ddx = 0.64
+    At iterate      1    nfg =    30    J =      0.098411    ddx = 0.64
+    At iterate      2    nfg =    59    J =      0.045409    ddx = 0.32
+    At iterate      3    nfg =    88    J =      0.038182    ddx = 0.16
+    At iterate      4    nfg =   117    J =      0.037362    ddx = 0.08
+    At iterate      5    nfg =   150    J =      0.037087    ddx = 0.02
+    At iterate      6    nfg =   183    J =      0.036800    ddx = 0.02
+    At iterate      7    nfg =   216    J =      0.036763    ddx = 0.01
     CONVERGENCE: DDX < 0.01
 
 Get the simulated discharges
 
 >>> %(model_example_response)s.response.q
-array([[1.9826430e-03, 1.3466669e-07, 6.7617895e-12, ..., 3.2273201e+01,
-        3.2118713e+01, 3.1965160e+01],
-       [2.3777038e-04, 7.3761623e-09, 1.7551447e-13, ..., 7.9022121e+00,
-        7.8704414e+00, 7.8388391e+00],
-       [2.9721676e-05, 5.4272520e-10, 8.4623445e-15, ..., 2.0933011e+00,
-        2.0847433e+00, 2.0762112e+00]], dtype=float32)
+array([[5.8217382e-04, 4.7552516e-04, 3.5390016e-04, ..., 1.9439360e+01,
+        1.9214035e+01, 1.8993553e+01],
+       [1.2144950e-04, 6.6219603e-05, 3.0706105e-05, ..., 4.8059664e+00,
+        4.7563825e+00, 4.7077618e+00],
+       [1.9631827e-05, 6.9778653e-06, 2.2202073e-06, ..., 1.2523955e+00,
+        1.2394531e+00, 1.2267693e+00]], dtype=float32)
 """
 )
 
@@ -926,12 +928,12 @@ def _gen_docstring_from_base_doc(
 Get the simulated discharges
 
 >>> %(model_example_response)s.response.q
-array([[9.4456947e-05, 9.3808041e-05, 9.3033530e-05, ..., 1.7851851e+01,
-        1.7691467e+01, 1.7528925e+01],
-       [2.5728115e-05, 2.4717448e-05, 2.3537395e-05, ..., 3.9668279e+00,
-        3.9301457e+00, 3.8938913e+00],
-       [5.9592148e-06, 5.0831463e-06, 4.2117308e-06, ..., 9.8274386e-01,
-        9.7398454e-01, 9.6534210e-01]], dtype=float32)
+array([[9.4571486e-05, 9.3920688e-05, 9.3143637e-05, ..., 1.7423288e+01,
+        1.7193638e+01, 1.6963835e+01],
+       [2.5758292e-05, 2.4744744e-05, 2.3561088e-05, ..., 3.6616585e+00,
+        3.6165960e+00, 3.5724759e+00],
+       [5.9654208e-06, 5.0872231e-06, 4.2139386e-06, ..., 9.0600485e-01,
+        8.9525825e-01, 8.8473159e-01]], dtype=float32)
 """
 )
 
@@ -1038,26 +1040,24 @@ def _gen_docstring_from_base_doc(
 
 >>> %(model_example_func)s
 </> Bayesian Optimize
-    At iterate      0    nfg =     1    J =     26.510803    ddx = 0.64
-    At iterate      1    nfg =    68    J =      2.536702    ddx = 0.64
-    At iterate      2    nfg =   136    J =      2.402311    ddx = 0.16
-    At iterate      3    nfg =   202    J =      2.329653    ddx = 0.16
-    At iterate      4    nfg =   270    J =      2.277469    ddx = 0.04
-    At iterate      5    nfg =   343    J =      2.271495    ddx = 0.02
-    At iterate      6    nfg =   416    J =      2.270596    ddx = 0.01
-    At iterate      7    nfg =   488    J =      2.269927    ddx = 0.01
-    At iterate      8    nfg =   561    J =      2.269505    ddx = 0.01
+    At iterate      0    nfg =     1    J =     77.049133    ddx = 0.64
+    At iterate      1    nfg =    68    J =      2.584603    ddx = 0.64
+    At iterate      2    nfg =   135    J =      2.324317    ddx = 0.32
+    At iterate      3    nfg =   202    J =      2.304130    ddx = 0.08
+    At iterate      4    nfg =   269    J =      2.262191    ddx = 0.08
+    At iterate      5    nfg =   343    J =      2.260251    ddx = 0.01
+    At iterate      6    nfg =   416    J =      2.258220    ddx = 0.00
     CONVERGENCE: DDX < 0.01
 
 Get the simulated discharges:
 
 >>> %(model_example_response)s.response.q
-array([[2.8263923e-04, 2.6972688e-04, 2.5154054e-04, ..., 2.2405909e+01,
-        2.2121796e+01, 2.1843649e+01],
-       [7.0431546e-05, 5.6386390e-05, 4.1599422e-05, ..., 5.3495941e+00,
-        5.2908549e+00, 5.2333012e+00],
-       [1.3507083e-05, 7.9051424e-06, 4.3033124e-06, ..., 1.3771975e+00,
-        1.3625814e+00, 1.3482267e+00]], dtype=float32)
+array([[3.8725851e-04, 3.5436003e-04, 3.0995562e-04, ..., 1.9623451e+01,
+        1.9391096e+01, 1.9163759e+01],
+       [9.0669761e-05, 6.3609077e-05, 3.9684928e-05, ..., 4.7896295e+00,
+        4.7395453e+00, 4.6904192e+00],
+       [1.6137006e-05, 7.8192916e-06, 3.4578904e-06, ..., 1.2418083e+00,
+        1.2288600e+00, 1.2161492e+00]], dtype=float32)
 """
 )
 
@@ -1128,7 +1128,7 @@ def _gen_docstring_from_base_doc(
 Get the cost values through multiple forward runs
 
 >>> mfr.cost
-array([1.17112  , 1.0390087, 1.2135248, 1.2491335, 1.2172333], dtype=float32)
+array([1.2170078, 1.0733036, 1.2239422, 1.2506678, 1.2261102], dtype=float32)
     """
 )
 
@@ -1250,7 +1250,7 @@ def _gen_docstring_from_base_doc(
 Get the cost values through multiple runs of optimization
 
 >>> mopt.cost
-array([0.5622911 , 0.0809496 , 0.16873538], dtype=float32)
+array([0.51374453, 0.0528878 , 0.15056956], dtype=float32)
 """
 )
 
@@ -1632,13 +1632,12 @@ def _gen_docstring_from_base_doc(
         },
     )
 </> Bayesian Optimize
-    At iterate      0    nfg =     1    J =     29.988619    ddx = 0.64
-    At iterate      1    nfg =    68    J =      2.749650    ddx = 0.64
-    At iterate      2    nfg =   135    J =      2.580721    ddx = 0.32
-    At iterate      3    nfg =   202    J =      2.473113    ddx = 0.16
-    At iterate      4    nfg =   269    J =      2.432616    ddx = 0.08
-    At iterate      5    nfg =   342    J =      2.430271    ddx = 0.02
-    At iterate      6    nfg =   416    J =      2.429243    ddx = 0.01
+    At iterate      0    nfg =     1    J =     80.526947    ddx = 0.64
+    At iterate      1    nfg =    68    J =      3.029253    ddx = 0.64
+    At iterate      2    nfg =   135    J =      2.764919    ddx = 0.32
+    At iterate      3    nfg =   203    J =      2.760564    ddx = 0.04
+    At iterate      4    nfg =   271    J =      2.755039    ddx = 0.02
+    At iterate      5    nfg =   344    J =      2.754198    ddx = 0.01
     CONVERGENCE: DDX < 0.01
 """
 )
diff --git a/smash/core/simulation/options.py b/smash/core/simulation/options.py
index 3ce6d44a..cee75531 100644
--- a/smash/core/simulation/options.py
+++ b/smash/core/simulation/options.py
@@ -102,12 +102,14 @@ def default_optimize_options(
 
     >>> model_u = smash.optimize(model, mapping="uniform", optimize_options=opt_u)
     </> Optimize
-        At iterate      0    nfg =     1    J =      0.643190    ddx = 0.64
-        At iterate      1    nfg =    30    J =      0.097397    ddx = 0.64
-        At iterate      2    nfg =    59    J =      0.052158    ddx = 0.32
-        At iterate      3    nfg =    88    J =      0.043086    ddx = 0.08
-        At iterate      4    nfg =   118    J =      0.040684    ddx = 0.02
-        At iterate      5    nfg =   152    J =      0.040604    ddx = 0.01
+        At iterate      0    nfg =     1    J =      0.695010    ddx = 0.64
+        At iterate      1    nfg =    30    J =      0.098411    ddx = 0.64
+        At iterate      2    nfg =    59    J =      0.045409    ddx = 0.32
+        At iterate      3    nfg =    88    J =      0.038182    ddx = 0.16
+        At iterate      4    nfg =   117    J =      0.037362    ddx = 0.08
+        At iterate      5    nfg =   150    J =      0.037087    ddx = 0.02
+        At iterate      6    nfg =   183    J =      0.036800    ddx = 0.02
+        At iterate      7    nfg =   216    J =      0.036763    ddx = 0.01
         CONVERGENCE: DDX < 0.01
 
     Customize the optimize options dictionary by removing ``'kexc'`` from the optimized parameters
@@ -158,20 +160,16 @@ def default_optimize_options(
 
     >>> model_u = smash.optimize(model, mapping="uniform", optimize_options=opt_u)
     </> Optimize
-        At iterate      0    nfg =     1    J =      0.643190    ddx = 0.64
-        At iterate      1    nfg =    17    J =      0.111053    ddx = 0.64
-        At iterate      2    nfg =    33    J =      0.081869    ddx = 0.32
-        At iterate      3    nfg =    49    J =      0.048932    ddx = 0.16
-        At iterate      4    nfg =    65    J =      0.048322    ddx = 0.04
-        At iterate      5    nfg =    85    J =      0.046548    ddx = 0.02
-        At iterate      6    nfg =   105    J =      0.046372    ddx = 0.01
-        At iterate      7    nfg =   123    J =      0.046184    ddx = 0.01
-        At iterate      8    nfg =   142    J =      0.046049    ddx = 0.01
-        At iterate      9    nfg =   161    J =      0.045862    ddx = 0.01
-        At iterate     10    nfg =   180    J =      0.045674    ddx = 0.01
-        At iterate     11    nfg =   200    J =      0.045624    ddx = 0.01
-        At iterate     12    nfg =   219    J =      0.045537    ddx = 0.00
-        CONVERGENCE: DDX < 0.01
+        At iterate      0    nfg =     1    J =      0.695010    ddx = 0.64
+        At iterate      1    nfg =    17    J =      0.128863    ddx = 0.64
+        At iterate      2    nfg =    32    J =      0.069483    ddx = 0.32
+        At iterate      3    nfg =    49    J =      0.045072    ddx = 0.16
+        At iterate      4    nfg =    65    J =      0.044047    ddx = 0.08
+        At iterate      5    nfg =    84    J =      0.043528    ddx = 0.04
+        At iterate      6    nfg =   102    J =      0.042690    ddx = 0.02
+        At iterate      7    nfg =   122    J =      0.042665    ddx = 0.01
+        At iterate      8    nfg =   140    J =      0.042606    ddx = 0.00
+        CONVERGENCE: DDX < 0.0
 
     Get the default optimization options for a different mapping
 
@@ -212,19 +210,19 @@ def default_optimize_options(
     >>> opt_ann["random_state"] = 11
     >>> model.optimize(mapping="ann", optimize_options=opt_ann)
     </> Optimize
-        At epoch      1    J =  1.150821    |proj g| =  0.002341
-        At epoch      2    J =  1.142061    |proj g| =  0.002610
+        At epoch      1    J =  1.208756    |proj g| =  0.001343
+        At epoch      2    J =  1.204409    |proj g| =  0.001567
         ...
-        At epoch     26    J =  0.103324    |proj g| =  0.005455
-        At epoch     27    J =  0.089769    |proj g| =  0.016618
-        At epoch     28    J =  0.104150    |proj g| =  0.019718
+        At epoch     25    J =  0.099455    |proj g| =  0.009175
+        At epoch     26    J =  0.078379    |proj g| =  0.013312
+        At epoch     27    J =  0.080410    |proj g| =  0.022210
         ...
-        At epoch     36    J =  0.185123    |proj g| =  0.011936
-        At epoch     37    J =  0.179911    |proj g| =  0.011819
-        Training:  92%|██████████████████████████▊  | 37/40 [00:30<00:02,  1.23it/s]
+        At epoch     35    J =  0.176823    |proj g| =  0.015846
+        At epoch     36    J =  0.172950    |proj g| =  0.015864
+        Training:  90%|████████████████▌     | 36/40 [00:21<00:02,  1.65it/s]
 
-    The training process was terminated after 37 epochs, where the loss did not decrease below the minimal
-    value at epoch 27 for 10 consecutive epochs. The optimal parameters are thus recorded at epoch 27
+    The training process was terminated after 36 epochs, where the loss did not decrease below the minimal
+    value at epoch 26 for 10 consecutive epochs. The optimal parameters are thus recorded at epoch 26
     """
     mapping, optimizer = _standardize_default_optimize_options_args(mapping, optimizer)
 
@@ -310,15 +308,13 @@ def default_bayesian_optimize_options(
 
     >>> model_u = smash.bayesian_optimize(model, mapping="uniform", optimize_options=opt_u)
     </> Bayesian Optimize
-        At iterate      0    nfg =     1    J =     26.510803    ddx = 0.64
-        At iterate      1    nfg =    68    J =      2.536702    ddx = 0.64
-        At iterate      2    nfg =   136    J =      2.402311    ddx = 0.16
-        At iterate      3    nfg =   202    J =      2.329653    ddx = 0.16
-        At iterate      4    nfg =   270    J =      2.277469    ddx = 0.04
-        At iterate      5    nfg =   343    J =      2.271495    ddx = 0.02
-        At iterate      6    nfg =   416    J =      2.270596    ddx = 0.01
-        At iterate      7    nfg =   488    J =      2.269927    ddx = 0.01
-        At iterate      8    nfg =   561    J =      2.269505    ddx = 0.01
+        At iterate      0    nfg =     1    J =     77.049133    ddx = 0.64
+        At iterate      1    nfg =    68    J =      2.584603    ddx = 0.64
+        At iterate      2    nfg =   135    J =      2.324317    ddx = 0.32
+        At iterate      3    nfg =   202    J =      2.304130    ddx = 0.08
+        At iterate      4    nfg =   269    J =      2.262191    ddx = 0.08
+        At iterate      5    nfg =   343    J =      2.260251    ddx = 0.01
+        At iterate      6    nfg =   416    J =      2.258220    ddx = 0.00
         CONVERGENCE: DDX < 0.01
 
     Get the default bayesian optimization options for a different mapping
@@ -351,19 +347,19 @@ def default_bayesian_optimize_options(
     >>> opt_ml["bounds"]["sg0"] = (1e-3, 100)
     >>> opt_ml["descriptor"]["cp"] = "slope"
     >>> opt_ml["termination_crit"]["maxiter"] = 10
-    >>> model.optimize(mapping="multi-linear", optimize_options=opt_ml)
+    >>> model.bayesian_optimize(mapping="multi-linear", optimize_options=opt_ml)
     </> Bayesian Optimize
-        At iterate      0    nfg =     1    J =     26.510801    |proj g| =     51.156681
-        At iterate      1    nfg =     2    J =      4.539526    |proj g| =      2.556828
-        At iterate      2    nfg =     4    J =      4.361516    |proj g| =      2.209623
-        At iterate      3    nfg =     5    J =      4.102887    |proj g| =      1.631524
-        At iterate      4    nfg =     6    J =      3.968790    |proj g| =      1.408150
-        At iterate      5    nfg =     7    J =      3.944855    |proj g| =      1.280883
-        At iterate      6    nfg =     8    J =      3.926999    |proj g| =      1.158015
-        At iterate      7    nfg =     9    J =      3.909628    |proj g| =      1.107825
-        At iterate      8    nfg =    10    J =      3.800612    |proj g| =      0.996645
-        At iterate      9    nfg =    11    J =      3.663789    |proj g| =      0.984840
-        At iterate     10    nfg =    12    J =      3.505583    |proj g| =      0.504640
+        At iterate      0    nfg =     1    J =     77.049095    |proj g| =    147.958771
+        At iterate      1    nfg =     2    J =      6.694370    |proj g| =      4.301311
+        At iterate      2    nfg =     3    J =      6.527157    |proj g| =      3.754591
+        At iterate      3    nfg =     4    J =      5.088758    |proj g| =      2.261085
+        At iterate      4    nfg =     5    J =      4.736641    |proj g| =      1.637916
+        At iterate      5    nfg =     6    J =      4.421250    |proj g| =      1.300914
+        At iterate      6    nfg =     7    J =      4.284939    |proj g| =      1.345769
+        At iterate      7    nfg =     8    J =      4.196455    |proj g| =      1.179024
+        At iterate      8    nfg =     9    J =      4.139528    |proj g| =      1.159659
+        At iterate      9    nfg =    10    J =      4.099973    |proj g| =      1.042880
+        At iterate     10    nfg =    11    J =      4.027408    |proj g| =      1.142047
         STOP: TOTAL NO. OF ITERATION EXCEEDS LIMIT
 
     The optimization process was terminated after 10 iterations, the maximal value we defined.
diff --git a/smash/fcore/forward/forward_db.f90 b/smash/fcore/forward/forward_db.f90
index 03cf3817..5b16a14a 100644
--- a/smash/fcore/forward/forward_db.f90
+++ b/smash/fcore/forward/forward_db.f90
@@ -12696,7 +12696,7 @@ SUBROUTINE GR4_TIME_STEP_D(setup, mesh, input_data, options, time_step&
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_d
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     REAL(sp) :: pn_d, en_d, pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d
     INTRINSIC MAX
     REAL(sp) :: temp
@@ -12706,12 +12706,15 @@ SUBROUTINE GR4_TIME_STEP_D(setup, mesh, input_data, options, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp_d = ac_mlt_d
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht&
-!$OMP&, ac_qt), SHARED(ac_prcp_d, ac_ci_d, ac_cp_d, ac_ct_d, ac_kexc_d, &
-!$OMP&ac_hi_d, ac_hp_d, ac_ht_d, ac_qt_d), PRIVATE(row, col, k, pn, en, &
-!$OMP&pr, perc, l, prr, prd, qr, qd), PRIVATE(pn_d, en_d, pr_d, perc_d, &
-!$OMP&l_d, prr_d, prd_d, qr_d, qd_d), PRIVATE(temp), SCHEDULE(static)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_hi, ac_hp&
+!$OMP&, ac_ht, ac_qt), SHARED(ac_prcp_d, ac_ci_d, ac_cp_d, ac_ct_d, &
+!$OMP&ac_kexc_d, ac_hi_d, ac_hp_d, ac_ht_d, ac_qt_d), PRIVATE(row, col, &
+!$OMP&k, pn, en, pr, perc, l, prr, prd, qr, qd), PRIVATE(pn_d, en_d, &
+!$OMP&pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d), PRIVATE(temp), &
+!$OMP&                                             SCHEDULE(static)
     DO col=1,mesh%ncol
       DO row=1,mesh%nrow
         IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%&
@@ -12722,8 +12725,8 @@ SUBROUTINE GR4_TIME_STEP_D(setup, mesh, input_data, options, time_step&
 &                            ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)&
 &                            , pn, pn_d, en, en_d)
             CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k&
-&                          ), 9._sp/4._sp, ac_hp(k), ac_hp_d(k), pr, &
-&                          pr_d, perc, perc_d)
+&                          ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc&
+&                          , perc_d)
             CALL GR_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), ac_ht(k), &
 &                        ac_ht_d(k), l, l_d)
           ELSE
@@ -12787,7 +12790,7 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, time_step&
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_b
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b
     INTRINSIC MAX
     INTEGER :: branch
@@ -12798,10 +12801,12 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, time_step&
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht&
-!$OMP&, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, &
-!$OMP&qd), PRIVATE(chunk_start, chunk_end)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_hi, ac_hp&
+!$OMP&, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, prr, &
+!$OMP&prd, qr, qd), PRIVATE(chunk_start, chunk_end)
     CALL GETSTATICSCHEDULE(1, mesh%ncol, 1, chunk_start, chunk_end)
     DO col=chunk_start,chunk_end
       DO row=1,mesh%nrow
@@ -12817,8 +12822,8 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, time_step&
             CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(&
 &                          k), pn, en)
             CALL PUSHREAL4(ac_hp(k))
-            CALL GR_PRODUCTION(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, &
+&                        perc)
             CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l)
             CALL PUSHCONTROL1B(1)
           ELSE
@@ -12848,12 +12853,12 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, time_step&
 !$OMP END PARALLEL
     ac_prcp_b = 0.0_4
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht&
-!$OMP&, ac_qt), SHARED(ac_prcp_b, ac_ci_b, ac_cp_b, ac_ct_b, ac_kexc_b, &
-!$OMP&ac_hi_b, ac_hp_b, ac_ht_b, ac_qt_b), PRIVATE(row, col, k, pn, en, &
-!$OMP&pr, perc, l, prr, prd, qr, qd), PRIVATE(pn_b, en_b, pr_b, perc_b, &
-!$OMP&l_b, prr_b, prd_b, qr_b, qd_b), PRIVATE(branch, chunk_end, &
-!$OMP&chunk_start)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_hi, ac_hp&
+!$OMP&, ac_ht, ac_qt), SHARED(ac_prcp_b, ac_ci_b, ac_cp_b, ac_ct_b, &
+!$OMP&ac_kexc_b, ac_hi_b, ac_hp_b, ac_ht_b, ac_qt_b), PRIVATE(row, col, &
+!$OMP&k, pn, en, pr, perc, l, prr, prd, qr, qd), PRIVATE(pn_b, en_b, &
+!$OMP&pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b), PRIVATE(branch, &
+!$OMP&chunk_end, chunk_start)
     CALL POPREAL4(en)
     CALL POPREAL4(prr)
     CALL POPREAL4(pn)
@@ -12898,8 +12903,8 @@ SUBROUTINE GR4_TIME_STEP_B(setup, mesh, input_data, options, time_step&
 &                        ac_ht_b(k), l, l_b)
             CALL POPREAL4(ac_hp(k))
             CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k&
-&                          ), 9._sp/4._sp, ac_hp(k), ac_hp_b(k), pr, &
-&                          pr_b, perc, perc_b)
+&                          ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc&
+&                          , perc_b)
             CALL POPREAL4(ac_hi(k))
             CALL POPREAL4(pn)
             CALL POPREAL4(en)
@@ -12929,17 +12934,19 @@ SUBROUTINE GR4_TIME_STEP(setup, mesh, input_data, options, time_step, &
     REAL(sp), DIMENSION(mesh%nac), INTENT(INOUT) :: ac_qt
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     INTRINSIC MAX
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'prcp', ac_prcp)
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht&
-!$OMP&, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, &
-!$OMP&qd), SCHEDULE(static)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_hi, ac_hp&
+!$OMP&, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, prr, &
+!$OMP&prd, qr, qd), SCHEDULE(static)
     DO col=1,mesh%ncol
       DO row=1,mesh%nrow
         IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%&
@@ -12948,8 +12955,8 @@ SUBROUTINE GR4_TIME_STEP(setup, mesh, input_data, options, time_step, &
           IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN
             CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(&
 &                          k), pn, en)
-            CALL GR_PRODUCTION(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, &
+&                        perc)
             CALL GR_EXCHANGE(ac_kexc(k), ac_ht(k), l)
           ELSE
             pr = 0._sp
@@ -13002,7 +13009,7 @@ SUBROUTINE GR5_TIME_STEP_D(setup, mesh, input_data, options, time_step&
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_d
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     REAL(sp) :: pn_d, en_d, pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d
     INTRINSIC MAX
     REAL(sp) :: temp
@@ -13012,11 +13019,13 @@ SUBROUTINE GR5_TIME_STEP_D(setup, mesh, input_data, options, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp_d = ac_mlt_d
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_aexc, ac_hi, &
-!$OMP&ac_hp, ac_ht, ac_qt), SHARED(ac_prcp_d, ac_ci_d, ac_cp_d, ac_ct_d&
-!$OMP&, ac_kexc_d, ac_aexc_d, ac_hi_d, ac_hp_d, ac_ht_d, ac_qt_d), &
-!$OMP&PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd), &
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_aexc, &
+!$OMP&ac_hi, ac_hp, ac_ht, ac_qt), SHARED(ac_prcp_d, ac_ci_d, ac_cp_d, &
+!$OMP&ac_ct_d, ac_kexc_d, ac_aexc_d, ac_hi_d, ac_hp_d, ac_ht_d, ac_qt_d)&
+!$OMP&, PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd), &
 !$OMP&PRIVATE(pn_d, en_d, pr_d, perc_d, l_d, prr_d, prd_d, qr_d, qd_d), &
 !$OMP&PRIVATE(temp), SCHEDULE(static)
     DO col=1,mesh%ncol
@@ -13029,8 +13038,8 @@ SUBROUTINE GR5_TIME_STEP_D(setup, mesh, input_data, options, time_step&
 &                            ac_ci(k), ac_ci_d(k), ac_hi(k), ac_hi_d(k)&
 &                            , pn, pn_d, en, en_d)
             CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_cp(k), ac_cp_d(k&
-&                          ), 9._sp/4._sp, ac_hp(k), ac_hp_d(k), pr, &
-&                          pr_d, perc, perc_d)
+&                          ), beta, ac_hp(k), ac_hp_d(k), pr, pr_d, perc&
+&                          , perc_d)
             CALL GR_THRESHOLD_EXCHANGE_D(ac_kexc(k), ac_kexc_d(k), &
 &                                  ac_aexc(k), ac_aexc_d(k), ac_ht(k), &
 &                                  ac_ht_d(k), l, l_d)
@@ -13095,7 +13104,7 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, time_step&
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_b
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     REAL(sp) :: pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b
     INTRINSIC MAX
     INTEGER :: branch
@@ -13106,10 +13115,12 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, time_step&
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_aexc, ac_hi, &
-!$OMP&ac_hp, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, &
-!$OMP&prr, prd, qr, qd), PRIVATE(chunk_start, chunk_end)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_aexc, &
+!$OMP&ac_hi, ac_hp, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc&
+!$OMP&, l, prr, prd, qr, qd), PRIVATE(chunk_start, chunk_end)
     CALL GETSTATICSCHEDULE(1, mesh%ncol, 1, chunk_start, chunk_end)
     DO col=chunk_start,chunk_end
       DO row=1,mesh%nrow
@@ -13125,8 +13136,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, time_step&
             CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(&
 &                          k), pn, en)
             CALL PUSHREAL4(ac_hp(k))
-            CALL GR_PRODUCTION(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, &
+&                        perc)
             CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)&
 &                                , l)
             CALL PUSHCONTROL1B(1)
@@ -13157,10 +13168,10 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, time_step&
 !$OMP END PARALLEL
     ac_prcp_b = 0.0_4
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_aexc, ac_hi, &
-!$OMP&ac_hp, ac_ht, ac_qt), SHARED(ac_prcp_b, ac_ci_b, ac_cp_b, ac_ct_b&
-!$OMP&, ac_kexc_b, ac_aexc_b, ac_hi_b, ac_hp_b, ac_ht_b, ac_qt_b), &
-!$OMP&PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd), &
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_aexc, &
+!$OMP&ac_hi, ac_hp, ac_ht, ac_qt), SHARED(ac_prcp_b, ac_ci_b, ac_cp_b, &
+!$OMP&ac_ct_b, ac_kexc_b, ac_aexc_b, ac_hi_b, ac_hp_b, ac_ht_b, ac_qt_b)&
+!$OMP&, PRIVATE(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd), &
 !$OMP&PRIVATE(pn_b, en_b, pr_b, perc_b, l_b, prr_b, prd_b, qr_b, qd_b), &
 !$OMP&PRIVATE(branch, chunk_end, chunk_start)
     CALL POPREAL4(en)
@@ -13208,8 +13219,8 @@ SUBROUTINE GR5_TIME_STEP_B(setup, mesh, input_data, options, time_step&
 &                                  ac_ht_b(k), l, l_b)
             CALL POPREAL4(ac_hp(k))
             CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_cp(k), ac_cp_b(k&
-&                          ), 9._sp/4._sp, ac_hp(k), ac_hp_b(k), pr, &
-&                          pr_b, perc, perc_b)
+&                          ), beta, ac_hp(k), ac_hp_b(k), pr, pr_b, perc&
+&                          , perc_b)
             CALL POPREAL4(ac_hi(k))
             CALL POPREAL4(pn)
             CALL POPREAL4(en)
@@ -13240,17 +13251,19 @@ SUBROUTINE GR5_TIME_STEP(setup, mesh, input_data, options, time_step, &
     REAL(sp), DIMENSION(mesh%nac), INTENT(INOUT) :: ac_qt
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     INTEGER :: row, col, k
-    REAL(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+    REAL(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
     INTRINSIC MAX
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'prcp', ac_prcp)
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_aexc, ac_hi, &
-!$OMP&ac_hp, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc, l, &
-!$OMP&prr, prd, qr, qd), SCHEDULE(static)
+!$OMP&ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_aexc, &
+!$OMP&ac_hi, ac_hp, ac_ht, ac_qt), PRIVATE(row, col, k, pn, en, pr, perc&
+!$OMP&, l, prr, prd, qr, qd), SCHEDULE(static)
     DO col=1,mesh%ncol
       DO row=1,mesh%nrow
         IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%&
@@ -13259,8 +13272,8 @@ SUBROUTINE GR5_TIME_STEP(setup, mesh, input_data, options, time_step, &
           IF (ac_prcp(k) .GE. 0._sp .AND. ac_pet(k) .GE. 0._sp) THEN
             CALL GR_INTERCEPTION(ac_prcp(k), ac_pet(k), ac_ci(k), ac_hi(&
 &                          k), pn, en)
-            CALL GR_PRODUCTION(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_cp(k), beta, ac_hp(k), pr, &
+&                        perc)
             CALL GR_THRESHOLD_EXCHANGE(ac_kexc(k), ac_aexc(k), ac_ht(k)&
 &                                , l)
           ELSE
@@ -13611,7 +13624,7 @@ SUBROUTINE LOIEAU_TIME_STEP_D(setup, mesh, input_data, options, &
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_d
     INTEGER :: row, col, k
-    REAL(sp) :: ei, pn, en, pr, perc, prr, prd, qr, qd
+    REAL(sp) :: beta, ei, pn, en, pr, perc, prr, prd, qr, qd
     REAL(sp) :: ei_d, pn_d, en_d, pr_d, perc_d, prr_d, prd_d, qr_d, qd_d
     INTRINSIC MIN
     INTRINSIC MAX
@@ -13622,12 +13635,14 @@ SUBROUTINE LOIEAU_TIME_STEP_D(setup, mesh, input_data, options, &
 &                              , 'pet', ac_pet)
     ac_prcp_d = ac_mlt_d
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ca, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), SHARED&
-!$OMP&(ac_prcp_d, ac_ca_d, ac_cc_d, ac_kb_d, ac_ha_d, ac_hc_d, ac_qt_d)&
-!$OMP&, PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr, qd), &
-!$OMP&PRIVATE(ei_d, pn_d, en_d, pr_d, perc_d, prr_d, prd_d, qr_d, qd_d)&
-!$OMP&, PRIVATE(temp), SCHEDULE(static)
+!$OMP&ac_prcp, ac_pet, ac_ca, beta, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
+!$OMP&SHARED(ac_prcp_d, ac_ca_d, ac_cc_d, ac_kb_d, ac_ha_d, ac_hc_d, &
+!$OMP&ac_qt_d), PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr&
+!$OMP&, qd), PRIVATE(ei_d, pn_d, en_d, pr_d, perc_d, prr_d, prd_d, qr_d&
+!$OMP&, qd_d), PRIVATE(temp), SCHEDULE(static)
     DO col=1,mesh%ncol
       DO row=1,mesh%nrow
         IF (.NOT.(mesh%active_cell(row, col) .EQ. 0 .OR. mesh%&
@@ -13651,8 +13666,8 @@ SUBROUTINE LOIEAU_TIME_STEP_D(setup, mesh, input_data, options, &
             en_d = -ei_d
             en = ac_pet(k) - ei
             CALL GR_PRODUCTION_D(pn, pn_d, en, en_d, ac_ca(k), ac_ca_d(k&
-&                          ), 9._sp/4._sp, ac_ha(k), ac_ha_d(k), pr, &
-&                          pr_d, perc, perc_d)
+&                          ), beta, ac_ha(k), ac_ha_d(k), pr, pr_d, perc&
+&                          , perc_d)
           ELSE
             pr = 0._sp
             perc = 0._sp
@@ -13708,7 +13723,7 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, &
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp_b
     INTEGER :: row, col, k
-    REAL(sp) :: ei, pn, en, pr, perc, prr, prd, qr, qd
+    REAL(sp) :: beta, ei, pn, en, pr, perc, prr, prd, qr, qd
     REAL(sp) :: ei_b, pn_b, en_b, pr_b, perc_b, prr_b, prd_b, qr_b, qd_b
     INTRINSIC MIN
     INTRINSIC MAX
@@ -13720,8 +13735,10 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, &
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ca, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
+!$OMP&ac_prcp, ac_pet, ac_ca, beta, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
 !$OMP&PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr, qd), &
 !$OMP&PRIVATE(chunk_start, chunk_end)
     CALL GETSTATICSCHEDULE(1, mesh%ncol, 1, chunk_start, chunk_end)
@@ -13752,8 +13769,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, &
             CALL PUSHREAL4(en)
             en = ac_pet(k) - ei
             CALL PUSHREAL4(ac_ha(k))
-            CALL GR_PRODUCTION(pn, en, ac_ca(k), 9._sp/4._sp, ac_ha(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, &
+&                        perc)
             CALL PUSHCONTROL1B(1)
           ELSE
             CALL PUSHCONTROL1B(0)
@@ -13788,11 +13805,11 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, &
 !$OMP END PARALLEL
     ac_prcp_b = 0.0_4
 !$OMP PARALLEL NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ca, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), SHARED&
-!$OMP&(ac_prcp_b, ac_ca_b, ac_cc_b, ac_kb_b, ac_ha_b, ac_hc_b, ac_qt_b)&
-!$OMP&, PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr, qd), &
-!$OMP&PRIVATE(ei_b, pn_b, en_b, pr_b, perc_b, prr_b, prd_b, qr_b, qd_b)&
-!$OMP&, PRIVATE(branch, chunk_end, chunk_start)
+!$OMP&ac_prcp, ac_pet, ac_ca, beta, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
+!$OMP&SHARED(ac_prcp_b, ac_ca_b, ac_cc_b, ac_kb_b, ac_ha_b, ac_hc_b, &
+!$OMP&ac_qt_b), PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr&
+!$OMP&, qd), PRIVATE(ei_b, pn_b, en_b, pr_b, perc_b, prr_b, prd_b, qr_b&
+!$OMP&, qd_b), PRIVATE(branch, chunk_end, chunk_start)
     CALL POPREAL4(qd)
     CALL POPREAL4(qr)
     CALL POPREAL4(en)
@@ -13839,8 +13856,8 @@ SUBROUTINE LOIEAU_TIME_STEP_B(setup, mesh, input_data, options, &
           IF (branch .NE. 0) THEN
             CALL POPREAL4(ac_ha(k))
             CALL GR_PRODUCTION_B(pn, pn_b, en, en_b, ac_ca(k), ac_ca_b(k&
-&                          ), 9._sp/4._sp, ac_ha(k), ac_ha_b(k), pr, &
-&                          pr_b, perc, perc_b)
+&                          ), beta, ac_ha(k), ac_ha_b(k), pr, pr_b, perc&
+&                          , perc_b)
             CALL POPREAL4(en)
             ei_b = -en_b
             CALL POPCONTROL1B(branch)
@@ -13879,7 +13896,7 @@ SUBROUTINE LOIEAU_TIME_STEP(setup, mesh, input_data, options, &
     REAL(sp), DIMENSION(mesh%nac), INTENT(INOUT) :: ac_qt
     REAL(sp), DIMENSION(mesh%nac) :: ac_prcp, ac_pet
     INTEGER :: row, col, k
-    REAL(sp) :: ei, pn, en, pr, perc, prr, prd, qr, qd
+    REAL(sp) :: beta, ei, pn, en, pr, perc, prr, prd, qr, qd
     INTRINSIC MIN
     INTRINSIC MAX
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
@@ -13887,8 +13904,10 @@ SUBROUTINE LOIEAU_TIME_STEP(setup, mesh, input_data, options, &
     CALL GET_AC_ATMOS_DATA_TIME_STEP(setup, mesh, input_data, time_step&
 &                              , 'pet', ac_pet)
     ac_prcp = ac_prcp + ac_mlt
+! Beta percolation parameter is time step dependent
+    beta = 9._sp/4._sp*(86400._sp/setup%dt)**0.25_sp
 !$OMP PARALLEL DO NUM_THREADS(options%comm%ncpu), SHARED(setup, mesh, &
-!$OMP&ac_prcp, ac_pet, ac_ca, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
+!$OMP&ac_prcp, ac_pet, ac_ca, beta, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt), &
 !$OMP&PRIVATE(row, col, k, ei, pn, en, pr, perc, prr, prd, qr, qd), &
 !$OMP&                                              SCHEDULE(static)
     DO col=1,mesh%ncol
@@ -13908,8 +13927,8 @@ SUBROUTINE LOIEAU_TIME_STEP(setup, mesh, input_data, options, &
               pn = 0._sp
             END IF
             en = ac_pet(k) - ei
-            CALL GR_PRODUCTION(pn, en, ac_ca(k), 9._sp/4._sp, ac_ha(k), &
-&                        pr, perc)
+            CALL GR_PRODUCTION(pn, en, ac_ca(k), beta, ac_ha(k), pr, &
+&                        perc)
           ELSE
             pr = 0._sp
             perc = 0._sp
diff --git a/smash/fcore/operator/md_gr_operator.f90 b/smash/fcore/operator/md_gr_operator.f90
index 43ff841b..63e55dfa 100644
--- a/smash/fcore/operator/md_gr_operator.f90
+++ b/smash/fcore/operator/md_gr_operator.f90
@@ -151,15 +151,18 @@ subroutine gr4_time_step(setup, mesh, input_data, options, time_step, ac_mlt, ac
 
         real(sp), dimension(mesh%nac) :: ac_prcp, ac_pet
         integer :: row, col, k
-        real(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+        real(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
 
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "prcp", ac_prcp)
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "pet", ac_pet)
 
         ac_prcp = ac_prcp + ac_mlt
 
+        ! Beta percolation parameter is time step dependent
+        beta = (9._sp/4._sp)*(86400._sp/setup%dt)**0.25_sp
+
         !$OMP parallel do schedule(static) num_threads(options%comm%ncpu) &
-        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht, &
+        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_hi, ac_hp, ac_ht, &
         !$OMP& ac_qt) &
         !$OMP& private(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd)
         do col = 1, mesh%ncol
@@ -174,7 +177,7 @@ subroutine gr4_time_step(setup, mesh, input_data, options, time_step, ac_mlt, ac
                     call gr_interception(ac_prcp(k), ac_pet(k), ac_ci(k), &
                     & ac_hi(k), pn, en)
 
-                    call gr_production(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), pr, perc)
+                    call gr_production(pn, en, ac_cp(k), beta, ac_hp(k), pr, perc)
 
                     call gr_exchange(ac_kexc(k), ac_ht(k), l)
 
@@ -221,16 +224,19 @@ subroutine gr5_time_step(setup, mesh, input_data, options, time_step, ac_mlt, ac
 
         real(sp), dimension(mesh%nac) :: ac_prcp, ac_pet
         integer :: row, col, k
-        real(sp) :: pn, en, pr, perc, l, prr, prd, qr, qd
+        real(sp) :: beta, pn, en, pr, perc, l, prr, prd, qr, qd
 
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "prcp", ac_prcp)
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "pet", ac_pet)
 
         ac_prcp = ac_prcp + ac_mlt
 
+        ! Beta percolation parameter is time step dependent
+        beta = (9._sp/4._sp)*(86400._sp/setup%dt)**0.25_sp
+
         !$OMP parallel do schedule(static) num_threads(options%comm%ncpu) &
-        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ci, ac_cp, ac_ct, ac_kexc, ac_aexc, ac_hi, ac_hp, &
-        !$OMP& ac_ht, ac_qt) &
+        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ci, ac_cp, beta, ac_ct, ac_kexc, ac_aexc, ac_hi, &
+        !$OMP& ac_hp, ac_ht, ac_qt) &
         !$OMP& private(row, col, k, pn, en, pr, perc, l, prr, prd, qr, qd)
         do col = 1, mesh%ncol
             do row = 1, mesh%nrow
@@ -244,7 +250,7 @@ subroutine gr5_time_step(setup, mesh, input_data, options, time_step, ac_mlt, ac
                     call gr_interception(ac_prcp(k), ac_pet(k), ac_ci(k), &
                     & ac_hi(k), pn, en)
 
-                    call gr_production(pn, en, ac_cp(k), 9._sp/4._sp, ac_hp(k), pr, perc)
+                    call gr_production(pn, en, ac_cp(k), beta, ac_hp(k), pr, perc)
 
                     call gr_threshold_exchange(ac_kexc(k), ac_aexc(k), ac_ht(k), l)
 
@@ -357,15 +363,18 @@ subroutine loieau_time_step(setup, mesh, input_data, options, time_step, ac_mlt,
 
         real(sp), dimension(mesh%nac) :: ac_prcp, ac_pet
         integer :: row, col, k
-        real(sp) :: ei, pn, en, pr, perc, prr, prd, qr, qd
+        real(sp) :: beta, ei, pn, en, pr, perc, prr, prd, qr, qd
 
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "prcp", ac_prcp)
         call get_ac_atmos_data_time_step(setup, mesh, input_data, time_step, "pet", ac_pet)
 
         ac_prcp = ac_prcp + ac_mlt
 
+        ! Beta percolation parameter is time step dependent
+        beta = (9._sp/4._sp)*(86400._sp/setup%dt)**0.25_sp
+
         !$OMP parallel do schedule(static) num_threads(options%comm%ncpu) &
-        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ca, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt) &
+        !$OMP& shared(setup, mesh, ac_prcp, ac_pet, ac_ca, beta, ac_cc, ac_kb, ac_ha, ac_hc, ac_qt) &
         !$OMP& private(row, col, k, ei, pn, en, pr, perc, prr, prd, qr, qd)
         do col = 1, mesh%ncol
             do row = 1, mesh%nrow
@@ -382,7 +391,7 @@ subroutine loieau_time_step(setup, mesh, input_data, options, time_step, ac_mlt,
 
                     en = ac_pet(k) - ei
 
-                    call gr_production(pn, en, ac_ca(k), 9._sp/4._sp, ac_ha(k), pr, perc)
+                    call gr_production(pn, en, ac_ca(k), beta, ac_ha(k), pr, perc)
 
                 else
 
diff --git a/smash/tests/baseline.hdf5 b/smash/tests/baseline.hdf5
index 9a9cb031..6748cc58 100644
Binary files a/smash/tests/baseline.hdf5 and b/smash/tests/baseline.hdf5 differ
diff --git a/smash/tests/core/simulation/test_estimate.py b/smash/tests/core/simulation/test_estimate.py
index e241d532..d74614f9 100644
--- a/smash/tests/core/simulation/test_estimate.py
+++ b/smash/tests/core/simulation/test_estimate.py
@@ -16,10 +16,10 @@ def generic_multiset_estimate(model: smash.Model, **kwargs) -> dict:
     problem = {
         "num_vars": 3,
         "names": ["cp", "kexc", "hp"],
-        "bounds": [(1, 500), (-50, 50), (0, 1)],
+        "bounds": [(20, 1000), (-20, 5), (0.2, 0.8)],
     }
 
-    sample = smash.factory.generate_samples(problem, n=20, random_state=11)
+    sample = smash.factory.generate_samples(problem, n=30, random_state=11)
 
     multisets = {}
 
diff --git a/smash/tests/diff_baseline.csv b/smash/tests/diff_baseline.csv
index e6457974..963654e7 100644
--- a/smash/tests/diff_baseline.csv
+++ b/smash/tests/diff_baseline.csv
@@ -1,11 +1,8 @@
-commit 645bfc7bb0f197c3d511bb7e7818babed01668e5
-Merge: ebd49e4 c8267f0
-Author: Francois Colleoni <110899888+inoelloc@users.noreply.github.com>
-Date:   Tue Mar 5 11:05:45 2024 +0100
+commit 307eb806a3ba38a0b22e3318a42de96a94793400
+Author: Francois Colleoni <francois.colleoni@inrae.fr>
+Date:   Thu Mar 28 16:26:18 2024 +0100
 
-    Merge pull request #146 from DassHydro/fix-read-qobs-date-edge
-    
-    FIX: Fix a bug when the csv file of qobs start or end with start_time or end_time
+    FIX PR: Change test multiset estimate setup
 
 TEST NAME                                                 |STATUS
 bbox_mesh.active_cell                                     |NON MODIFIED
@@ -26,101 +23,101 @@ bbox_mesh.xmin                                            |NON MODIFIED
 bbox_mesh.xres                                            |NON MODIFIED
 bbox_mesh.ymax                                            |NON MODIFIED
 bbox_mesh.yres                                            |NON MODIFIED
-custom_bayesian_optimize.zero-gr4-lr.custom_set_1.sim_q   |NON MODIFIED
-custom_bayesian_optimize.zero-gr4-lr.custom_set_2.sim_q   |NON MODIFIED
-custom_bayesian_optimize.zero-gr4-lr.custom_set_3.sim_q   |NON MODIFIED
-custom_bayesian_optimize.zero-gr4-lr.custom_set_4.sim_q   |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_1.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_2.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_3.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_4.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_5.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_6.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_7.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_8.sim_q            |NON MODIFIED
-custom_optimize.zero-gr4-lr.custom_set_9.sim_q            |NON MODIFIED
-forward_run.zero-gr4-kw.cost                              |NON MODIFIED
-forward_run.zero-gr4-kw.jobs                              |NON MODIFIED
-forward_run.zero-gr4-kw.q_domain                          |NON MODIFIED
+custom_bayesian_optimize.zero-gr4-lr.custom_set_1.sim_q   |MODIFIED
+custom_bayesian_optimize.zero-gr4-lr.custom_set_2.sim_q   |MODIFIED
+custom_bayesian_optimize.zero-gr4-lr.custom_set_3.sim_q   |MODIFIED
+custom_bayesian_optimize.zero-gr4-lr.custom_set_4.sim_q   |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_1.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_2.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_3.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_4.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_5.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_6.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_7.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_8.sim_q            |MODIFIED
+custom_optimize.zero-gr4-lr.custom_set_9.sim_q            |MODIFIED
+forward_run.zero-gr4-kw.cost                              |MODIFIED
+forward_run.zero-gr4-kw.jobs                              |MODIFIED
+forward_run.zero-gr4-kw.q_domain                          |MODIFIED
 forward_run.zero-gr4-kw.rr_states.hi                      |NON MODIFIED
-forward_run.zero-gr4-kw.rr_states.hp                      |NON MODIFIED
-forward_run.zero-gr4-kw.rr_states.ht                      |NON MODIFIED
-forward_run.zero-gr4-kw.sim_q                             |NON MODIFIED
-forward_run.zero-gr4-lag0.cost                            |NON MODIFIED
-forward_run.zero-gr4-lag0.jobs                            |NON MODIFIED
-forward_run.zero-gr4-lag0.q_domain                        |NON MODIFIED
+forward_run.zero-gr4-kw.rr_states.hp                      |MODIFIED
+forward_run.zero-gr4-kw.rr_states.ht                      |MODIFIED
+forward_run.zero-gr4-kw.sim_q                             |MODIFIED
+forward_run.zero-gr4-lag0.cost                            |MODIFIED
+forward_run.zero-gr4-lag0.jobs                            |MODIFIED
+forward_run.zero-gr4-lag0.q_domain                        |MODIFIED
 forward_run.zero-gr4-lag0.rr_states.hi                    |NON MODIFIED
-forward_run.zero-gr4-lag0.rr_states.hp                    |NON MODIFIED
-forward_run.zero-gr4-lag0.rr_states.ht                    |NON MODIFIED
-forward_run.zero-gr4-lag0.sim_q                           |NON MODIFIED
-forward_run.zero-gr4-lr.cost                              |NON MODIFIED
-forward_run.zero-gr4-lr.jobs                              |NON MODIFIED
-forward_run.zero-gr4-lr.q_domain                          |NON MODIFIED
+forward_run.zero-gr4-lag0.rr_states.hp                    |MODIFIED
+forward_run.zero-gr4-lag0.rr_states.ht                    |MODIFIED
+forward_run.zero-gr4-lag0.sim_q                           |MODIFIED
+forward_run.zero-gr4-lr.cost                              |MODIFIED
+forward_run.zero-gr4-lr.jobs                              |MODIFIED
+forward_run.zero-gr4-lr.q_domain                          |MODIFIED
 forward_run.zero-gr4-lr.rr_states.hi                      |NON MODIFIED
 forward_run.zero-gr4-lr.rr_states.hlr                     |NON MODIFIED
-forward_run.zero-gr4-lr.rr_states.hp                      |NON MODIFIED
-forward_run.zero-gr4-lr.rr_states.ht                      |NON MODIFIED
-forward_run.zero-gr4-lr.sim_q                             |NON MODIFIED
-forward_run.zero-gr5-kw.cost                              |NON MODIFIED
-forward_run.zero-gr5-kw.jobs                              |NON MODIFIED
-forward_run.zero-gr5-kw.q_domain                          |NON MODIFIED
+forward_run.zero-gr4-lr.rr_states.hp                      |MODIFIED
+forward_run.zero-gr4-lr.rr_states.ht                      |MODIFIED
+forward_run.zero-gr4-lr.sim_q                             |MODIFIED
+forward_run.zero-gr5-kw.cost                              |MODIFIED
+forward_run.zero-gr5-kw.jobs                              |MODIFIED
+forward_run.zero-gr5-kw.q_domain                          |MODIFIED
 forward_run.zero-gr5-kw.rr_states.hi                      |NON MODIFIED
-forward_run.zero-gr5-kw.rr_states.hp                      |NON MODIFIED
-forward_run.zero-gr5-kw.rr_states.ht                      |NON MODIFIED
-forward_run.zero-gr5-kw.sim_q                             |NON MODIFIED
-forward_run.zero-gr5-lag0.cost                            |NON MODIFIED
-forward_run.zero-gr5-lag0.jobs                            |NON MODIFIED
-forward_run.zero-gr5-lag0.q_domain                        |NON MODIFIED
+forward_run.zero-gr5-kw.rr_states.hp                      |MODIFIED
+forward_run.zero-gr5-kw.rr_states.ht                      |MODIFIED
+forward_run.zero-gr5-kw.sim_q                             |MODIFIED
+forward_run.zero-gr5-lag0.cost                            |MODIFIED
+forward_run.zero-gr5-lag0.jobs                            |MODIFIED
+forward_run.zero-gr5-lag0.q_domain                        |MODIFIED
 forward_run.zero-gr5-lag0.rr_states.hi                    |NON MODIFIED
-forward_run.zero-gr5-lag0.rr_states.hp                    |NON MODIFIED
-forward_run.zero-gr5-lag0.rr_states.ht                    |NON MODIFIED
-forward_run.zero-gr5-lag0.sim_q                           |NON MODIFIED
-forward_run.zero-gr5-lr.cost                              |NON MODIFIED
-forward_run.zero-gr5-lr.jobs                              |NON MODIFIED
-forward_run.zero-gr5-lr.q_domain                          |NON MODIFIED
+forward_run.zero-gr5-lag0.rr_states.hp                    |MODIFIED
+forward_run.zero-gr5-lag0.rr_states.ht                    |MODIFIED
+forward_run.zero-gr5-lag0.sim_q                           |MODIFIED
+forward_run.zero-gr5-lr.cost                              |MODIFIED
+forward_run.zero-gr5-lr.jobs                              |MODIFIED
+forward_run.zero-gr5-lr.q_domain                          |MODIFIED
 forward_run.zero-gr5-lr.rr_states.hi                      |NON MODIFIED
 forward_run.zero-gr5-lr.rr_states.hlr                     |NON MODIFIED
-forward_run.zero-gr5-lr.rr_states.hp                      |NON MODIFIED
-forward_run.zero-gr5-lr.rr_states.ht                      |NON MODIFIED
-forward_run.zero-gr5-lr.sim_q                             |NON MODIFIED
-forward_run.zero-grd-kw.cost                              |MODIFIED
-forward_run.zero-grd-kw.jobs                              |MODIFIED
-forward_run.zero-grd-kw.q_domain                          |MODIFIED
-forward_run.zero-grd-kw.rr_states.hp                      |MODIFIED
-forward_run.zero-grd-kw.rr_states.ht                      |MODIFIED
-forward_run.zero-grd-kw.sim_q                             |MODIFIED
-forward_run.zero-grd-lag0.cost                            |MODIFIED
-forward_run.zero-grd-lag0.jobs                            |MODIFIED
-forward_run.zero-grd-lag0.q_domain                        |MODIFIED
-forward_run.zero-grd-lag0.rr_states.hp                    |MODIFIED
-forward_run.zero-grd-lag0.rr_states.ht                    |MODIFIED
-forward_run.zero-grd-lag0.sim_q                           |MODIFIED
-forward_run.zero-grd-lr.cost                              |MODIFIED
-forward_run.zero-grd-lr.jobs                              |MODIFIED
-forward_run.zero-grd-lr.q_domain                          |MODIFIED
+forward_run.zero-gr5-lr.rr_states.hp                      |MODIFIED
+forward_run.zero-gr5-lr.rr_states.ht                      |MODIFIED
+forward_run.zero-gr5-lr.sim_q                             |MODIFIED
+forward_run.zero-grd-kw.cost                              |NON MODIFIED
+forward_run.zero-grd-kw.jobs                              |NON MODIFIED
+forward_run.zero-grd-kw.q_domain                          |NON MODIFIED
+forward_run.zero-grd-kw.rr_states.hp                      |NON MODIFIED
+forward_run.zero-grd-kw.rr_states.ht                      |NON MODIFIED
+forward_run.zero-grd-kw.sim_q                             |NON MODIFIED
+forward_run.zero-grd-lag0.cost                            |NON MODIFIED
+forward_run.zero-grd-lag0.jobs                            |NON MODIFIED
+forward_run.zero-grd-lag0.q_domain                        |NON MODIFIED
+forward_run.zero-grd-lag0.rr_states.hp                    |NON MODIFIED
+forward_run.zero-grd-lag0.rr_states.ht                    |NON MODIFIED
+forward_run.zero-grd-lag0.sim_q                           |NON MODIFIED
+forward_run.zero-grd-lr.cost                              |NON MODIFIED
+forward_run.zero-grd-lr.jobs                              |NON MODIFIED
+forward_run.zero-grd-lr.q_domain                          |NON MODIFIED
 forward_run.zero-grd-lr.rr_states.hlr                     |NON MODIFIED
-forward_run.zero-grd-lr.rr_states.hp                      |MODIFIED
-forward_run.zero-grd-lr.rr_states.ht                      |MODIFIED
-forward_run.zero-grd-lr.sim_q                             |MODIFIED
-forward_run.zero-loieau-kw.cost                           |NON MODIFIED
-forward_run.zero-loieau-kw.jobs                           |NON MODIFIED
-forward_run.zero-loieau-kw.q_domain                       |NON MODIFIED
-forward_run.zero-loieau-kw.rr_states.ha                   |NON MODIFIED
-forward_run.zero-loieau-kw.rr_states.hc                   |NON MODIFIED
-forward_run.zero-loieau-kw.sim_q                          |NON MODIFIED
-forward_run.zero-loieau-lag0.cost                         |NON MODIFIED
-forward_run.zero-loieau-lag0.jobs                         |NON MODIFIED
-forward_run.zero-loieau-lag0.q_domain                     |NON MODIFIED
-forward_run.zero-loieau-lag0.rr_states.ha                 |NON MODIFIED
-forward_run.zero-loieau-lag0.rr_states.hc                 |NON MODIFIED
-forward_run.zero-loieau-lag0.sim_q                        |NON MODIFIED
-forward_run.zero-loieau-lr.cost                           |NON MODIFIED
-forward_run.zero-loieau-lr.jobs                           |NON MODIFIED
-forward_run.zero-loieau-lr.q_domain                       |NON MODIFIED
-forward_run.zero-loieau-lr.rr_states.ha                   |NON MODIFIED
-forward_run.zero-loieau-lr.rr_states.hc                   |NON MODIFIED
+forward_run.zero-grd-lr.rr_states.hp                      |NON MODIFIED
+forward_run.zero-grd-lr.rr_states.ht                      |NON MODIFIED
+forward_run.zero-grd-lr.sim_q                             |NON MODIFIED
+forward_run.zero-loieau-kw.cost                           |MODIFIED
+forward_run.zero-loieau-kw.jobs                           |MODIFIED
+forward_run.zero-loieau-kw.q_domain                       |MODIFIED
+forward_run.zero-loieau-kw.rr_states.ha                   |MODIFIED
+forward_run.zero-loieau-kw.rr_states.hc                   |MODIFIED
+forward_run.zero-loieau-kw.sim_q                          |MODIFIED
+forward_run.zero-loieau-lag0.cost                         |MODIFIED
+forward_run.zero-loieau-lag0.jobs                         |MODIFIED
+forward_run.zero-loieau-lag0.q_domain                     |MODIFIED
+forward_run.zero-loieau-lag0.rr_states.ha                 |MODIFIED
+forward_run.zero-loieau-lag0.rr_states.hc                 |MODIFIED
+forward_run.zero-loieau-lag0.sim_q                        |MODIFIED
+forward_run.zero-loieau-lr.cost                           |MODIFIED
+forward_run.zero-loieau-lr.jobs                           |MODIFIED
+forward_run.zero-loieau-lr.q_domain                       |MODIFIED
+forward_run.zero-loieau-lr.rr_states.ha                   |MODIFIED
+forward_run.zero-loieau-lr.rr_states.hc                   |MODIFIED
 forward_run.zero-loieau-lr.rr_states.hlr                  |NON MODIFIED
-forward_run.zero-loieau-lr.sim_q                          |NON MODIFIED
+forward_run.zero-loieau-lr.sim_q                          |MODIFIED
 forward_run.zero-vic3l-kw.cost                            |NON MODIFIED
 forward_run.zero-vic3l-kw.jobs                            |NON MODIFIED
 forward_run.zero-vic3l-kw.q_domain                        |NON MODIFIED
@@ -227,16 +224,16 @@ model_io.setup.end_time                                   |NON MODIFIED
 model_io.setup.start_time                                 |NON MODIFIED
 model_io.setup.structure                                  |NON MODIFIED
 multiset_estimate.mfwr.lcurve_multiset.alpha              |NON MODIFIED
-multiset_estimate.mfwr.lcurve_multiset.alpha_opt          |NON MODIFIED
-multiset_estimate.mfwr.lcurve_multiset.cost               |NON MODIFIED
-multiset_estimate.mfwr.lcurve_multiset.mahal_dist         |NON MODIFIED
-multiset_estimate.mfwr.sim_q                              |NON MODIFIED
-multiset_estimate.mopt_ml.sim_q                           |NON MODIFIED
+multiset_estimate.mfwr.lcurve_multiset.alpha_opt          |MODIFIED
+multiset_estimate.mfwr.lcurve_multiset.cost               |MODIFIED
+multiset_estimate.mfwr.lcurve_multiset.mahal_dist         |MODIFIED
+multiset_estimate.mfwr.sim_q                              |MODIFIED
+multiset_estimate.mopt_ml.sim_q                           |MODIFIED
 multiset_estimate.mopt_unf.lcurve_multiset.alpha          |NON MODIFIED
-multiset_estimate.mopt_unf.lcurve_multiset.alpha_opt      |NON MODIFIED
-multiset_estimate.mopt_unf.lcurve_multiset.cost           |NON MODIFIED
-multiset_estimate.mopt_unf.lcurve_multiset.mahal_dist     |NON MODIFIED
-multiset_estimate.mopt_unf.sim_q                          |NON MODIFIED
+multiset_estimate.mopt_unf.lcurve_multiset.alpha_opt      |MODIFIED
+multiset_estimate.mopt_unf.lcurve_multiset.cost           |MODIFIED
+multiset_estimate.mopt_unf.lcurve_multiset.mahal_dist     |MODIFIED
+multiset_estimate.mopt_unf.sim_q                          |MODIFIED
 net_init.bias_layer_1                                     |NON MODIFIED
 net_init.bias_layer_2                                     |NON MODIFIED
 net_init.bias_layer_3                                     |NON MODIFIED
@@ -246,162 +243,162 @@ net_init.weight_layer_1                                   |NON MODIFIED
 net_init.weight_layer_2                                   |NON MODIFIED
 net_init.weight_layer_3                                   |NON MODIFIED
 net_init.weight_layer_4                                   |NON MODIFIED
-optimize.zero-gr4-kw.ann.sim_q                            |NON MODIFIED
-optimize.zero-gr4-kw.distributed.control_vector           |NON MODIFIED
-optimize.zero-gr4-kw.distributed.iter_cost                |NON MODIFIED
-optimize.zero-gr4-kw.distributed.sim_q                    |NON MODIFIED
-optimize.zero-gr4-kw.multi-linear.control_vector          |NON MODIFIED
-optimize.zero-gr4-kw.multi-linear.iter_cost               |NON MODIFIED
-optimize.zero-gr4-kw.multi-linear.sim_q                   |NON MODIFIED
-optimize.zero-gr4-kw.multi-polynomial.control_vector      |NON MODIFIED
-optimize.zero-gr4-kw.multi-polynomial.iter_cost           |NON MODIFIED
-optimize.zero-gr4-kw.multi-polynomial.sim_q               |NON MODIFIED
-optimize.zero-gr4-kw.uniform.control_vector               |NON MODIFIED
-optimize.zero-gr4-kw.uniform.iter_cost                    |NON MODIFIED
-optimize.zero-gr4-kw.uniform.sim_q                        |NON MODIFIED
-optimize.zero-gr4-lag0.ann.sim_q                          |NON MODIFIED
-optimize.zero-gr4-lag0.distributed.control_vector         |NON MODIFIED
-optimize.zero-gr4-lag0.distributed.iter_cost              |NON MODIFIED
-optimize.zero-gr4-lag0.distributed.sim_q                  |NON MODIFIED
-optimize.zero-gr4-lag0.multi-linear.control_vector        |NON MODIFIED
-optimize.zero-gr4-lag0.multi-linear.iter_cost             |NON MODIFIED
-optimize.zero-gr4-lag0.multi-linear.sim_q                 |NON MODIFIED
-optimize.zero-gr4-lag0.multi-polynomial.control_vector    |NON MODIFIED
-optimize.zero-gr4-lag0.multi-polynomial.iter_cost         |NON MODIFIED
-optimize.zero-gr4-lag0.multi-polynomial.sim_q             |NON MODIFIED
-optimize.zero-gr4-lag0.uniform.control_vector             |NON MODIFIED
-optimize.zero-gr4-lag0.uniform.iter_cost                  |NON MODIFIED
-optimize.zero-gr4-lag0.uniform.sim_q                      |NON MODIFIED
-optimize.zero-gr4-lr.ann.sim_q                            |NON MODIFIED
-optimize.zero-gr4-lr.distributed.control_vector           |NON MODIFIED
-optimize.zero-gr4-lr.distributed.iter_cost                |NON MODIFIED
-optimize.zero-gr4-lr.distributed.sim_q                    |NON MODIFIED
-optimize.zero-gr4-lr.multi-linear.control_vector          |NON MODIFIED
-optimize.zero-gr4-lr.multi-linear.iter_cost               |NON MODIFIED
-optimize.zero-gr4-lr.multi-linear.sim_q                   |NON MODIFIED
-optimize.zero-gr4-lr.multi-polynomial.control_vector      |NON MODIFIED
-optimize.zero-gr4-lr.multi-polynomial.iter_cost           |NON MODIFIED
-optimize.zero-gr4-lr.multi-polynomial.sim_q               |NON MODIFIED
-optimize.zero-gr4-lr.uniform.control_vector               |NON MODIFIED
-optimize.zero-gr4-lr.uniform.iter_cost                    |NON MODIFIED
-optimize.zero-gr4-lr.uniform.sim_q                        |NON MODIFIED
-optimize.zero-gr5-kw.ann.sim_q                            |NON MODIFIED
-optimize.zero-gr5-kw.distributed.control_vector           |NON MODIFIED
-optimize.zero-gr5-kw.distributed.iter_cost                |NON MODIFIED
-optimize.zero-gr5-kw.distributed.sim_q                    |NON MODIFIED
-optimize.zero-gr5-kw.multi-linear.control_vector          |NON MODIFIED
-optimize.zero-gr5-kw.multi-linear.iter_cost               |NON MODIFIED
-optimize.zero-gr5-kw.multi-linear.sim_q                   |NON MODIFIED
-optimize.zero-gr5-kw.multi-polynomial.control_vector      |NON MODIFIED
-optimize.zero-gr5-kw.multi-polynomial.iter_cost           |NON MODIFIED
-optimize.zero-gr5-kw.multi-polynomial.sim_q               |NON MODIFIED
-optimize.zero-gr5-kw.uniform.control_vector               |NON MODIFIED
-optimize.zero-gr5-kw.uniform.iter_cost                    |NON MODIFIED
-optimize.zero-gr5-kw.uniform.sim_q                        |NON MODIFIED
-optimize.zero-gr5-lag0.ann.sim_q                          |NON MODIFIED
-optimize.zero-gr5-lag0.distributed.control_vector         |NON MODIFIED
-optimize.zero-gr5-lag0.distributed.iter_cost              |NON MODIFIED
-optimize.zero-gr5-lag0.distributed.sim_q                  |NON MODIFIED
-optimize.zero-gr5-lag0.multi-linear.control_vector        |NON MODIFIED
-optimize.zero-gr5-lag0.multi-linear.iter_cost             |NON MODIFIED
-optimize.zero-gr5-lag0.multi-linear.sim_q                 |NON MODIFIED
-optimize.zero-gr5-lag0.multi-polynomial.control_vector    |NON MODIFIED
-optimize.zero-gr5-lag0.multi-polynomial.iter_cost         |NON MODIFIED
-optimize.zero-gr5-lag0.multi-polynomial.sim_q             |NON MODIFIED
-optimize.zero-gr5-lag0.uniform.control_vector             |NON MODIFIED
-optimize.zero-gr5-lag0.uniform.iter_cost                  |NON MODIFIED
-optimize.zero-gr5-lag0.uniform.sim_q                      |NON MODIFIED
-optimize.zero-gr5-lr.ann.sim_q                            |NON MODIFIED
-optimize.zero-gr5-lr.distributed.control_vector           |NON MODIFIED
-optimize.zero-gr5-lr.distributed.iter_cost                |NON MODIFIED
-optimize.zero-gr5-lr.distributed.sim_q                    |NON MODIFIED
-optimize.zero-gr5-lr.multi-linear.control_vector          |NON MODIFIED
-optimize.zero-gr5-lr.multi-linear.iter_cost               |NON MODIFIED
-optimize.zero-gr5-lr.multi-linear.sim_q                   |NON MODIFIED
-optimize.zero-gr5-lr.multi-polynomial.control_vector      |NON MODIFIED
-optimize.zero-gr5-lr.multi-polynomial.iter_cost           |NON MODIFIED
-optimize.zero-gr5-lr.multi-polynomial.sim_q               |NON MODIFIED
-optimize.zero-gr5-lr.uniform.control_vector               |NON MODIFIED
-optimize.zero-gr5-lr.uniform.iter_cost                    |NON MODIFIED
-optimize.zero-gr5-lr.uniform.sim_q                        |NON MODIFIED
-optimize.zero-grd-kw.ann.sim_q                            |MODIFIED
-optimize.zero-grd-kw.distributed.control_vector           |MODIFIED
-optimize.zero-grd-kw.distributed.iter_cost                |MODIFIED
-optimize.zero-grd-kw.distributed.sim_q                    |MODIFIED
-optimize.zero-grd-kw.multi-linear.control_vector          |MODIFIED
-optimize.zero-grd-kw.multi-linear.iter_cost               |MODIFIED
-optimize.zero-grd-kw.multi-linear.sim_q                   |MODIFIED
-optimize.zero-grd-kw.multi-polynomial.control_vector      |MODIFIED
-optimize.zero-grd-kw.multi-polynomial.iter_cost           |MODIFIED
-optimize.zero-grd-kw.multi-polynomial.sim_q               |MODIFIED
-optimize.zero-grd-kw.uniform.control_vector               |MODIFIED
-optimize.zero-grd-kw.uniform.iter_cost                    |MODIFIED
-optimize.zero-grd-kw.uniform.sim_q                        |MODIFIED
-optimize.zero-grd-lag0.ann.sim_q                          |MODIFIED
-optimize.zero-grd-lag0.distributed.control_vector         |MODIFIED
-optimize.zero-grd-lag0.distributed.iter_cost              |MODIFIED
-optimize.zero-grd-lag0.distributed.sim_q                  |MODIFIED
-optimize.zero-grd-lag0.multi-linear.control_vector        |MODIFIED
-optimize.zero-grd-lag0.multi-linear.iter_cost             |MODIFIED
-optimize.zero-grd-lag0.multi-linear.sim_q                 |MODIFIED
-optimize.zero-grd-lag0.multi-polynomial.control_vector    |MODIFIED
-optimize.zero-grd-lag0.multi-polynomial.iter_cost         |MODIFIED
-optimize.zero-grd-lag0.multi-polynomial.sim_q             |MODIFIED
-optimize.zero-grd-lag0.uniform.control_vector             |MODIFIED
-optimize.zero-grd-lag0.uniform.iter_cost                  |MODIFIED
-optimize.zero-grd-lag0.uniform.sim_q                      |MODIFIED
-optimize.zero-grd-lr.ann.sim_q                            |MODIFIED
-optimize.zero-grd-lr.distributed.control_vector           |MODIFIED
-optimize.zero-grd-lr.distributed.iter_cost                |MODIFIED
-optimize.zero-grd-lr.distributed.sim_q                    |MODIFIED
-optimize.zero-grd-lr.multi-linear.control_vector          |MODIFIED
-optimize.zero-grd-lr.multi-linear.iter_cost               |MODIFIED
-optimize.zero-grd-lr.multi-linear.sim_q                   |MODIFIED
-optimize.zero-grd-lr.multi-polynomial.control_vector      |MODIFIED
-optimize.zero-grd-lr.multi-polynomial.iter_cost           |MODIFIED
-optimize.zero-grd-lr.multi-polynomial.sim_q               |MODIFIED
-optimize.zero-grd-lr.uniform.control_vector               |MODIFIED
-optimize.zero-grd-lr.uniform.iter_cost                    |MODIFIED
-optimize.zero-grd-lr.uniform.sim_q                        |MODIFIED
-optimize.zero-loieau-kw.ann.sim_q                         |NON MODIFIED
-optimize.zero-loieau-kw.distributed.control_vector        |NON MODIFIED
-optimize.zero-loieau-kw.distributed.iter_cost             |NON MODIFIED
-optimize.zero-loieau-kw.distributed.sim_q                 |NON MODIFIED
-optimize.zero-loieau-kw.multi-linear.control_vector       |NON MODIFIED
-optimize.zero-loieau-kw.multi-linear.iter_cost            |NON MODIFIED
-optimize.zero-loieau-kw.multi-linear.sim_q                |NON MODIFIED
-optimize.zero-loieau-kw.multi-polynomial.control_vector   |NON MODIFIED
-optimize.zero-loieau-kw.multi-polynomial.iter_cost        |NON MODIFIED
-optimize.zero-loieau-kw.multi-polynomial.sim_q            |NON MODIFIED
-optimize.zero-loieau-kw.uniform.control_vector            |NON MODIFIED
-optimize.zero-loieau-kw.uniform.iter_cost                 |NON MODIFIED
-optimize.zero-loieau-kw.uniform.sim_q                     |NON MODIFIED
-optimize.zero-loieau-lag0.ann.sim_q                       |NON MODIFIED
-optimize.zero-loieau-lag0.distributed.control_vector      |NON MODIFIED
-optimize.zero-loieau-lag0.distributed.iter_cost           |NON MODIFIED
-optimize.zero-loieau-lag0.distributed.sim_q               |NON MODIFIED
-optimize.zero-loieau-lag0.multi-linear.control_vector     |NON MODIFIED
-optimize.zero-loieau-lag0.multi-linear.iter_cost          |NON MODIFIED
-optimize.zero-loieau-lag0.multi-linear.sim_q              |NON MODIFIED
-optimize.zero-loieau-lag0.multi-polynomial.control_vector |NON MODIFIED
-optimize.zero-loieau-lag0.multi-polynomial.iter_cost      |NON MODIFIED
-optimize.zero-loieau-lag0.multi-polynomial.sim_q          |NON MODIFIED
-optimize.zero-loieau-lag0.uniform.control_vector          |NON MODIFIED
-optimize.zero-loieau-lag0.uniform.iter_cost               |NON MODIFIED
-optimize.zero-loieau-lag0.uniform.sim_q                   |NON MODIFIED
-optimize.zero-loieau-lr.ann.sim_q                         |NON MODIFIED
-optimize.zero-loieau-lr.distributed.control_vector        |NON MODIFIED
-optimize.zero-loieau-lr.distributed.iter_cost             |NON MODIFIED
-optimize.zero-loieau-lr.distributed.sim_q                 |NON MODIFIED
-optimize.zero-loieau-lr.multi-linear.control_vector       |NON MODIFIED
-optimize.zero-loieau-lr.multi-linear.iter_cost            |NON MODIFIED
-optimize.zero-loieau-lr.multi-linear.sim_q                |NON MODIFIED
-optimize.zero-loieau-lr.multi-polynomial.control_vector   |NON MODIFIED
-optimize.zero-loieau-lr.multi-polynomial.iter_cost        |NON MODIFIED
-optimize.zero-loieau-lr.multi-polynomial.sim_q            |NON MODIFIED
-optimize.zero-loieau-lr.uniform.control_vector            |NON MODIFIED
-optimize.zero-loieau-lr.uniform.iter_cost                 |NON MODIFIED
-optimize.zero-loieau-lr.uniform.sim_q                     |NON MODIFIED
+optimize.zero-gr4-kw.ann.sim_q                            |MODIFIED
+optimize.zero-gr4-kw.distributed.control_vector           |MODIFIED
+optimize.zero-gr4-kw.distributed.iter_cost                |MODIFIED
+optimize.zero-gr4-kw.distributed.sim_q                    |MODIFIED
+optimize.zero-gr4-kw.multi-linear.control_vector          |MODIFIED
+optimize.zero-gr4-kw.multi-linear.iter_cost               |MODIFIED
+optimize.zero-gr4-kw.multi-linear.sim_q                   |MODIFIED
+optimize.zero-gr4-kw.multi-polynomial.control_vector      |MODIFIED
+optimize.zero-gr4-kw.multi-polynomial.iter_cost           |MODIFIED
+optimize.zero-gr4-kw.multi-polynomial.sim_q               |MODIFIED
+optimize.zero-gr4-kw.uniform.control_vector               |MODIFIED
+optimize.zero-gr4-kw.uniform.iter_cost                    |MODIFIED
+optimize.zero-gr4-kw.uniform.sim_q                        |MODIFIED
+optimize.zero-gr4-lag0.ann.sim_q                          |MODIFIED
+optimize.zero-gr4-lag0.distributed.control_vector         |MODIFIED
+optimize.zero-gr4-lag0.distributed.iter_cost              |MODIFIED
+optimize.zero-gr4-lag0.distributed.sim_q                  |MODIFIED
+optimize.zero-gr4-lag0.multi-linear.control_vector        |MODIFIED
+optimize.zero-gr4-lag0.multi-linear.iter_cost             |MODIFIED
+optimize.zero-gr4-lag0.multi-linear.sim_q                 |MODIFIED
+optimize.zero-gr4-lag0.multi-polynomial.control_vector    |MODIFIED
+optimize.zero-gr4-lag0.multi-polynomial.iter_cost         |MODIFIED
+optimize.zero-gr4-lag0.multi-polynomial.sim_q             |MODIFIED
+optimize.zero-gr4-lag0.uniform.control_vector             |MODIFIED
+optimize.zero-gr4-lag0.uniform.iter_cost                  |MODIFIED
+optimize.zero-gr4-lag0.uniform.sim_q                      |MODIFIED
+optimize.zero-gr4-lr.ann.sim_q                            |MODIFIED
+optimize.zero-gr4-lr.distributed.control_vector           |MODIFIED
+optimize.zero-gr4-lr.distributed.iter_cost                |MODIFIED
+optimize.zero-gr4-lr.distributed.sim_q                    |MODIFIED
+optimize.zero-gr4-lr.multi-linear.control_vector          |MODIFIED
+optimize.zero-gr4-lr.multi-linear.iter_cost               |MODIFIED
+optimize.zero-gr4-lr.multi-linear.sim_q                   |MODIFIED
+optimize.zero-gr4-lr.multi-polynomial.control_vector      |MODIFIED
+optimize.zero-gr4-lr.multi-polynomial.iter_cost           |MODIFIED
+optimize.zero-gr4-lr.multi-polynomial.sim_q               |MODIFIED
+optimize.zero-gr4-lr.uniform.control_vector               |MODIFIED
+optimize.zero-gr4-lr.uniform.iter_cost                    |MODIFIED
+optimize.zero-gr4-lr.uniform.sim_q                        |MODIFIED
+optimize.zero-gr5-kw.ann.sim_q                            |MODIFIED
+optimize.zero-gr5-kw.distributed.control_vector           |MODIFIED
+optimize.zero-gr5-kw.distributed.iter_cost                |MODIFIED
+optimize.zero-gr5-kw.distributed.sim_q                    |MODIFIED
+optimize.zero-gr5-kw.multi-linear.control_vector          |MODIFIED
+optimize.zero-gr5-kw.multi-linear.iter_cost               |MODIFIED
+optimize.zero-gr5-kw.multi-linear.sim_q                   |MODIFIED
+optimize.zero-gr5-kw.multi-polynomial.control_vector      |MODIFIED
+optimize.zero-gr5-kw.multi-polynomial.iter_cost           |MODIFIED
+optimize.zero-gr5-kw.multi-polynomial.sim_q               |MODIFIED
+optimize.zero-gr5-kw.uniform.control_vector               |MODIFIED
+optimize.zero-gr5-kw.uniform.iter_cost                    |MODIFIED
+optimize.zero-gr5-kw.uniform.sim_q                        |MODIFIED
+optimize.zero-gr5-lag0.ann.sim_q                          |MODIFIED
+optimize.zero-gr5-lag0.distributed.control_vector         |MODIFIED
+optimize.zero-gr5-lag0.distributed.iter_cost              |MODIFIED
+optimize.zero-gr5-lag0.distributed.sim_q                  |MODIFIED
+optimize.zero-gr5-lag0.multi-linear.control_vector        |MODIFIED
+optimize.zero-gr5-lag0.multi-linear.iter_cost             |MODIFIED
+optimize.zero-gr5-lag0.multi-linear.sim_q                 |MODIFIED
+optimize.zero-gr5-lag0.multi-polynomial.control_vector    |MODIFIED
+optimize.zero-gr5-lag0.multi-polynomial.iter_cost         |MODIFIED
+optimize.zero-gr5-lag0.multi-polynomial.sim_q             |MODIFIED
+optimize.zero-gr5-lag0.uniform.control_vector             |MODIFIED
+optimize.zero-gr5-lag0.uniform.iter_cost                  |MODIFIED
+optimize.zero-gr5-lag0.uniform.sim_q                      |MODIFIED
+optimize.zero-gr5-lr.ann.sim_q                            |MODIFIED
+optimize.zero-gr5-lr.distributed.control_vector           |MODIFIED
+optimize.zero-gr5-lr.distributed.iter_cost                |MODIFIED
+optimize.zero-gr5-lr.distributed.sim_q                    |MODIFIED
+optimize.zero-gr5-lr.multi-linear.control_vector          |MODIFIED
+optimize.zero-gr5-lr.multi-linear.iter_cost               |MODIFIED
+optimize.zero-gr5-lr.multi-linear.sim_q                   |MODIFIED
+optimize.zero-gr5-lr.multi-polynomial.control_vector      |MODIFIED
+optimize.zero-gr5-lr.multi-polynomial.iter_cost           |MODIFIED
+optimize.zero-gr5-lr.multi-polynomial.sim_q               |MODIFIED
+optimize.zero-gr5-lr.uniform.control_vector               |MODIFIED
+optimize.zero-gr5-lr.uniform.iter_cost                    |MODIFIED
+optimize.zero-gr5-lr.uniform.sim_q                        |MODIFIED
+optimize.zero-grd-kw.ann.sim_q                            |NON MODIFIED
+optimize.zero-grd-kw.distributed.control_vector           |NON MODIFIED
+optimize.zero-grd-kw.distributed.iter_cost                |NON MODIFIED
+optimize.zero-grd-kw.distributed.sim_q                    |NON MODIFIED
+optimize.zero-grd-kw.multi-linear.control_vector          |NON MODIFIED
+optimize.zero-grd-kw.multi-linear.iter_cost               |NON MODIFIED
+optimize.zero-grd-kw.multi-linear.sim_q                   |NON MODIFIED
+optimize.zero-grd-kw.multi-polynomial.control_vector      |NON MODIFIED
+optimize.zero-grd-kw.multi-polynomial.iter_cost           |NON MODIFIED
+optimize.zero-grd-kw.multi-polynomial.sim_q               |NON MODIFIED
+optimize.zero-grd-kw.uniform.control_vector               |NON MODIFIED
+optimize.zero-grd-kw.uniform.iter_cost                    |NON MODIFIED
+optimize.zero-grd-kw.uniform.sim_q                        |NON MODIFIED
+optimize.zero-grd-lag0.ann.sim_q                          |NON MODIFIED
+optimize.zero-grd-lag0.distributed.control_vector         |NON MODIFIED
+optimize.zero-grd-lag0.distributed.iter_cost              |NON MODIFIED
+optimize.zero-grd-lag0.distributed.sim_q                  |NON MODIFIED
+optimize.zero-grd-lag0.multi-linear.control_vector        |NON MODIFIED
+optimize.zero-grd-lag0.multi-linear.iter_cost             |NON MODIFIED
+optimize.zero-grd-lag0.multi-linear.sim_q                 |NON MODIFIED
+optimize.zero-grd-lag0.multi-polynomial.control_vector    |NON MODIFIED
+optimize.zero-grd-lag0.multi-polynomial.iter_cost         |NON MODIFIED
+optimize.zero-grd-lag0.multi-polynomial.sim_q             |NON MODIFIED
+optimize.zero-grd-lag0.uniform.control_vector             |NON MODIFIED
+optimize.zero-grd-lag0.uniform.iter_cost                  |NON MODIFIED
+optimize.zero-grd-lag0.uniform.sim_q                      |NON MODIFIED
+optimize.zero-grd-lr.ann.sim_q                            |NON MODIFIED
+optimize.zero-grd-lr.distributed.control_vector           |NON MODIFIED
+optimize.zero-grd-lr.distributed.iter_cost                |NON MODIFIED
+optimize.zero-grd-lr.distributed.sim_q                    |NON MODIFIED
+optimize.zero-grd-lr.multi-linear.control_vector          |NON MODIFIED
+optimize.zero-grd-lr.multi-linear.iter_cost               |NON MODIFIED
+optimize.zero-grd-lr.multi-linear.sim_q                   |NON MODIFIED
+optimize.zero-grd-lr.multi-polynomial.control_vector      |NON MODIFIED
+optimize.zero-grd-lr.multi-polynomial.iter_cost           |NON MODIFIED
+optimize.zero-grd-lr.multi-polynomial.sim_q               |NON MODIFIED
+optimize.zero-grd-lr.uniform.control_vector               |NON MODIFIED
+optimize.zero-grd-lr.uniform.iter_cost                    |NON MODIFIED
+optimize.zero-grd-lr.uniform.sim_q                        |NON MODIFIED
+optimize.zero-loieau-kw.ann.sim_q                         |MODIFIED
+optimize.zero-loieau-kw.distributed.control_vector        |MODIFIED
+optimize.zero-loieau-kw.distributed.iter_cost             |MODIFIED
+optimize.zero-loieau-kw.distributed.sim_q                 |MODIFIED
+optimize.zero-loieau-kw.multi-linear.control_vector       |MODIFIED
+optimize.zero-loieau-kw.multi-linear.iter_cost            |MODIFIED
+optimize.zero-loieau-kw.multi-linear.sim_q                |MODIFIED
+optimize.zero-loieau-kw.multi-polynomial.control_vector   |MODIFIED
+optimize.zero-loieau-kw.multi-polynomial.iter_cost        |MODIFIED
+optimize.zero-loieau-kw.multi-polynomial.sim_q            |MODIFIED
+optimize.zero-loieau-kw.uniform.control_vector            |MODIFIED
+optimize.zero-loieau-kw.uniform.iter_cost                 |MODIFIED
+optimize.zero-loieau-kw.uniform.sim_q                     |MODIFIED
+optimize.zero-loieau-lag0.ann.sim_q                       |MODIFIED
+optimize.zero-loieau-lag0.distributed.control_vector      |MODIFIED
+optimize.zero-loieau-lag0.distributed.iter_cost           |MODIFIED
+optimize.zero-loieau-lag0.distributed.sim_q               |MODIFIED
+optimize.zero-loieau-lag0.multi-linear.control_vector     |MODIFIED
+optimize.zero-loieau-lag0.multi-linear.iter_cost          |MODIFIED
+optimize.zero-loieau-lag0.multi-linear.sim_q              |MODIFIED
+optimize.zero-loieau-lag0.multi-polynomial.control_vector |MODIFIED
+optimize.zero-loieau-lag0.multi-polynomial.iter_cost      |MODIFIED
+optimize.zero-loieau-lag0.multi-polynomial.sim_q          |MODIFIED
+optimize.zero-loieau-lag0.uniform.control_vector          |MODIFIED
+optimize.zero-loieau-lag0.uniform.iter_cost               |MODIFIED
+optimize.zero-loieau-lag0.uniform.sim_q                   |MODIFIED
+optimize.zero-loieau-lr.ann.sim_q                         |MODIFIED
+optimize.zero-loieau-lr.distributed.control_vector        |MODIFIED
+optimize.zero-loieau-lr.distributed.iter_cost             |MODIFIED
+optimize.zero-loieau-lr.distributed.sim_q                 |MODIFIED
+optimize.zero-loieau-lr.multi-linear.control_vector       |MODIFIED
+optimize.zero-loieau-lr.multi-linear.iter_cost            |MODIFIED
+optimize.zero-loieau-lr.multi-linear.sim_q                |MODIFIED
+optimize.zero-loieau-lr.multi-polynomial.control_vector   |MODIFIED
+optimize.zero-loieau-lr.multi-polynomial.iter_cost        |MODIFIED
+optimize.zero-loieau-lr.multi-polynomial.sim_q            |MODIFIED
+optimize.zero-loieau-lr.uniform.control_vector            |MODIFIED
+optimize.zero-loieau-lr.uniform.iter_cost                 |MODIFIED
+optimize.zero-loieau-lr.uniform.sim_q                     |MODIFIED
 optimize.zero-vic3l-kw.ann.sim_q                          |NON MODIFIED
 optimize.zero-vic3l-kw.distributed.control_vector         |NON MODIFIED
 optimize.zero-vic3l-kw.distributed.iter_cost              |NON MODIFIED