diff --git a/include/evalhyd/probabilistic.hpp b/include/evalhyd/probabilistic.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..880742fdef49b27e68775690e12b076813987f84
--- /dev/null
+++ b/include/evalhyd/probabilistic.hpp
@@ -0,0 +1,115 @@
+#ifndef EVALHYD_PROBABILISTIC_HPP
+#define EVALHYD_PROBABILISTIC_HPP
+
+#include <xtensor/xexpression.hpp>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xindex_view.hpp>
+#include <xtensor/xio.hpp>
+
+#include "utils.hpp"
+
+namespace eh = evalhyd;
+
+namespace evalhyd
+{
+ /// Compute the Brier Score (BS).
+ ///
+ /// \param [in] q_obs:
+ /// 2D array of streamflow observations.
+ /// \param [in] q_frc:
+ /// 2D array of streamflow forecasts.
+ /// \param [in] q_thr:
+ /// 1D array of streamflow exceedance threshold(s).
+ /// \param [in] axis (optional):
+ /// Array axis along which the temporal dimension is for both
+ /// *frc* and *obs* (either 0 or 1). If not provided, set to
+ /// default value 0.
+ /// \return
+ /// 1D array of Brier Score for each threshold.
+ xt::xtensor<double, 1> bs(
+ const xt::xtensor<double, 2>& q_obs,
+ const xt::xtensor<double, 2>& q_frc,
+ const xt::xtensor<double, 1>& q_thr,
+ int axis = 0
+ )
+ {
+ // check that axis is a valid index for a 2D array
+ if ((axis != 0) && (axis != 1))
+ {
+ throw std::out_of_range(
+ "array axis must be 0 or 1 for 2D arrays"
+ );
+ }
+
+ // determine observed realisation of threshold(s) exceedance
+ xt::xtensor<double, 2> p_obs =
+ xt::squeeze(eh::utils::is_above_threshold(q_obs, q_thr));
+
+ // determine forecast probability of threshold(s) exceedance
+ // > determine if members have exceeded threshold(s)
+ xt::xtensor<double, 3> e_frc =
+ eh::utils::is_above_threshold(q_frc, q_thr);
+
+ // > calculate how many members have exceeded threshold(s)
+ int axis_mbr = (axis == 0) ? 1 : 0;
+ xt::xtensor<int, 2> n_frc = xt::sum(e_frc, axis_mbr);
+
+ // > determine correspondence between number of exceeding members
+ // and forecast probabilities of exceedance assuming members are
+ // equiprobable
+ std::size_t n_mbr = q_frc.shape(axis_mbr);
+ xt::xtensor<double, 1> p_mbr =
+ xt::arange<double>(double (n_mbr + 1)) / (n_mbr + 1);
+
+ // > replace number exceeding members with corresponding probability
+ xt::xtensor<double, 2> p_frc = xt::zeros<double>(n_frc.shape());
+ for (int i = 0; i < q_thr.size(); i++)
+ {
+ xt::col(p_frc, i) = xt::index_view(p_mbr, xt::col(n_frc, i));
+ }
+
+ // return computed BS
+ return xt::mean(xt::square(p_obs - p_frc), std::vector<int> { axis });
+ }
+
+ /// Compute the Brier Skill Score (BSS).
+ ///
+ /// \param [in] q_obs:
+ /// 2D array of streamflow observations.
+ /// \param [in] q_frc:
+ /// 2D array of streamflow forecasts.
+ /// \param [in] q_thr:
+ /// 1D array of streamflow exceedance threshold(s).
+ /// \param [in] axis (optional):
+ /// Array axis along which the temporal dimension is for both
+ /// *frc* and *obs* (either 0 or 1). If not provided, set to
+ /// default value 0.
+ /// \return
+ /// 1D array of Brier Skill Score for each threshold.
+ xt::xtensor<double, 1> bss(
+ const xt::xtensor<double, 2>& q_obs,
+ const xt::xtensor<double, 2>& q_frc,
+ const xt::xtensor<double, 1>& q_thr,
+ int axis = 0
+ )
+ {
+ // determine observed realisation of threshold(s) exceedance
+ xt::xtensor<double, 2> p_obs =
+ xt::squeeze(eh::utils::is_above_threshold(q_obs, q_thr));
+
+ // calculate mean observed realisation of threshold(s) exceedance
+ xt::xtensor<double, 1> p_obs_mean =
+ xt::mean(p_obs, std::vector<int> { axis });
+
+ // calculate reference Brier Score
+ xt::xtensor<double, 1> bs_ref = xt::mean(
+ xt::square(p_obs - p_obs_mean), std::vector<int> { axis }
+ );
+
+ // return computed BSS
+ return 1 - (bs(q_obs, q_frc, q_thr, axis) / bs_ref);
+ }
+
+}
+
+#endif //EVALHYD_PROBABILISTIC_HPP
diff --git a/include/evalhyd/utils.hpp b/include/evalhyd/utils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..016e986fbc51adba8f7192e85cd240147381a27b
--- /dev/null
+++ b/include/evalhyd/utils.hpp
@@ -0,0 +1,80 @@
+#ifndef EVALHYD_UTILS_HPP
+#define EVALHYD_UTILS_HPP
+
+#include <xtensor/xtensor.hpp>
+#include <xtensor/xbroadcast.hpp>
+#include <xtensor/xadapt.hpp>
+#include <xtensor/xview.hpp>
+#include <xtensor/xmath.hpp>
+
+namespace evalhyd
+{
+ namespace utils
+ {
+ /// Determine whether flows are strictly greater than given threshold(s)
+ ///
+ /// Streamflow data is always expected as a 2D array (even if only for
+ /// one time series). Threshold(s) given is (are) always expected as a
+ /// 1D array (even if only one threshold is given). A 3D array is
+ /// returned whose first two dimensions are of the same sizes as the
+ /// streamflow data and whose third dimension is of size equal to the
+ /// number of thresholds given. The returned array contains ones where
+ /// the threshold is strictly exceeded, and zeros otherwise.
+ ///
+ /// \param [in] q:
+ /// 2D array of streamflow data
+ /// \param [in] thr:
+ /// 1D array of streamflow threshold(s)
+ /// \return
+ /// 3D array of ones and zeros
+ xt::xtensor<double, 3> is_above_threshold(
+ const xt::xtensor<double, 2>& q,
+ const xt::xtensor<double, 1>& thr
+ )
+ {
+ // add one dimension to the streamflow data
+ auto new_shape = {q.shape(0), q.shape(1), thr.size()};
+ auto q2 = xt::broadcast(
+ xt::view(q, xt::all(), xt::all(), xt::newaxis()),
+ new_shape
+ );
+
+ // return a boolean-like array
+ return q2 > thr;
+ }
+
+ /// Determine whether flows are strictly lower than given threshold(s)
+ ///
+ /// Streamflow data is always expected as a 2D array (even if only for
+ /// one time series). Threshold(s) given is (are) always expected as a
+ /// 1D array (even if only one threshold is given). A 3D array is
+ /// returned whose first two dimensions are of the same sizes as the
+ /// streamflow data and whose third dimension is of size equal to the
+ /// number of thresholds given. The returned array contains ones where
+ /// the threshold is strictly not exceeded, and zeros otherwise.
+ ///
+ /// \param [in] q:
+ /// 2D array of streamflow data
+ /// \param [in] thr:
+ /// 1D array of streamflow threshold(s)
+ /// \return
+ /// 3D array of ones and zeros
+ xt::xtensor<double, 3> is_below_threshold(
+ const xt::xtensor<double, 2>& q,
+ const xt::xtensor<double, 1>& thr
+ )
+ {
+ // add one dimension to the streamflow data
+ auto new_shape = {q.shape(0), q.shape(1), thr.size()};
+ auto q2 = xt::broadcast(
+ xt::view(q, xt::all(), xt::all(), xt::newaxis()),
+ new_shape
+ );
+
+ // return a boolean-like array
+ return q2 < thr;
+ }
+ }
+}
+
+#endif //EVALHYD_UTILS_HPP