diff --git a/deps/evalhyd b/deps/evalhyd
index a2957cb3f6e674df6cf5e148870f74c9d53251b4..25ca41bdfcb85cbbef072b635727077df0326f01 160000
--- a/deps/evalhyd
+++ b/deps/evalhyd
@@ -1 +1 @@
-Subproject commit a2957cb3f6e674df6cf5e148870f74c9d53251b4
+Subproject commit 25ca41bdfcb85cbbef072b635727077df0326f01
diff --git a/src/evalhyd-python.cpp b/src/evalhyd-python.cpp
index 9635b12800989f6735251098a71fed4e87603b01..2b822e1a28dc25b28c23141021fa38c914dbf11d 100644
--- a/src/evalhyd-python.cpp
+++ b/src/evalhyd-python.cpp
@@ -166,11 +166,19 @@ PYBIND11_MODULE(evalhyd, m)
             :Parameters:
 
                 q_obs: `numpy.ndarray`
-                    2D array of streamflow observations.
+                    2D array of streamflow observations. Time steps with
+                    missing observation must be assigned `numpy.nan`
+                    values. Those time steps will be pairwise ignored in
+                    the observations and the predictions before the
+                    *metrics* are computed.
                     shape: (sites, time)
 
                 q_prd: `numpy.ndarray`
-                    4D array of streamflow predictions.
+                    4D array of streamflow predictions. Time steps with
+                    missing prediction must be assigned `numpy.nan`
+                    values. Those time steps will be pairwise ignored in
+                    the observations and the predictions before the
+                    *metrics* are computed.
                     shape: (sites, lead times, members, time)
 
                 metrics: `List[str]`
diff --git a/tests/test_probabilist.py b/tests/test_probabilist.py
index 6dfc37f2dad9afabb30b1783dbb6dd2c21470da9..a59c9b86a664a4d184cd23adc6d7fad0e0c774c2 100644
--- a/tests/test_probabilist.py
+++ b/tests/test_probabilist.py
@@ -87,6 +87,34 @@ class TestMasking(unittest.TestCase):
         )
 
 
+class TestMissingData(unittest.TestCase):
+
+    def test_nan(self):
+        thr = numpy.array([[690, 534, 445, numpy.nan]])
+        for metric in ("BS", "BSS", "BS_CRD", "BS_LBD", "QS", "CRPS"):
+            with self.subTest(metric=metric):
+                numpy.testing.assert_almost_equal(
+                    # missing data flagged as NaN
+                    evalhyd.evalp(
+                        [[4.7, numpy.nan, 5.5, 2.7, 4.1]],
+                        [[[[5.3, 4.2, 5.7, 2.3, numpy.nan],
+                           [4.3, 4.2, 4.7, 4.3, numpy.nan],
+                           [5.3, 5.2, 5.7, 2.3, numpy.nan]]]],
+                        [metric],
+                        thr
+                    )[0],
+                    # missing data pairwise deleted from series
+                    evalhyd.evalp(
+                        [[4.7, 5.5, 2.7]],
+                        [[[[5.3, 5.7, 2.3],
+                           [4.3, 4.7, 4.3],
+                           [5.3, 5.7, 2.3]]]],
+                        [metric],
+                        thr
+                    )[0]
+                )
+
+
 if __name__ == '__main__':
     test_loader = unittest.TestLoader()
     test_suite = unittest.TestSuite()
@@ -100,6 +128,9 @@ if __name__ == '__main__':
     test_suite.addTests(
         test_loader.loadTestsFromTestCase(TestMasking)
     )
+    test_suite.addTests(
+        test_loader.loadTestsFromTestCase(TestMissingData)
+    )
 
     runner = unittest.TextTestRunner(verbosity=2)
     result = runner.run(test_suite)