diff --git a/py/qsotools/mocklib.py b/py/qsotools/mocklib.py
index 67652c3..9dc4598 100644
--- a/py/qsotools/mocklib.py
+++ b/py/qsotools/mocklib.py
@@ -221,6 +221,10 @@ class DLASampler():
     settled on resampling the input Gaussian field, which reduced the excess
     clustering of DLAs.
 
+    This can be further tuned with ``sigma_kms`` option, which adds random
+    velocity separation to DLA redshift obtained from abundance matching. It
+    adds random ``dz = (1 + z_dla) * Normal(0, sigma_kms) / c`` to DLA redshift
+
     Other methods tested:
     - Optical depth mapping on the resampled (0.2 A) skewers. The PDF turns
     out non-trivial when smoothed; and you can't control excess clustering.
@@ -245,13 +249,14 @@ def fn_z_evo(z, zpivot=2.4, gamma_l=1.5):
         return ((1 + z) / (1 + zpivot))**gamma_l
 
     def __init__(
-            self, wide_pix=2**10., nmin=19., nmax=23.,
+            self, wide_pix=2**10., sigma_kms=0, nmin=19., nmax=23.,
             zmin=0, zmax=20., nzbins=5000
     ):
         self.nmin = nmin
         self.nmax = nmax
 
         self.wide_pix = None
+        self.sigma_kms = sigma_kms
         self.sigma_scale = None
         self.set_var_gauss(wide_pix)
 
@@ -406,6 +411,13 @@ def _downsample(x):
             num_dlas = w.sum()
 
             z_dlas = refac_z[w]
+            if self.sigma_kms > 0:
+                # Store and restore random state for backwards reproducibility
+                state = RNST.bit_generator.state
+                z_dlas += (1 + z_dlas) * RNST.normal(
+                    0, self.sigma_kms, size=z_dlas.size) / LIGHT_SPEED
+                RNST.bit_generator.state = state
+
             Nhi_dlas = self.get_random_NHi(num_dlas, RNST)
             id_dlas = np.array([
                 hash(f"{mockids[jj]:020d}{x:02d}") for x in np.arange(num_dlas)
diff --git a/py/qsotools/scripts/generate_mocks.py b/py/qsotools/scripts/generate_mocks.py
index 7bee6b2..a5a2bba 100644
--- a/py/qsotools/scripts/generate_mocks.py
+++ b/py/qsotools/scripts/generate_mocks.py
@@ -101,6 +101,9 @@ def get_parser():
     parser.add_argument(
         "--invcdf-nz", default=PKG_ICDF_Z_TABLE,
         help="Table for inverse cdf of n(z).")
+    parser.add_argument(
+        "--dla-sigma-kms", default=0, type=float,
+        help="Random shift central DLA redshift.")
 
     parser.add_argument(
         "--chunk-dyn", action="store_true",
@@ -590,7 +593,7 @@ def main():
         args.z_forest_min = 0
         args.keep_nolya_pixels = True
         args.save_full_flux = True
-        dla_sampler = lm.DLASampler()
+        dla_sampler = lm.DLASampler(sigma_kms=args.dla_sigma_kms)
     else:
         dla_sampler = None