diff --git a/src/determinist/evaluator.hpp b/src/determinist/evaluator.hpp
index da85e3fbe4228f80e5e21d399cb041f010062528..ad130f9ab7dd050849d3ffacf87c43d3bd24035e 100644
--- a/src/determinist/evaluator.hpp
+++ b/src/determinist/evaluator.hpp
@@ -6,6 +6,19 @@
 
 namespace evalhyd
 {
+
+    namespace detail {
+        // function to calculate standard deviation on last axis
+        //
+        // TODO: substitute with `xt::stddev` when fixed for `xt::rtensor`
+        //       (see https://gitlab.irstea.fr/HYCAR-Hydro/evalhyd/evalhyd-r/-/issues/1)
+        template <class A1, class A2>
+        inline auto std_dev(A1&& arr, A2&& mean_arr)
+        {
+            return xt::sqrt(xt::mean(xt::square(arr - mean_arr), -1, xt::keep_dims));
+        }
+    }
+
     namespace determinist
     {
         template <class A>
@@ -255,9 +268,8 @@ namespace evalhyd
         void Evaluator<A>::calc_KGE()
         {
             // calculate error in spread of flow $alpha$
-            auto alpha =
-                    xt::stddev(q_prd, {q_prd.dimension() - 1}, xt::keep_dims)
-                    / xt::stddev(q_obs, {q_obs.dimension() - 1}, xt::keep_dims);
+            auto alpha = detail::std_dev(q_prd, mean_prd)
+                         / detail::std_dev(q_obs, mean_obs);
 
             // compute KGE
             KGE = 1 - xt::sqrt(
@@ -294,9 +306,8 @@ namespace evalhyd
         void Evaluator<A>::calc_KGEPRIME()
         {
             // calculate error in spread of flow $gamma$
-            auto gamma =
-                    (xt::stddev(q_prd, {q_prd.dimension() - 1}, xt::keep_dims) / mean_prd)
-                    / (xt::stddev(q_obs, {q_obs.dimension() - 1}, xt::keep_dims) / mean_obs);
+            auto gamma = (detail::std_dev(q_prd, mean_prd) / mean_prd)
+                         / (detail::std_dev(q_obs, mean_obs) / mean_obs);
 
             // compute KGE
             KGEPRIME = 1 - xt::sqrt(