diff --git a/tests/expected/evald/CONT_TBL.csv b/tests/expected/evald/CONT_TBL.csv
new file mode 100644
index 0000000000000000000000000000000000000000..2a6370d276de47340b13af813dc3d952d6139c34
--- /dev/null
+++ b/tests/expected/evald/CONT_TBL.csv
@@ -0,0 +1,204 @@
+157.,19.,15.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+157.,19.,15.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+157.,19.,15.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+157.,19.,15.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,19.,14.,120.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,21.,14.,118.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,21.,14.,118.
+200.,15.,8.,88.
+220.,21.,6.,64.
+nan,nan,nan,nan
+158.,21.,14.,118.
+200.,16.,8.,87.
+221.,21.,5.,64.
+nan,nan,nan,nan
diff --git a/tests/test_determinist.py b/tests/test_determinist.py
index 980493a4c076b5a1336819af136966fe6c5f494f..7757d64602e11f6fa189033107695102b9de9041 100644
--- a/tests/test_determinist.py
+++ b/tests/test_determinist.py
@@ -8,12 +8,19 @@ import evalhyd
 _prd = numpy.genfromtxt("./data/q_prd.csv", delimiter=',')[:, :]
 _obs = numpy.genfromtxt("./data/q_obs.csv", delimiter=',')[numpy.newaxis, :]
 
+_thr = numpy.repeat(
+    numpy.array([[690, 534, 445, numpy.nan]]), repeats=_prd.shape[0], axis=0
+)
+_events = "high"
+
 # list all available deterministic metrics
 _all_metrics = (
     # errors-based
     'MAE', 'MARE', 'MSE', 'RMSE',
     # efficiencies-based
     'NSE', 'KGE', 'KGE_D', 'KGEPRIME', 'KGEPRIME_D',
+    # contingency table-based
+    'CONT_TBL'
 )
 
 # list all available deterministic diagnostics
@@ -30,12 +37,18 @@ class TestMetrics(unittest.TestCase):
             [:, numpy.newaxis, numpy.newaxis]
         ) for metric in _all_metrics
     }
+    # /!\ stacked-up thresholds in CSV file for CONT_TBL
+    #     because 5D metric so need to reshape array
+    expected['CONT_TBL'] = (
+        expected['CONT_TBL'].reshape(expected['NSE'].shape + (4, 4))
+    )
 
     def test_metrics_2d(self):
         for metric in self.expected.keys():
             with self.subTest(metric=metric):
                 numpy.testing.assert_almost_equal(
-                    evalhyd.evald(_obs, _prd, [metric])[0],
+                    evalhyd.evald(_obs, _prd, [metric],
+                                  q_thr=_thr, events=_events)[0],
                     self.expected[metric]
                 )
 
@@ -43,7 +56,8 @@ class TestMetrics(unittest.TestCase):
         for metric in self.expected.keys():
             with self.subTest(metric=metric):
                 numpy.testing.assert_almost_equal(
-                    evalhyd.evald(_obs[0], _prd[0], [metric])[0],
+                    evalhyd.evald(_obs[0], _prd[0], [metric],
+                                  q_thr=_thr[0], events=_events)[0],
                     [self.expected[metric][0]]
                 )
 
@@ -51,31 +65,48 @@ class TestMetrics(unittest.TestCase):
 class TestTransform(unittest.TestCase):
 
     def test_transform_sqrt(self):
-        numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], transform="sqrt")[0],
-            evalhyd.evald(_obs ** 0.5, _prd ** 0.5, ["NSE"])[0]
-        )
+        for metric in _all_metrics:
+            with self.subTest(metric=metric):
+                numpy.testing.assert_almost_equal(
+                    evalhyd.evald(_obs, _prd, ["NSE"], transform="sqrt",
+                                  q_thr=_thr, events=_events)[0],
+                    evalhyd.evald(_obs ** 0.5, _prd ** 0.5, ["NSE"],
+                                  q_thr=_thr ** 0.5, events=_events)[0]
+                )
 
     def test_transform_inv(self):
         eps = 0.01 * numpy.mean(_obs)
-        numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], transform="inv")[0],
-            evalhyd.evald(1 / (_obs + eps), 1 / (_prd + eps), ["NSE"])[0]
-        )
+        for metric in _all_metrics:
+            with self.subTest(metric=metric):
+                numpy.testing.assert_almost_equal(
+                    evalhyd.evald(_obs, _prd, ["NSE"], transform="inv",
+                                  q_thr=_thr, events=_events)[0],
+                    evalhyd.evald(1 / (_obs + eps), 1 / (_prd + eps), ["NSE"],
+                                  q_thr=1 / (_thr + eps), events=_events)[0]
+                )
 
     def test_transform_log(self):
         eps = 0.01 * numpy.mean(_obs)
-        numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], transform="log")[0],
-            evalhyd.evald(numpy.log(_obs + eps), numpy.log(_prd + eps),
-                          ["NSE"])[0]
-        )
+        for metric in _all_metrics:
+            with self.subTest(metric=metric):
+                numpy.testing.assert_almost_equal(
+                    evalhyd.evald(_obs, _prd, ["NSE"], transform="log",
+                                  q_thr=_thr, events=_events)[0],
+                    evalhyd.evald(numpy.log(_obs + eps), numpy.log(_prd + eps),
+                                  ["NSE"],
+                                  q_thr=numpy.log(_thr + eps), events=_events)[0]
+                )
 
     def test_transform_pow(self):
-        numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], transform="pow", exponent=0.3)[0],
-            evalhyd.evald(_obs ** 0.3, _prd ** 0.3, ["NSE"])[0]
-        )
+        for metric in _all_metrics:
+            with self.subTest(metric=metric):
+                numpy.testing.assert_almost_equal(
+                    evalhyd.evald(_obs, _prd, ["NSE"],
+                                  q_thr=_thr, events=_events,
+                                  transform="pow", exponent=0.3)[0],
+                    evalhyd.evald(_obs ** 0.3, _prd ** 0.3, ["NSE"],
+                                  q_thr=_thr ** 0.3, events=_events)[0]
+                )
 
 
 class TestMasking(unittest.TestCase):
@@ -90,8 +121,10 @@ class TestMasking(unittest.TestCase):
         prd = _prd[..., 99:].copy()
 
         numpy.testing.assert_almost_equal(
-            evalhyd.evald(_obs, _prd, ["NSE"], t_msk=msk)[0],
-            evalhyd.evald(obs, prd, ["NSE"])[0]
+            evalhyd.evald(_obs, _prd, ["NSE"],
+                          q_thr=_thr, events=_events, t_msk=msk)[0],
+            evalhyd.evald(obs, prd, ["NSE"],
+                          q_thr=_thr, events=_events)[0]
         )
 
     def test_conditions(self):
@@ -106,8 +139,10 @@ class TestMasking(unittest.TestCase):
             prd = _prd[..., msk].copy()
 
             numpy.testing.assert_almost_equal(
-                evalhyd.evald(_obs, _prd, ["NSE"], m_cdt=cdt)[0],
-                evalhyd.evald(obs, prd, ["NSE"])[0]
+                evalhyd.evald(_obs, _prd, ["NSE"],
+                              q_thr=_thr, events=_events, m_cdt=cdt)[0],
+                evalhyd.evald(obs, prd, ["NSE"],
+                              q_thr=_thr, events=_events)[0]
             )
 
         with self.subTest(conditions="observed streamflow statistics"):
@@ -121,8 +156,10 @@ class TestMasking(unittest.TestCase):
             prd = _prd[..., msk].copy()
 
             numpy.testing.assert_almost_equal(
-                evalhyd.evald(_obs, _prd, ["NSE"], m_cdt=cdt)[0],
-                evalhyd.evald(obs, prd, ["NSE"])[0]
+                evalhyd.evald(_obs, _prd, ["NSE"],
+                              q_thr=_thr, events=_events, m_cdt=cdt)[0],
+                evalhyd.evald(obs, prd, ["NSE"],
+                              q_thr=_thr, events=_events)[0]
             )
 
         with self.subTest(conditions="time indices"):
@@ -138,15 +175,17 @@ class TestMasking(unittest.TestCase):
             prd = _prd[..., 20:].copy()
 
             numpy.testing.assert_almost_equal(
-                evalhyd.evald(_obs, _prd, ["NSE"], m_cdt=cdt)[0],
-                evalhyd.evald(obs, prd, ["NSE"])[0]
+                evalhyd.evald(_obs, _prd, ["NSE"],
+                              q_thr=_thr, events=_events, m_cdt=cdt)[0],
+                evalhyd.evald(obs, prd, ["NSE"],
+                              q_thr=_thr, events=_events)[0]
             )
 
 
 class TestMissingData(unittest.TestCase):
 
     def test_nan(self):
-        for metric in ('RMSE', 'NSE', 'KGE', 'KGEPRIME'):
+        for metric in _all_metrics:
             obs = numpy.array(
                 [[4.7, numpy.nan, 5.5, 2.7, 4.1]]
             )
@@ -155,9 +194,16 @@ class TestMissingData(unittest.TestCase):
                  [numpy.nan, 4.2, 4.7, 4.3, 3.3],
                  [5.3, 5.2, 5.7, numpy.nan, 3.9]]
             )
+            thr = numpy.array(
+                [[4., 5.],
+                 [4., 5.],
+                 [4., 5.]]
+            )
+            events = "low"
 
             with self.subTest(metric=metric):
-                res = evalhyd.evald(obs, prd, [metric])[0]
+                res = evalhyd.evald(obs, prd, [metric],
+                                    q_thr=thr, events=events)[0]
 
                 for i in range(prd.shape[0]):
                     msk = ~numpy.isnan(obs[0]) & ~numpy.isnan(prd[i])
@@ -169,7 +215,8 @@ class TestMissingData(unittest.TestCase):
                         evalhyd.evald(
                             obs[:, msk],
                             prd[i, msk][numpy.newaxis],
-                            [metric]
+                            [metric],
+                            q_thr=thr[i, :][numpy.newaxis], events=events
                         )[0]
                     )
 
@@ -190,21 +237,30 @@ class TestUncertainty(unittest.TestCase):
         obs_3yrs = numpy.hstack((obs_1yr,) * 3)
         prd_3yrs = numpy.hstack((prd_1yr,) * 3)
 
-        for metric in ('RMSE', 'NSE', 'KGE', 'KGEPRIME'):
+        thr = numpy.repeat(
+            numpy.array([[690, 534, 445, numpy.nan]]),
+            repeats=prd_1yr.shape[0], axis=0
+        )
+        events = "low"
+
+        for metric in _all_metrics:
             with self.subTest(metric=metric):
                 numpy.testing.assert_almost_equal(
                     # bootstrap with only one year of data
                     # (compare last sample only to have matching dimensions)
                     evalhyd.evald(
                         obs_1yr, prd_1yr, [metric],
+                        q_thr=thr,
+                        events=events,
                         bootstrap={
                             "n_samples": 10, "len_sample": 3, "summary": 0
                         },
                         dts=dts_1yr
-                    )[0][..., [0]],
+                    )[0][:, :, [0]],
                     # repeat year of data three times to correspond to a
                     # bootstrap sample of length 3
-                    evalhyd.evald(obs_3yrs, prd_3yrs, [metric])[0]
+                    evalhyd.evald(obs_3yrs, prd_3yrs, [metric],
+                                  q_thr=thr, events=events)[0]
                 )