diff --git a/include/evalhyd/evald.hpp b/include/evalhyd/evald.hpp
index 248071dcdffbc5f58044813a1ef54848a57be116..cb25cb8ad0d23ba606e318b481a82bd29b5d8c27 100644
--- a/include/evalhyd/evald.hpp
+++ b/include/evalhyd/evald.hpp
@@ -38,15 +38,19 @@ namespace evalhyd
     ///       Streamflow observations. Time steps with missing observations
     ///       must be assigned `NAN` values. Those time steps will be ignored
     ///       both in the observations and the predictions before the
-    ///       *metrics* are computed.
-    ///       shape: ({1+, ...}, time)
+    ///       *metrics* are computed. Observations and predictions must be
+    ///       feature the same number of dimensions and they must be
+    ///       broadcastable.
+    ///       shape: ({... ,} time)
     ///
     ///    q_prd: ``xt::xexpression<A>``
     ///       Streamflow predictions. Time steps with missing predictions
     ///       must be assigned `NAN` values. Those time steps will be ignored
     ///       both in the observations and the predictions before the
-    ///       *metrics* are computed.
-    ///       shape: ({1+, ...}, time)
+    ///       *metrics* are computed. Observations and predictions must be
+    ///       feature the same number of dimensions and they must be
+    ///       broadcastable.
+    ///       shape: ({... ,} time)
     ///
     ///    metrics: ``std::vector<std::string>``
     ///       The sequence of evaluation metrics to be computed.
@@ -130,39 +134,66 @@ namespace evalhyd
             double epsilon=-9
     )
     {
+        auto q_obs_ = q_obs.derived_cast();
+        auto q_prd_ = q_prd.derived_cast();
+
+        // check that observations and predictions dimensions are compatible
+        if (q_prd_.dimension() != q_obs_.dimension())
+        {
+            throw std::runtime_error(
+                    "observations and predictions feature different numbers "
+                    "of dimensions"
+            );
+        }
+
+        auto obs_shp = q_obs_.shape();
+        auto prd_shp = q_prd_.shape();
+
+        if (obs_shp != prd_shp)
+        {
+            // check observations and predictions are broadcastable
+            for (int a = 0; a < obs_shp.size(); a++)
+            {
+                if (obs_shp[a] != prd_shp[a])
+                    if ((obs_shp[a] != 1) and (prd_shp[a] != 1))
+                        throw std::runtime_error(
+                                "observations and predictions are not "
+                                "broadcastable"
+                        );
+            }
+        }
+
+        // apply streamflow transformation if requested
         A obs;
         A prd;
 
-        // apply streamflow transformation if requested
         if ( transform == "none" || (transform == "pow" && exponent == 1))
         {
-            obs = std::move(
-                    q_obs.derived_cast());
-            prd = std::move(
-                    q_prd.derived_cast());
+            obs = std::move(q_obs_);
+            prd = std::move(q_prd_);
         }
         else if ( transform == "sqrt" )
         {
-            obs = xt::sqrt(q_obs.derived_cast());
-            prd = xt::sqrt(q_prd.derived_cast());
+            obs = xt::sqrt(q_obs_);
+            prd = xt::sqrt(q_prd_);
         }
         else if ( transform == "inv" )
         {
             if ( epsilon == -9 )
                 // determine an epsilon value to avoid zero divide
-                epsilon = xt::mean(q_obs.derived_cast())() * 0.01;
+                epsilon = xt::mean(q_obs_)() * 0.01;
 
-            obs = 1. / (q_obs.derived_cast() + epsilon);
-            prd = 1. / (q_prd.derived_cast() + epsilon);
+            obs = 1. / (q_obs_ + epsilon);
+            prd = 1. / (q_prd_ + epsilon);
         }
         else if ( transform == "log" )
         {
             if ( epsilon == -9 )
                 // determine an epsilon value to avoid log zero
-                epsilon = xt::mean(q_obs.derived_cast())() * 0.01;
+                epsilon = xt::mean(q_obs_)() * 0.01;
 
-            obs = xt::log(q_obs.derived_cast() + epsilon);
-            prd = xt::log(q_prd.derived_cast() + epsilon);
+            obs = xt::log(q_obs_ + epsilon);
+            prd = xt::log(q_prd_ + epsilon);
         }
         else if ( transform == "pow" )
         {
@@ -170,15 +201,15 @@ namespace evalhyd
             {
                 if ( epsilon == -9 )
                     // determine an epsilon value to avoid zero divide
-                    epsilon = xt::mean(q_obs.derived_cast())() * 0.01;
+                    epsilon = xt::mean(q_obs_)() * 0.01;
 
-                obs = xt::pow(q_obs.derived_cast() + epsilon, exponent);
-                prd = xt::pow(q_prd.derived_cast() + epsilon, exponent);
+                obs = xt::pow(q_obs_ + epsilon, exponent);
+                prd = xt::pow(q_prd_ + epsilon, exponent);
             }
             else
             {
-                obs = xt::pow(q_obs.derived_cast(), exponent);
-                prd = xt::pow(q_prd.derived_cast(), exponent);
+                obs = xt::pow(q_obs_, exponent);
+                prd = xt::pow(q_prd_, exponent);
             }
         }
         else
diff --git a/include/evalhyd/evalp.hpp b/include/evalhyd/evalp.hpp
index d4efbb9bf3a2ebea690048da10bab56350d36bee..ce9d8e02f2ad737d536c33171d904ab0bfeea100 100644
--- a/include/evalhyd/evalp.hpp
+++ b/include/evalhyd/evalp.hpp
@@ -102,6 +102,32 @@ namespace evalhyd
         // check that optional parameters are given as arguments
         eh::utils::check_optionals(metrics, q_thr);
 
+        // check that data dimensions are compatible
+        // > time
+        if (q_obs.shape(1) != q_prd.shape(3))
+            throw std::runtime_error(
+                    "observations and predictions feature different "
+                    "temporal lengths"
+            );
+        if (t_msk.size() > 0)
+            if (q_obs.shape(1) != t_msk.shape(2))
+                throw std::runtime_error(
+                        "observations and masks feature different "
+                        "temporal lengths"
+                );
+        // > sites
+        if (q_obs.shape(0) != q_prd.shape(0))
+            throw std::runtime_error(
+                    "observations and predictions feature different "
+                    "numbers of sites"
+            );
+        if (t_msk.size() > 0)
+            if (q_obs.shape(0) != t_msk.shape(0))
+                throw std::runtime_error(
+                        "observations and masks feature different "
+                        "numbers of sites"
+                );
+
         // retrieve dimensions
         std::size_t n_sit = q_prd.shape(0);
         std::size_t n_ltm = q_prd.shape(1);
diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index 124b9f06d3bc566a8bf7b761277b2426d0d77e9e..c83591d8e1c4cda9692a71bb5d8df7708a75298b 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -63,24 +63,7 @@ namespace evalhyd
                     // missing (i.e. NaN)
                     q_prd{xt::where(xt::isnan(obs.derived_cast()),
                                     NAN, prd.derived_cast())}
-            {
-                // check that data dimensions are compatible
-                if (_q_prd.dimension() != _q_obs.dimension())
-                {
-                    throw std::runtime_error(
-                            "number of dimensions of 'sim' and 'obs' "
-                            "must be identical"
-                    );
-                }
-                else if (_q_prd.shape(_q_prd.dimension() - 1)
-                         != _q_obs.shape( _q_obs.dimension() - 1))
-                {
-                    throw std::runtime_error(
-                            "length of time series of 'sim' and 'obs' "
-                            "must be identical"
-                    );
-                }
-            };
+            {};
 
             // members for evaluation metrics
             A RMSE;
@@ -108,10 +91,10 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \assign mean_obs:
         //     Mean observed streamflow.
-        //     shape: ({1, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_mean_obs()
         {
@@ -122,10 +105,10 @@ namespace evalhyd
         //
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \assign mean_prd:
         //     Mean predicted streamflow.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_mean_prd()
         {
@@ -136,13 +119,13 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \assign quad_err:
         //     Quadratic errors between observations and predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         template <class A>
         void Evaluator<A>::calc_quad_err()
         {
@@ -153,13 +136,13 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_obs:
         //     Mean observed streamflow.
-        //     shape: ({1, ...}, 1)
+        //     shape: ({... ,} 1)
         // \assign quad_obs:
         //     Quadratic errors between observations and mean observation.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         template <class A>
         void Evaluator<A>::calc_quad_obs()
         {
@@ -170,13 +153,13 @@ namespace evalhyd
         //
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_prd:
         //     Mean predicted streamflow.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \assign quad_prd:
         //     Quadratic errors between predictions and mean prediction.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         template <class A>
         void Evaluator<A>::calc_quad_prd()
         {
@@ -187,25 +170,25 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_obs:
         //     Mean observed streamflow.
-        //     shape: ({1, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_prd:
         //     Mean predicted streamflow.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require quad_obs:
         //     Quadratic errors between observations and mean observation.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require quad_prd:
         //     Quadratic errors between predictions and mean prediction.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \assign r_pearson:
         //     Pearson correlation coefficients.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         template <class A>
         void Evaluator<A>::calc_r_pearson()
         {
@@ -227,13 +210,13 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \assign bias:
         //     Biases.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_bias()
         {
@@ -246,10 +229,10 @@ namespace evalhyd
         //
         // \require quad_err:
         //     Quadratic errors between observations and predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \assign RMSE:
         //     Root-mean-square errors.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_RMSE()
         {
@@ -261,13 +244,13 @@ namespace evalhyd
         //
         // \require quad_err:
         //     Quadratic errors between observations and predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require quad_obs:
         //     Quadratic errors between observations and mean observation.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \assign NSE:
         //     Nash-Sutcliffe efficiencies.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_NSE()
         {
@@ -283,25 +266,25 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_obs:
         //     Mean observed streamflow.
-        //     shape: ({1, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_prd:
         //     Mean predicted streamflow.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require r_pearson:
         //     Pearson correlation coefficients.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require bias:
         //     Biases.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \assign KGE:
         //     Kling-Gupta efficiencies.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_KGE()
         {
@@ -321,25 +304,25 @@ namespace evalhyd
         //
         // \require q_obs:
         //     Streamflow observations.
-        //     shape: ({1, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_obs:
         //     Mean observed streamflow.
-        //     shape: ({1, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require q_prd:
         //     Streamflow predictions.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require mean_prd:
         //     Mean predicted streamflow.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \require r_pearson:
         //     Pearson correlation coefficients.
-        //     shape: ({1+, ...}, time)
+        //     shape: ({... ,} time)
         // \require bias:
         //     Biases.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         // \assign KGEPRIME:
         //     Modified Kling-Gupta efficiencies.
-        //     shape: ({1+, ...}, 1)
+        //     shape: ({... ,} 1)
         template <class A>
         void Evaluator<A>::calc_KGEPRIME()
         {
diff --git a/src/probabilist/evaluator.h b/src/probabilist/evaluator.h
index e7dad5f9ae6f899f26232f0ac5f8d77d228a2e47..c37d0741d2315048f75585bbdb4214c93404f52f 100644
--- a/src/probabilist/evaluator.h
+++ b/src/probabilist/evaluator.h
@@ -83,11 +83,9 @@ namespace evalhyd
                 auto sum_nan = xt::eval(xt::sum(prd_nan, -1));
 
                 if (xt::amin(sum_nan) != xt::amax(sum_nan))
-                {
                     throw std::runtime_error(
                             "predictions across members feature non-equal lengths"
                     );
-                }
 
                 auto msk_nan = xt::where(obs_nan | xt::row(prd_nan, 0))[0];