diff --git a/deps/evalhyd b/deps/evalhyd
index b20321058a877c029040ef5bad0d1298b2cd05b3..cbe1588dfbdfbfb3a2c66d49b5187955b57bc59f 160000
--- a/deps/evalhyd
+++ b/deps/evalhyd
@@ -1 +1 @@
-Subproject commit b20321058a877c029040ef5bad0d1298b2cd05b3
+Subproject commit cbe1588dfbdfbfb3a2c66d49b5187955b57bc59f
diff --git a/evalhyd/evald.py b/evalhyd/evald.py
index 05ba76ec8f6582b735dcd811b045765b6138e10e..80d56aab777b4eba5e765b801640461f16ad73a2 100644
--- a/evalhyd/evald.py
+++ b/evalhyd/evald.py
@@ -11,6 +11,8 @@ except ImportError:
 def evald(q_obs: NDArray[dtype('float64')],
           q_prd: NDArray[dtype('float64')],
           metrics: List[str],
+          q_thr: NDArray[dtype('float64')] = None,
+          events: str = None,
           transform: str = None,
           exponent: float = None,
           epsilon: float = None,
@@ -18,7 +20,8 @@ def evald(q_obs: NDArray[dtype('float64')],
           m_cdt: NDArray[dtype('|S32')] = None,
           bootstrap: Dict[str, int] = None,
           dts: NDArray[dtype('|S32')] = None,
-          seed: int = None) -> List[NDArray[dtype('float64')]]:
+          seed: int = None,
+          diagnostics: List[str] = None) -> List[NDArray[dtype('float64')]]:
     """Function to evaluate deterministic streamflow predictions"""
 
     # required arguments
@@ -30,6 +33,12 @@ def evald(q_obs: NDArray[dtype('float64')],
     }
 
     # optional arguments
+    if q_thr is not None:
+        kwargs['q_thr'] = (
+            q_thr.reshape(1, q_thr.size) if q_thr.ndim == 1 else q_thr
+        )
+    if events is not None:
+        kwargs['events'] = events
     if transform is not None:
         kwargs['transform'] = transform
     if exponent is not None:
@@ -46,11 +55,14 @@ def evald(q_obs: NDArray[dtype('float64')],
         kwargs['dts'] = dts
     if seed is not None:
         kwargs['seed'] = seed
+    if diagnostics is not None:
+        kwargs['diagnostics'] = diagnostics
 
     # check array ranks
     _expected = {
         'q_obs': 2,
         'q_prd': 2,
+        'q_thr': 2,
         't_msk': 3,
         'm_cdt': 2,
         'dts': 1
diff --git a/evalhyd/evalp.py b/evalhyd/evalp.py
index 17e2c5955deafa75e4a85409db7a95655e688145..4639bea84656b00b49dd3766a53443a058a5443c 100644
--- a/evalhyd/evalp.py
+++ b/evalhyd/evalp.py
@@ -18,7 +18,8 @@ def evalp(q_obs: NDArray[dtype('float64')],
           m_cdt: NDArray[dtype('|S32')] = None,
           bootstrap: Dict[str, int] = None,
           dts: NDArray[dtype('|S32')] = None,
-          seed: int = None) -> List[NDArray[dtype('float64')]]:
+          seed: int = None,
+          diagnostics: List[str] = None) -> List[NDArray[dtype('float64')]]:
     """Function to evaluate probabilistic streamflow predictions"""
 
     # required arguments
@@ -45,6 +46,8 @@ def evalp(q_obs: NDArray[dtype('float64')],
         kwargs['dts'] = dts
     if seed is not None:
         kwargs['seed'] = seed
+    if diagnostics is not None:
+        kwargs['diagnostics'] = diagnostics
 
     # check array ranks
     _expected = {
diff --git a/evalhyd/src/evalhyd.cpp b/evalhyd/src/evalhyd.cpp
index 720c682b3bf7a5f0c37ca30a7b1d73b40799651f..47c424d932c41b4dced0ad4bf248039701b64bad 100644
--- a/evalhyd/src/evalhyd.cpp
+++ b/evalhyd/src/evalhyd.cpp
@@ -21,6 +21,8 @@ auto evald(
     const xt::pytensor<double, 2>& q_obs,
     const xt::pytensor<double, 2>& q_prd,
     const std::vector<std::string>& metrics,
+    const xt::pytensor<double, 2>& q_thr,
+    std::optional<std::string> events,
     std::optional<std::string> transform,
     std::optional<double> exponent,
     std::optional<double> epsilon,
@@ -28,13 +30,16 @@ auto evald(
     const xt::pytensor<std::array<char, 32>, 2>& m_cdt,
     std::optional<std::unordered_map<std::string, int>> bootstrap,
     const std::vector<std::string>& dts,
-    std::optional<int> seed
+    std::optional<int> seed,
+    std::optional<std::vector<std::string>> diagnostics
 )
 {
     return evalhyd::evald(
         q_obs,
         q_prd,
         metrics,
+        q_thr,
+        (events.has_value()) ? events.value() : xtl::missing<std::string>(),
         (transform.has_value()) ? transform.value() : xtl::missing<std::string>(),
         (exponent.has_value()) ? exponent.value() : xtl::missing<double>(),
         (epsilon.has_value()) ? epsilon.value() : xtl::missing<double>(),
@@ -44,7 +49,10 @@ auto evald(
         ? bootstrap.value()
         : xtl::missing<std::unordered_map<std::string, int>>(),
         dts,
-        (seed.has_value()) ? seed.value() : xtl::missing<int>()
+        (seed.has_value()) ? seed.value() : xtl::missing<int>(),
+        (diagnostics.has_value())
+        ? diagnostics.value()
+        : xtl::missing<std::vector<std::string>>()
     );
 }
 
@@ -59,7 +67,8 @@ auto evalp(
     const xt::pytensor<std::array<char, 32>, 2>& m_cdt,
     std::optional<std::unordered_map<std::string, int>> bootstrap,
     const std::vector<std::string>& dts,
-    std::optional<int> seed
+    std::optional<int> seed,
+    std::optional<std::vector<std::string>> diagnostics
 )
 {
     return evalhyd::evalp(
@@ -75,7 +84,10 @@ auto evalp(
         ? bootstrap.value()
         : xtl::missing<std::unordered_map<std::string, int>>(),
         dts,
-        (seed.has_value()) ? seed.value() : xtl::missing<int>()
+        (seed.has_value()) ? seed.value() : xtl::missing<int>(),
+        (diagnostics.has_value())
+        ? diagnostics.value()
+        : xtl::missing<std::vector<std::string>>()
     );
 }
 
@@ -94,6 +106,8 @@ PYBIND11_MODULE(_evalhyd, m)
         py::arg("q_obs"),
         py::arg("q_prd"),
         py::arg("metrics"),
+        py::arg("q_thr") = xt::pytensor<double, 2>({0}),
+        py::arg("events") = py::none(),
         py::arg("transform") = py::none(),
         py::arg("exponent") = py::none(),
         py::arg("epsilon") = py::none(),
@@ -101,7 +115,8 @@ PYBIND11_MODULE(_evalhyd, m)
         py::arg("m_cdt") = xt::pytensor<std::array<char, 32>, 2>({0}),
         py::arg("bootstrap") = py::none(),
         py::arg("dts") = py::list(),
-        py::arg("seed") = py::none()
+        py::arg("seed") = py::none(),
+        py::arg("diagnostics") = py::none()
     );
 
     // probabilistic evaluation
@@ -119,7 +134,8 @@ PYBIND11_MODULE(_evalhyd, m)
         py::arg("m_cdt") = xt::pytensor<std::array<char, 32>, 2>({0}),
         py::arg("bootstrap") = py::none(),
         py::arg("dts") = py::list(),
-        py::arg("seed") = py::none()
+        py::arg("seed") = py::none(),
+        py::arg("diagnostics") = py::none()
     );
 
 #ifdef VERSION_INFO
diff --git a/tests/test_determinist.py b/tests/test_determinist.py
index a2750e9d529de955970fa9ecd793e2ef108e3338..8ff4605f8edc246f1af932ccfdb1116c50c0d9c1 100644
--- a/tests/test_determinist.py
+++ b/tests/test_determinist.py
@@ -47,28 +47,28 @@ class TestTransform(unittest.TestCase):
 
     def test_transform_sqrt(self):
         numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], "sqrt")[0],
+            evalhyd.evald(_obs, _prd, ["NSE"], transform="sqrt")[0],
             evalhyd.evald(_obs ** 0.5, _prd ** 0.5, ["NSE"])[0]
         )
 
     def test_transform_inv(self):
         eps = 0.01 * numpy.mean(_obs)
         numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], "inv")[0],
+            evalhyd.evald(_obs, _prd, ["NSE"], transform="inv")[0],
             evalhyd.evald(1 / (_obs + eps), 1 / (_prd + eps), ["NSE"])[0]
         )
 
     def test_transform_log(self):
         eps = 0.01 * numpy.mean(_obs)
         numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], "log")[0],
+            evalhyd.evald(_obs, _prd, ["NSE"], transform="log")[0],
             evalhyd.evald(numpy.log(_obs + eps), numpy.log(_prd + eps),
                           ["NSE"])[0]
         )
 
     def test_transform_pow(self):
         numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], "pow", exponent=0.3)[0],
+            evalhyd.evald(_obs, _prd, ["NSE"], transform="pow", exponent=0.3)[0],
             evalhyd.evald(_obs ** 0.3, _prd ** 0.3, ["NSE"])[0]
         )