diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py
index 7a0148cab..b381003e3 100644
--- a/pysteps/blending/steps.py
+++ b/pysteps/blending/steps.py
@@ -90,7 +90,7 @@ def forecast(
     conditional=False,
     probmatching_method="cdf",
     mask_method="incremental",
-    resample_distribution=False,
+    resample_distribution=True,
     smooth_radar_mask_range=0,
     callback=None,
     return_output=True,
@@ -222,7 +222,7 @@ def forecast(
         Method to resample the distribution from the extrapolation and NWP cascade as input
         for the probability matching. Not resampling these distributions may lead to losing
         some extremes when the weight of both the extrapolation and NWP cascade is similar.
-        Defaults to False.
+        Defaults to True.
     smooth_radar_mask_range: int, Default is 0.
       Method to smooth the transition between the radar-NWP-noise blend and the NWP-noise
       blend near the edge of the radar domain (radar mask), where the radar data is either
@@ -1550,9 +1550,9 @@ def worker(j):
                         # that for the probability matching
                         if probmatching_method is not None and resample_distribution:
                             R_pm_resampled = probmatching.resample_distributions(
-                                extrapolation_cascade=R_pm_ep[t_index],
-                                model_cascade=precip_models_pm_temp[j],
-                                weight_extrapolation=weights_pm_normalized[0],
+                                first_array=R_pm_ep[t_index],
+                                second_array=precip_models_pm_temp[j],
+                                probability_first_array=weights_pm_normalized[0],
                             )
                         else:
                             R_pm_resampled = R_pm_blended.copy()
diff --git a/pysteps/postprocessing/probmatching.py b/pysteps/postprocessing/probmatching.py
index b650a261f..12d670cc2 100644
--- a/pysteps/postprocessing/probmatching.py
+++ b/pysteps/postprocessing/probmatching.py
@@ -260,42 +260,42 @@ def _get_error(scale):
     return shift, scale, R.reshape(shape)
 
 
-def resample_distributions(extrapolation_cascade, model_cascade, weight_extrapolation):
+def resample_distributions(first_array, second_array, probability_first_array):
     """
     Merges two distributions (e.g. from the extrapolation nowcast and NWP in the blending module)
     to effectively combine two distributions for the probability matching without losing extremes.
     Parameters
     ----------
-    extrapolation_cascade: array_like
+    first_array: array_like
         One of the two arrays from which the distribution should be samples (e.g. the extrapolation
         cascade).
-    model_cascade: array_like
+    second_array: array_like
         One of the two arrays from which the distribution should be samples (e.g. the NWP (model)
         cascade).
-    weight_extrapolation: float
-        The weight that extrapolation_cascade should get (as a value between 0 and 1). This is
+    probability_first_array: float
+        The weight that first_array should get (as a value between 0 and 1). This is
         typically based on cascade level 2 (as is done for all the probability matching steps in `blending/steps.py`)
 
     Returns
     ----------
     csort: array_like
-        The output distribution as extrapolation_cascade drawn binomial distribution from the input arrays extrapolation_cascade and model_cascade.
+        The output distribution as first_array drawn binomial distribution from the input arrays first_array and second_array.
     """
     # First make sure there are no nans in the input data
     # Where the extrapolation cascade is nan (outside the radar domain), we fill it up with the model data
-    nan_indices = np.isnan(model_cascade)
-    model_cascade[nan_indices] = np.nanmin(model_cascade)
-    nan_indices = np.isnan(extrapolation_cascade)
-    extrapolation_cascade[nan_indices] = model_cascade[nan_indices]
+    nan_indices = np.isnan(second_array)
+    second_array[nan_indices] = np.nanmin(second_array)
+    nan_indices = np.isnan(first_array)
+    first_array[nan_indices] = second_array[nan_indices]
 
     # Prepare the input distributions
-    assert extrapolation_cascade.size == model_cascade.size
-    asort = np.sort(extrapolation_cascade.flatten())[::-1]
-    bsort = np.sort(model_cascade.flatten())[::-1]
+    assert first_array.size == second_array.size
+    asort = np.sort(first_array.flatten())[::-1]
+    bsort = np.sort(second_array.flatten())[::-1]
     n = asort.shape[0]
 
     # Resample the distributions
-    idxsamples = np.random.binomial(1, weight_extrapolation, n).astype(bool)
+    idxsamples = np.random.binomial(1, probability_first_array, n).astype(bool)
     csort = bsort.copy()
     csort[idxsamples] = asort[idxsamples]
     csort = np.sort(csort)[::-1]