diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py
index 48fe7e25..bc19701e 100755
--- a/src/gstools/field/base.py
+++ b/src/gstools/field/base.py
@@ -578,7 +578,7 @@ def pos(self, pos):
             self._pos, self._field_shape = format_struct_pos_dim(pos, self.dim)
         # prepend dimension if we have a vector field
         if self.value_type == "vector":
-            self._field_shape = (self.dim,) + self._field_shape
+            self._field_shape = (self.generator.vec_dim,) + self._field_shape
             if self.latlon:
                 raise ValueError("Field: Vector fields not allowed for latlon")
 
diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py
index 60e327b3..549f0da2 100644
--- a/src/gstools/field/generator.py
+++ b/src/gstools/field/generator.py
@@ -413,6 +413,8 @@ class IncomprRandMeth(RandMeth):
         the mean velocity in x-direction
     mode_no : :class:`int`, optional
         number of Fourier modes. Default: ``1000``
+    vec_dim : :class:`int`, optional
+        vector dimension, in case it mismatches the model dimension
     seed : :class:`int` or :any:`None`, optional
         the seed of the random number generator.
         If "None", a random seed is used. Default: :any:`None`
@@ -465,18 +467,27 @@ def __init__(
         model,
         mean_velocity=1.0,
         mode_no=1000,
+        vec_dim=None,
         seed=None,
         verbose=False,
         sampling="auto",
         **kwargs,
     ):
-        if model.dim < 2 or model.dim > 3:
+        if vec_dim is None and (model.dim < 2 or model.dim > 3):
             raise ValueError(
-                "Only 2D and 3D incompressible fields can be generated."
+                "Only 2D and 3D incompressible vectors can be generated."
+            )
+        if vec_dim is not None and (vec_dim < 2 or vec_dim > 3):
+            raise ValueError(
+                "Only 2D and 3D incompressible vectors can be generated."
             )
         super().__init__(model, mode_no, seed, verbose, sampling, **kwargs)
 
         self.mean_u = mean_velocity
+        if vec_dim is None:
+            self.vec_dim = model.dim
+        else:
+            self.vec_dim = vec_dim
         self._value_type = "vector"
 
     def __call__(self, pos, add_nugget=True):
@@ -501,7 +512,7 @@ def __call__(self, pos, add_nugget=True):
         """
         pos = np.asarray(pos, dtype=np.double)
         summed_modes = summate_incompr(
-            self._cov_sample, self._z_1, self._z_2, pos
+            self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos
         )
         nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0
         e1 = self._create_unit_vector(summed_modes.shape)
@@ -532,7 +543,7 @@ def _create_unit_vector(self, broadcast_shape, axis=0):
             the unit vector
         """
         shape = np.ones(len(broadcast_shape), dtype=int)
-        shape[0] = self.model.dim
+        shape[0] = self.vec_dim
 
         e1 = np.zeros(shape)
         e1[axis] = 1.0
diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx
index ecd7ea58..aaf76170 100644
--- a/src/gstools/field/summator.pyx
+++ b/src/gstools/field/summator.pyx
@@ -50,6 +50,7 @@ cdef (double) abs_square(const double[:] vec) nogil:
 
 
 def summate_incompr(
+    const int vec_dim,
     const double[:, :] cov_samples,
     const double[:] z_1,
     const double[:] z_2,
@@ -58,24 +59,24 @@ def summate_incompr(
     cdef int i, j, d
     cdef double phase
     cdef double k_2
-    cdef int dim = pos.shape[0]
+    cdef int field_dim = pos.shape[0]
 
-    cdef double[:] e1 = np.zeros(dim, dtype=float)
+    cdef double[:] e1 = np.zeros(vec_dim, dtype=float)
     e1[0] = 1.
-    cdef double[:] proj = np.empty(dim)
+    cdef double[:] proj = np.empty(vec_dim)
 
     cdef int X_len = pos.shape[1]
     cdef int N = cov_samples.shape[1]
 
-    cdef double[:, :] summed_modes = np.zeros((dim, X_len), dtype=float)
+    cdef double[:, :] summed_modes = np.zeros((vec_dim, X_len), dtype=float)
 
     for i in range(X_len):
         for j in range(N):
-            k_2 = abs_square(cov_samples[:, j])
+            k_2 = abs_square(cov_samples[:vec_dim, j])
             phase = 0.
-            for d in range(dim):
+            for d in range(field_dim):
                 phase += cov_samples[d, j] * pos[d, i]
-            for d in range(dim):
+            for d in range(vec_dim):
                 proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2
                 summed_modes[d, i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))