diff --git a/rex/utilities/bc_utils.py b/rex/utilities/bc_utils.py
index daf32ca2..220928d9 100644
--- a/rex/utilities/bc_utils.py
+++ b/rex/utilities/bc_utils.py
@@ -107,7 +107,8 @@ class QuantileDeltaMapping:
 
     def __init__(self, params_oh, params_mh, params_mf, dist='empirical',
                  relative=True, sampling='linear', log_base=10,
-                 delta_denom_min=None, delta_denom_zero=None):
+                 delta_denom_min=None, delta_denom_zero=None,
+                 delta_range=None):
         """
         Parameters
         ----------
@@ -167,6 +168,11 @@ def __init__(self, params_oh, params_mh, params_mf, dist='empirical',
             division by a very small number making delta blow up and resulting
             in very large output bias corrected values. See equation 4 of
             Cannon et al., 2015 for the delta term.
+        delta_range : tuple | None
+            Option to set a (min, max) on the delta term in QDM. This can help
+            prevent QDM from making non-realistic increases/decreases in
+            otherwise physical values. See equation 4 of Cannon et al., 2015
+            for the delta term.
         """
 
         self.params_oh = params_oh
@@ -179,6 +185,7 @@ def __init__(self, params_oh, params_mh, params_mf, dist='empirical',
         self.scipy_dist = None
         self.delta_denom_min = delta_denom_min
         self.delta_denom_zero = delta_denom_zero
+        self.delta_range = delta_range
 
         if self.dist_name != 'empirical':
             self.scipy_dist = getattr(scipy.stats, self.dist_name, None)
@@ -303,7 +310,7 @@ def ppf(cls, p, params, scipy_dist, sampling, log_base):
     @classmethod
     def run_qdm(cls, arr, params_oh, params_mh, params_mf,
                 scipy_dist, relative, sampling, log_base, delta_denom_min,
-                delta_denom_zero):
+                delta_denom_zero, delta_range):
         """Run the actual QDM operation from args without initializing the
         ``QuantileDeltaMapping`` object
 
@@ -364,6 +371,11 @@ def run_qdm(cls, arr, params_oh, params_mh, params_mf,
             division by a very small number making delta blow up and resulting
             in very large output bias corrected values. See equation 4 of
             Cannon et al., 2015 for the delta term.
+        delta_range : tuple | None
+            Option to set a (min, max) on the delta term in QDM. This can help
+            prevent QDM from making non-realistic increases/decreases in
+            otherwise physical values. See equation 4 of Cannon et al., 2015
+            for the delta term.
 
         Returns
         -------
@@ -401,9 +413,16 @@ def run_qdm(cls, arr, params_oh, params_mh, params_mf,
             if delta_denom_min is not None:
                 x_mh_mf = np.maximum(x_mh_mf, delta_denom_min)
             delta = arr / x_mh_mf  # Eq.4: x_m_p / F-1_m_h(Tau_m_p)
+            if delta_range is not None:
+                delta = np.maximum(delta, np.min(delta_range))
+                delta = np.minimum(delta, np.max(delta_range))
             arr_bc = x_oh * delta  # Eq.6: x^_m_p = x^_o:m_h:p * delta
+
         else:
             delta = arr - x_mh_mf  # Eq.4: x_m_p - F-1_m_h(Tau_m_p)
+            if delta_range is not None:
+                delta = np.maximum(delta, np.min(delta_range))
+                delta = np.minimum(delta, np.max(delta_range))
             arr_bc = x_oh + delta  # Eq.6: x^_m_p = x^_o:m_h:p + delta
 
         return arr_bc
@@ -432,7 +451,8 @@ def __call__(self, arr, max_workers=1):
             arr_bc = self.run_qdm(arr, self.params_oh, self.params_mh,
                                   self.params_mf, self.scipy_dist,
                                   self.relative, self.sampling, self.log_base,
-                                  self.delta_denom_min, self.delta_denom_zero)
+                                  self.delta_denom_min, self.delta_denom_zero,
+                                  self.delta_range)
         else:
             max_workers = max_workers or os.cpu_count()
             sslices = np.array_split(np.arange(arr.shape[1]), arr.shape[1])
@@ -448,7 +468,7 @@ def __call__(self, arr, max_workers=1):
                                      self.params_mf[idx], self.scipy_dist,
                                      self.relative, self.sampling,
                                      self.log_base, self.delta_denom_min,
-                                     self.delta_denom_zero)
+                                     self.delta_denom_zero, self.delta_range)
                     futures[fut] = idx
                 for future in as_completed(futures):
                     idx = futures[future]
diff --git a/rex/version.py b/rex/version.py
index 6be23b51..3d0eb7c4 100644
--- a/rex/version.py
+++ b/rex/version.py
@@ -1,3 +1,3 @@
 """rex Version number"""
 
-__version__ = "0.2.90"
+__version__ = "0.2.91"
diff --git a/tests/test_bc.py b/tests/test_bc.py
index 2e210633..12b99eed 100644
--- a/tests/test_bc.py
+++ b/tests/test_bc.py
@@ -159,3 +159,11 @@ def test_difficult_qdm():
                                        delta_denom_min=0.01)
     arr_bc = qdm_rel_fut(arr[:, np.newaxis])
     assert arr_bc.max() < 40  # fixed result
+
+    # test that delta_range fixes this
+    qdm_rel_fut = QuantileDeltaMapping(params_oh, params_mh, params_mf,
+                                       dist='empirical', relative=True,
+                                       sampling='invlog', log_base=10,
+                                       delta_range=(0.5, 1.5))
+    arr_bc = qdm_rel_fut(arr[:, np.newaxis])
+    assert arr_bc.max() < (arr.max() * 1.5)  # fixed result