diff --git a/README.md b/README.md
index c786d34b21cc55efab7d95ab74d061908b09ab63..2800e9fc24fb3b00f9385618af28395296c99ca2 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,3 @@
 # evalhyd
 
-Utility to evaluate streamflow simulations/forecasts
+Utility to evaluate deterministic and probabilistic streamflow predictions.
diff --git a/include/evalhyd/determinist.hpp b/include/evalhyd/determinist.hpp
index 6b4bc9f64c66550b6fe9f57b1bd8bcbd87e43b9f..19dce0afb95ecb5a12fb2c403dab2f7e44b42bd9 100644
--- a/include/evalhyd/determinist.hpp
+++ b/include/evalhyd/determinist.hpp
@@ -12,14 +12,14 @@ namespace eh = evalhyd;
 
 namespace evalhyd
 {
-    /// Function allowing the evaluation of streamflow simulations using a
-    /// range of relevant metrics.
+    /// Function allowing the evaluation of deterministic streamflow
+    /// predictions using a range of relevant metrics.
     ///
     /// \param [in] q_obs:
     ///     2D array of streamflow observations.
     ///     shape: (1+, time)
-    /// \param [in] q_sim:
-    ///     2D array of streamflow simulations.
+    /// \param [in] q_prd:
+    ///     2D array of streamflow predictions.
     ///     shape: (1+, time)
     /// \param [in] metrics:
     ///     Vector of strings for the metric(s) to be computed.
@@ -28,12 +28,12 @@ namespace evalhyd
     template <class A>
     std::vector<A> evald(
             const xt::xexpression<A>& q_obs,
-            const xt::xexpression<A>& q_sim,
+            const xt::xexpression<A>& q_prd,
             const std::vector<std::string>& metrics
     )
     {
         const A& obs = q_obs.derived_cast();
-        const A& sim = q_sim.derived_cast();
+        const A& prd = q_prd.derived_cast();
 
         // check that the metrics to be computed are valid
         utils::check_metrics(
@@ -42,7 +42,7 @@ namespace evalhyd
         );
 
         // instantiate determinist evaluator
-        eh::determinist::Evaluator<A> evaluator(obs, sim);
+        eh::determinist::Evaluator<A> evaluator(obs, prd);
 
         // declare maps for memoisation purposes
         std::unordered_map<std::string, std::vector<std::string>> elt;
diff --git a/include/evalhyd/determinist/evaluator.hpp b/include/evalhyd/determinist/evaluator.hpp
index e2e42e12ebbde202f7ca349089366b9839f3f707..d31abdd03ae68d02ec651660e4d0383a81ffd599 100644
--- a/include/evalhyd/determinist/evaluator.hpp
+++ b/include/evalhyd/determinist/evaluator.hpp
@@ -13,15 +13,15 @@ namespace evalhyd
         private:
             // members for input data
             const A& q_obs;
-            const A& q_sim;
+            const A& q_prd;
             // members for computational elements
             // TODO
 
         public:
             // constructor method
-            Evaluator(const A& obs, const A& sim) : q_obs{obs}, q_sim{sim} {
+            Evaluator(const A& obs, const A& prd) : q_obs{obs}, q_prd{prd} {
                 // check that data dimensions are compatible
-                if (q_sim.dimension() != q_obs.dimension())
+                if (q_prd.dimension() != q_obs.dimension())
                 {
                     throw std::runtime_error(
                             "number of dimensions of 'sim' and 'obs' "
@@ -39,7 +39,7 @@ namespace evalhyd
 
         // Compute the Nash-Sutcliffe Efficiency (NSE).
         //
-        // If multi-dimensional arrays are provided, the arrays of simulations
+        // If multi-dimensional arrays are provided, the arrays of predictions
         // and observations must feature the same number of dimensions, they must
         // be broadcastable, and their temporal dimensions must be along their
         // last axis.
@@ -47,8 +47,8 @@ namespace evalhyd
         // \require q_obs:
         //     Array of streamflow observations.
         //     shape: (1, time)
-        // \require q_sim:
-        //     Array of streamflow simulations.
+        // \require q_prd:
+        //     Array of streamflow predictions.
         //     shape: (..., time)
         // \assign nse
         //     Array of computed Nash-Sutcliffe efficiencies.
@@ -60,7 +60,7 @@ namespace evalhyd
             A q_avg = xt::mean(q_obs, -1, xt::keep_dims);
 
             // compute squared errors operands
-            A f_num = xt::sum(xt::square(q_obs - q_sim), -1, xt::keep_dims);
+            A f_num = xt::sum(xt::square(q_obs - q_prd), -1, xt::keep_dims);
             A f_den = xt::sum(xt::square(q_obs - q_avg), -1, xt::keep_dims);
 
             // return computed NSE
diff --git a/include/evalhyd/probabilist.hpp b/include/evalhyd/probabilist.hpp
index 8bf5e1921a0abe143d89ded7609257971724a9c1..98c68a911648c4413d79176e346104b88fe703e0 100644
--- a/include/evalhyd/probabilist.hpp
+++ b/include/evalhyd/probabilist.hpp
@@ -14,14 +14,14 @@ namespace eh = evalhyd;
 
 namespace evalhyd
 {
-    /// Function allowing the evaluation of streamflow forecasts using a
-    /// range of relevant metrics.
+    /// Function allowing the evaluation of probabilistic streamflow
+    /// predictions using a range of relevant metrics.
     ///
     /// \param [in] q_obs:
     ///     2D array of streamflow observations.
     ///     shape: (1, time)
-    /// \param [in] q_frc:
-    ///     2D array of streamflow forecasts.
+    /// \param [in] q_prd:
+    ///     2D array of streamflow predictions.
     ///     shape: (members, time)
     /// \param [in] metrics:
     ///     Vector of strings for the metric(s) to be computed.
@@ -32,7 +32,7 @@ namespace evalhyd
     ///     Vector of 2D array of metrics for each threshold.
     std::vector<xt::xtensor<double, 2>> evalp(
             const xt::xtensor<double, 2>& q_obs,
-            const xt::xtensor<double, 2>& q_frc,
+            const xt::xtensor<double, 2>& q_prd,
             const std::vector<std::string>& metrics,
             const xt::xtensor<double, 1>& q_thr = {}
     )
@@ -47,7 +47,7 @@ namespace evalhyd
         eh::utils::check_optionals(metrics, q_thr);
 
         // instantiate probabilist evaluator
-        eh::probabilist::Evaluator evaluator(q_obs, q_frc, q_thr);
+        eh::probabilist::Evaluator evaluator(q_obs, q_prd, q_thr);
 
         // declare maps for memoisation purposes
         std::unordered_map<std::string, std::vector<std::string>> elt;
diff --git a/include/evalhyd/probabilist/evaluator.h b/include/evalhyd/probabilist/evaluator.h
index 5fe8af7deafc6263c87463c4151a617d50f7cc12..ccd60af7ec42280a9ded16b4fec52249beda84c9 100644
--- a/include/evalhyd/probabilist/evaluator.h
+++ b/include/evalhyd/probabilist/evaluator.h
@@ -12,7 +12,7 @@ namespace evalhyd
         private:
             // members for input data
             const xt::xtensor<double, 2>& q_obs;
-            const xt::xtensor<double, 2>& q_frc;
+            const xt::xtensor<double, 2>& q_prd;
             const xt::xtensor<double, 1>& q_thr;
 
             // members for computational elements
@@ -24,9 +24,9 @@ namespace evalhyd
         public:
             // constructor method
             Evaluator(const xt::xtensor<double, 2>& obs,
-                      const xt::xtensor<double, 2>& frc,
+                      const xt::xtensor<double, 2>& prd,
                       const xt::xtensor<double, 1>& thr) :
-                    q_obs{obs}, q_frc{frc}, q_thr{thr} {};
+                    q_obs{obs}, q_prd{prd}, q_thr{thr} {};
 
             // members for intermediate evaluation metrics
             // (i.e. before the reduction along the temporal axis)
diff --git a/include/evalhyd/probabilist/evaluator_brier.cpp b/include/evalhyd/probabilist/evaluator_brier.cpp
index 339dcb7e1cb7f18825c9848fdea20b1a2b53cc1a..dfe611a842e29ec50764aa2e2b1da80bc58bcc42 100644
--- a/include/evalhyd/probabilist/evaluator_brier.cpp
+++ b/include/evalhyd/probabilist/evaluator_brier.cpp
@@ -81,7 +81,7 @@ namespace evalhyd
             // define some dimensions
             std::size_t n = o_k.shape(1);
             std::size_t n_thr = o_k.shape(0);
-            std::size_t n_mbr = q_frc.shape(0);
+            std::size_t n_mbr = q_prd.shape(0);
             std::size_t n_cmp = 3;
 
             // initialise output variable
diff --git a/include/evalhyd/probabilist/evaluator_elements.cpp b/include/evalhyd/probabilist/evaluator_elements.cpp
index 4b7bc731aa6ec3987a2c0526a806dc98b8af6978..1618d8f0fa194bdd6a22a2e25490aa5fe6674239 100644
--- a/include/evalhyd/probabilist/evaluator_elements.cpp
+++ b/include/evalhyd/probabilist/evaluator_elements.cpp
@@ -41,8 +41,8 @@ namespace evalhyd
 
         // Determine forecast probability of threshold(s) exceedance to occur.
         //
-        // \require q_frc:
-        //     2D array of streamflow forecasts.
+        // \require q_prd:
+        //     2D array of streamflow predictions.
         //     shape: (members, time)
         // \require q_thr:
         //     1D array of streamflow exceedance threshold(s).
@@ -54,7 +54,7 @@ namespace evalhyd
         {
             // determine if members have exceeded threshold(s)
             xt::xtensor<double, 3> e_frc =
-                    is_above_threshold(q_frc, q_thr);
+                    is_above_threshold(q_prd, q_thr);
 
             // calculate how many members have exceeded threshold(s)
             xt::xtensor<double, 2> n_frc = xt::sum(e_frc, 1);
@@ -63,7 +63,7 @@ namespace evalhyd
             // and forecast probabilities of exceedance
             // /!\ probability calculation dividing by n (the number of
             //     members), not n+1 (the number of ranks) like in other metrics
-            std::size_t n_mbr = q_frc.shape(0);
+            std::size_t n_mbr = q_prd.shape(0);
 
             // determine probability of threshold(s) exceedance
             y_k = n_frc / n_mbr;
@@ -71,15 +71,15 @@ namespace evalhyd
 
         // Compute the forecast quantiles from the ensemble members.
         //
-        // \require q_frc:
-        //     2D array of streamflow forecasts.
+        // \require q_prd:
+        //     2D array of streamflow predictions.
         //     shape: (members, time)
         // \assign q_qnt:
         //     2D array of streamflow forecast quantiles.
         //     shape: (quantiles, time)
         void Evaluator::calc_q_qnt()
         {
-            q_qnt = xt::sort(q_frc, 0);
+            q_qnt = xt::sort(q_prd, 0);
         }
     }
 }
diff --git a/include/evalhyd/probabilist/evaluator_quantiles.cpp b/include/evalhyd/probabilist/evaluator_quantiles.cpp
index 6b227ac081e91f8dc504596d24202638c6ce8f46..89b54ce453207c03655c24cce6520eeb7a0878ce 100644
--- a/include/evalhyd/probabilist/evaluator_quantiles.cpp
+++ b/include/evalhyd/probabilist/evaluator_quantiles.cpp
@@ -17,7 +17,7 @@ namespace evalhyd
         //     2D array of streamflow observations.
         //     shape: (1, time)
         // \require q_qnt:
-        //     2D array of streamflow forecasts.
+        //     2D array of streamflow quantiles.
         //     shape: (quantiles, time)
         // \assign qs:
         //     2D array of the quantile scores for each time step.
diff --git a/tests/data/q_frc.csv b/tests/data/q_prd.csv
similarity index 100%
rename from tests/data/q_frc.csv
rename to tests/data/q_prd.csv
diff --git a/tests/test_determinist.cpp b/tests/test_determinist.cpp
index 50aca4387b8817d98aa3be5b6d780d11ef7d80d2..f9bc1b9728fe8e7d5a38ba118fba29ca004bc2a6 100644
--- a/tests/test_determinist.cpp
+++ b/tests/test_determinist.cpp
@@ -18,23 +18,23 @@ TEST(DeterministTests, TestNSE) {
 
     xt::xtensor<double, 1> observed_1d = xt::squeeze(observed_2d);
 
-    ifs.open("./data/q_frc.csv");
-    xt::xtensor<double, 2> forecast_2d = xt::view(
+    ifs.open("./data/q_prd.csv");
+    xt::xtensor<double, 2> predicted_2d = xt::view(
             xt::transpose(xt::load_csv<double>(ifs)), xt::range(0, 5), xt::all()
     );
     ifs.close();
 
-    xt::xtensor<double, 1> forecast_1d = xt::row(forecast_2d, 0);
+    xt::xtensor<double, 1> predicted_1d = xt::row(predicted_2d, 0);
 
     // compute scores (both with 2D and 1D tensors)
     std::vector<xt::xtensor<double, 2>> metrics_2d =
             evalhyd::evald<xt::xtensor<double, 2>>(
-                    observed_2d, forecast_2d, {"NSE"}
+                    observed_2d, predicted_2d, {"NSE"}
             );
 
     std::vector<xt::xtensor<double, 1>> metrics_1d =
             evalhyd::evald<xt::xtensor<double, 1>>(
-                    observed_1d, forecast_1d, {"NSE"}
+                    observed_1d, predicted_1d, {"NSE"}
             );
 
     // check results (both with 2D and 1D tensors)
diff --git a/tests/test_probabilist.cpp b/tests/test_probabilist.cpp
index 49c7a749b8a593b5027af95fca2d589338b820e1..f63e4096af27c18727dc290eec933a22e039da90 100644
--- a/tests/test_probabilist.cpp
+++ b/tests/test_probabilist.cpp
@@ -15,8 +15,8 @@ TEST(ProbabilistTests, TestBrier) {
     xt::xtensor<double, 2> observed = xt::load_csv<int>(ifs);
     ifs.close();
 
-    ifs.open("./data/q_frc.csv");
-    xt::xtensor<double, 2> forecast = xt::load_csv<double>(ifs);
+    ifs.open("./data/q_prd.csv");
+    xt::xtensor<double, 2> predicted = xt::load_csv<double>(ifs);
     ifs.close();
 
     // compute scores
@@ -24,7 +24,7 @@ TEST(ProbabilistTests, TestBrier) {
 
     std::vector<xt::xtensor<double, 2>> metrics =
             evalhyd::evalp(
-                    xt::transpose(observed), xt::transpose(forecast),
+                    xt::transpose(observed), xt::transpose(predicted),
                     {"BS", "BSS", "BS_CRD", "BS_LBD"},
                     thresholds
             );