diff --git a/src/gstools/field/generator.py b/src/gstools/field/generator.py
index 549f0da2..50dd56d0 100644
--- a/src/gstools/field/generator.py
+++ b/src/gstools/field/generator.py
@@ -10,6 +10,8 @@
    Generator
    RandMeth
    IncomprRandMeth
+   IncomprRandZeroVelMeth
+   GenericRandVectorFieldMeth
 """
 # pylint: disable=C0103, W0222, C0412, W0231
 import warnings
@@ -26,9 +28,11 @@
     # pylint: disable=E0401
     from gstools_core import summate, summate_incompr
 else:
-    from gstools.field.summator import summate, summate_incompr
+    from gstools.field.summator import summate, summate_incompr, summate_incompr_zero_vel, summate_generic_vector_field
+
+
+__all__ = ["RandMeth", "IncomprRandMeth", "IncomprRandZeroVelMeth", "GenericRandVectorFieldMeth", ]
 
-__all__ = ["Generator", "RandMeth", "IncomprRandMeth"]
 
 
 SAMPLING = ["auto", "inversion", "mcmc"]
@@ -176,10 +180,10 @@ def __init__(
     ):
         if kwargs:
             warnings.warn("gstools.RandMeth: **kwargs are ignored")
-        # initialize attributes
+        # initialize atributes
         self._mode_no = int(mode_no)
         self._verbose = bool(verbose)
-        # initialize private attributes
+        # initialize private atributes
         self._model = None
         self._seed = None
         self._rng = None
@@ -252,6 +256,7 @@ def update(self, model=None, seed=np.nan):
             If :any:`None`, a random seed is used. If :any:`numpy.nan`,
             the actual seed will be kept. Default: :any:`numpy.nan`
         """
+
         # check if a new model is given
         if isinstance(model, CovModel):
             if self.model != model:
@@ -311,12 +316,14 @@ def reset_seed(self, seed=np.nan):
         if seed is None or not np.isnan(seed):
             self._seed = seed
         self._rng = RNG(self._seed)
+
         # normal distributed samples for randmeth
         self._z_1 = self._rng.random.normal(size=self._mode_no)
         self._z_2 = self._rng.random.normal(size=self._mode_no)
+            
         # sample uniform on a sphere
         sphere_coord = self._rng.sample_sphere(self.model.dim, self._mode_no)
-        # sample radii according to radial spectral density of the model
+        # sample radii acording to radial spectral density of the model
         if self.sampling == "inversion" or (
             self.sampling == "auto" and self.model.has_ppf
         ):
@@ -324,12 +331,14 @@ def reset_seed(self, seed=np.nan):
             rad = self._rng.sample_dist(
                 size=self._mode_no, pdf=pdf, cdf=cdf, ppf=ppf, a=0
             )
+
         else:
             rad = self._rng.sample_ln_pdf(
                 ln_pdf=self.model.ln_spectral_rad_pdf,
                 size=self._mode_no,
                 sample_around=1.0 / self.model.len_rescaled,
             )
+
         # get fully spatial samples by multiplying sphere samples and radii
         self._cov_sample = rad * sphere_coord
 
@@ -474,12 +483,14 @@ def __init__(
         **kwargs,
     ):
         if vec_dim is None and (model.dim < 2 or model.dim > 3):
+
             raise ValueError(
                 "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)
 
@@ -517,7 +528,7 @@ def __call__(self, pos, add_nugget=True):
         nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0
         e1 = self._create_unit_vector(summed_modes.shape)
         return (
-            self.mean_u * e1
+            #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields
             + self.mean_u
             * np.sqrt(self.model.var / self._mode_no)
             * summed_modes
@@ -548,3 +559,302 @@ def _create_unit_vector(self, broadcast_shape, axis=0):
         e1 = np.zeros(shape)
         e1[axis] = 1.0
         return e1
+    
+class IncomprRandZeroVelMeth(RandMeth):
+    r"""RandMeth for incompressible random vector fields with zero velocity, using Eq. 20 of Kraichnan (1970)
+
+    Parameters
+    ----------
+    model : :any:`CovModel`
+        covariance model
+    mean_velocity : :class:`float`, optional
+        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`
+    verbose : :class:`bool`, optional
+        State if there should be output during the generation.
+        Default: :any:`False`
+    sampling : :class:`str`, optional
+        Sampling strategy. Either
+
+            * "auto": select best strategy depending on given model
+            * "inversion": use inversion method
+            * "mcmc": use mcmc sampling
+
+    **kwargs
+        Placeholder for keyword-args
+
+    Notes
+    -----
+    The Randomization method is used to generate isotropic
+    spatial incompressible random vector fields characterized
+    by a given covariance model. The equation is [Kraichnan1970]_:
+
+    .. math::
+       u\left(x\right)= 
+       \bar{u}\sqrt{\frac{\sigma^{2}}{N}}\cdot
+       \sum_{j=1}^{N}\left(
+       W_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+
+       W_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right)
+       \right)
+
+    where:
+
+        * :math:`\bar u` : mean velocity in :math:`e_1` direction
+        * :math:`N` : fourier mode number
+        * :math:`Z_{k,j}` : random samples from a normal distribution, each of size (N,vec_dim)
+        * :math:`k_j` : samples from the spectral density distribution of
+          the covariance model
+        * :math:`W_{k,j}` : cross-product of the random normal vector Z_{k,j} with the wave-vector k_j
+          ensuring the incompressibility
+
+    References
+    ----------
+    .. [Kraichnan1970] Kraichnan, R. H.,
+           "Diffusion by a random velocity field.",
+           The physics of fluids, 13(1), 22-31., (1970)
+    """
+
+    def __init__(
+        self,
+        model,
+        mean_velocity=1.0,
+        mode_no=1000,
+        vec_dim=None,
+        seed=None,
+        verbose=False,
+        sampling="auto",
+        periodic_bc=False,
+        box_len=None,
+        **kwargs,
+    ):
+        if vec_dim is None and (model.dim < 2 or model.dim > 3):
+            raise ValueError(
+                "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"
+        
+        self.periodic_bc=periodic_bc
+        self.box_len=box_len
+
+        # for using the Kraichnan method for zero-velocity fluid, z_1 and z_2 must instead contain N=self._mode_no independent realizations of normal vectors of size vec_dim,
+        mean = np.zeros(self.vec_dim)
+        cov = np.identity(self.vec_dim)
+
+        self._z_1 = self._rng.random.multivariate_normal(mean, cov, size=self._mode_no) # shape (_mode_no, self.vec_dim)
+        self._z_2 = self._rng.random.multivariate_normal(mean, cov, size=self._mode_no)
+        print("shape of z_1: ",np.shape(self._z_1))
+        
+    def __call__(self, pos):
+        """Calculate the random modes for the randomization method.
+
+        This method  calls the `summate_incompr_*` Cython methods,
+        which are the heart of the randomization method.
+        In this class the method contains a projector to
+        ensure the incompressibility of the vector field.
+
+        Parameters
+        ----------
+        pos : (d, n), :class:`numpy.ndarray`
+            the position tuple with d dimensions and n points.
+
+        Returns
+        -------
+        :class:`numpy.ndarray`
+            the random modes
+        """
+        pos = np.asarray(pos, dtype=np.double)
+        
+        print("\nStarting summate_incompr_zero_vel")
+
+        if self.periodic_bc:
+            if not self.box_len==None:
+                print("\nimposing periodic boundary conditions on spatial coordinates!")
+                print("\nrounding first 'vec_dim' vector components in cov_sample to multiples of 2pi/box_length")
+                fac = 2*np.pi/np.array(self.box_len)
+                self._cov_sample[:self.vec_dim,:] = fac[:,None]*np.round(self._cov_sample[:self.vec_dim,:]/fac[:,None])
+            else:
+                raise ValueError(
+                "For periodic boundary conditions on spatial coordinates, specify parameter box_len, an array of lengths of the box in each spatial dimension. The length of box_len must be equal to vec_dim."
+                )
+                
+        
+        summed_modes = summate_incompr_zero_vel(
+            self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos
+        )
+        print("\nFinished mode summation!")
+
+        nugget = self.get_nugget(summed_modes.shape)
+
+        return (
+            #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields
+            + self.mean_u
+            * np.sqrt(self.model.var / self._mode_no)
+            * summed_modes
+            + nugget
+        )
+
+
+    
+class GenericRandVectorFieldMeth(RandMeth):
+    r"""RandMeth for incompressible random vector fields.
+
+    Parameters
+    ----------
+    model : :any:`CovModel`
+        covariance model
+    mean_velocity : :class:`float`, optional
+        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`
+    verbose : :class:`bool`, optional
+        State if there should be output during the generation.
+        Default: :any:`False`
+    sampling : :class:`str`, optional
+        Sampling strategy. Either
+
+            * "auto": select best strategy depending on given model
+            * "inversion": use inversion method
+            * "mcmc": use mcmc sampling
+
+    **kwargs
+        Placeholder for keyword-args
+
+    Notes
+    -----
+    The Randomization method is used to generate isotropic
+    spatial incompressible random vector fields characterized
+    by a given covariance model. The equation is [Kraichnan1970]_:
+
+    .. math::
+       u_i\left(x\right)= \bar{u_i} \delta_{i1} +
+       \bar{u_i}\sqrt{\frac{\sigma^{2}}{N}}\cdot
+       \sum_{j=1}^{N}p_i(k_{j})\left(
+       Z_{1,j}\cdot\cos\left(\left\langle k_{j},x\right\rangle \right)+
+       Z_{2,j}\cdot\sin\left(\left\langle k_{j},x\right\rangle \right)
+       \right)
+
+    where:
+
+        * :math:`\bar u` : mean velocity in :math:`e_1` direction
+        * :math:`N` : fourier mode number
+        * :math:`Z_{k,j}` : random samples from a normal distribution
+        * :math:`k_j` : samples from the spectral density distribution of
+          the covariance model
+        * :math:`p_i(k_j) = e_1 - \frac{k_i k_1}{k^2}` : the projector
+          ensuring the incompressibility
+
+    References
+    ----------
+    .. [Kraichnan1970] Kraichnan, R. H.,
+           "Diffusion by a random velocity field.",
+           The physics of fluids, 13(1), 22-31., (1970)
+    """
+
+    def __init__(
+        self,
+        model,
+        mean_velocity=1.0,
+        mode_no=1000,
+        vec_dim=None,
+        seed=None,
+        verbose=False,
+        sampling="auto",
+        **kwargs,
+    ):
+        if vec_dim is None and (model.dim < 2 or model.dim > 3):
+            raise ValueError(
+                "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):
+        """Calculate the random modes for the randomization method.
+
+        This method  calls the `summate_incompr_*` Cython methods,
+        which are the heart of the randomization method.
+        In this class the method contains a projector to
+        ensure the incompressibility of the vector field.
+
+        Parameters
+        ----------
+        pos : (d, n), :class:`numpy.ndarray`
+            the position tuple with d dimensions and n points.
+
+        Returns
+        -------
+        :class:`numpy.ndarray`
+            the random modes
+        """
+        pos = np.asarray(pos, dtype=np.double)
+        summed_modes = summate_generic_vector_field(
+            self.vec_dim, self._cov_sample, self._z_1, self._z_2, pos
+        )
+        nugget = self.get_nugget(summed_modes.shape)
+
+        e1 = self._create_unit_vector(summed_modes.shape)
+
+        return (
+            #self.mean_u * e1 #!!! Joshua has commented this out to get zero-velocity fields
+            + self.mean_u
+            * np.sqrt(self.model.var / self._mode_no)
+            * summed_modes
+            + nugget
+        )
+
+    def _create_unit_vector(self, broadcast_shape, axis=0):
+        """Create a unit vector.
+
+        Can be multiplied with a vector of shape broadcast_shape
+
+        Parameters
+        ----------
+        broadcast_shape : :class:`tuple`
+            the shape of the array with which
+            the unit vector is to be multiplied
+        axis : :class:`int`, optional
+            the direction of the unit vector. Default: ``0``
+
+        Returns
+        -------
+        :class:`numpy.ndarray`
+            the unit vector
+        """
+        shape = np.ones(len(broadcast_shape), dtype=int)
+        shape[0] = self.vec_dim
+
+        e1 = np.zeros(shape)
+        e1[axis] = 1.0
+        return e1
\ No newline at end of file
diff --git a/src/gstools/field/srf.py b/src/gstools/field/srf.py
index 06839779..270173ec 100644
--- a/src/gstools/field/srf.py
+++ b/src/gstools/field/srf.py
@@ -14,7 +14,7 @@
 import numpy as np
 
 from gstools.field.base import Field
-from gstools.field.generator import Generator, IncomprRandMeth, RandMeth
+from gstools.field.generator import IncomprRandMeth, IncomprRandZeroVelMeth, RandMeth, GenericRandVectorFieldMeth 
 from gstools.field.upscaling import var_coarse_graining, var_no_scaling
 
 __all__ = ["SRF"]
@@ -22,8 +22,10 @@
 GENERATOR = {
     "RandMeth": RandMeth,
     "IncomprRandMeth": IncomprRandMeth,
+    "IncomprRandZeroVelMeth": IncomprRandZeroVelMeth,
     "VectorField": IncomprRandMeth,
     "VelocityField": IncomprRandMeth,
+    "GenericRandVectorFieldMeth": GenericRandVectorFieldMeth,
 }
 """dict: Standard generators for spatial random fields."""
 
@@ -120,7 +122,7 @@ def __call__(
             the position tuple, containing main direction and transversal
             directions
         seed : :class:`int`, optional
-            seed for RNG for resetting. Default: keep seed from generator
+            seed for RNG for reseting. Default: keep seed from generator
         point_volumes : :class:`float` or :class:`numpy.ndarray`
             If your evaluation points for the field are coming from a mesh,
             they are probably representing a certain element volume.
diff --git a/src/gstools/field/summator.pyx b/src/gstools/field/summator.pyx
index aaf76170..b49da3b3 100644
--- a/src/gstools/field/summator.pyx
+++ b/src/gstools/field/summator.pyx
@@ -81,3 +81,84 @@ def summate_incompr(
                 summed_modes[d, i] += proj[d] * (z_1[j] * cos(phase) + z_2[j] * sin(phase))
 
     return np.asarray(summed_modes)
+
+
+# function makes incompressible random vector field with zero velocity, unlike summate_incompr which makes the x-axis a preferential direction
+# z_1 has shape (vec_dim,) such that z_1[i] for i =0,...,vec_dim-1 is a 1D array of length N = cov_samples.shape[1]
+def summate_incompr_zero_vel(
+    const int vec_dim,
+    const double[:, :] cov_samples,
+    const double[:,:] z_1,
+    const double[:,:] z_2,
+    const double[:, :] pos
+    ):
+    cdef int i, j, d
+    cdef double phase
+    cdef int field_dim = pos.shape[0]
+
+    cdef int X_len = pos.shape[1]
+    cdef int N = cov_samples.shape[1]
+
+    cdef double[:, :] summed_modes = np.zeros((vec_dim, X_len), dtype=float)
+    cdef double[:] w_1 = np.empty(vec_dim, dtype=float) # !!! Joshua hopefully initialized w, v correctly, this is just for memory allocation right?
+    cdef double[:] w_2 = np.empty(vec_dim, dtype=float)
+
+    for i in range(X_len):
+        for j in range(N):
+            # compute cross-product of wavevectors with random vectors z_1,z_2 of size vec_dim drawn from (vec_dim)-dimensional Gaussian with unit covariance matrix
+            
+            # !!! using the following two lines is cleaner, but code takes forever (type instability?) WHY?????
+            
+            #w_1 = np.cross(z_1[j,:], cov_samples[:vec_dim, j]) # 1D array of size vec_dim
+            #w_2 = np.cross(z_2[j,:], cov_samples[:vec_dim, j])
+            
+            # Pending trying to get above two lines to work, I am hard coding the cross-product for vec_dim =3
+            # !!! WARNING: this only works for vec_dim = 3 
+            w_1[0] = z_1[j,1] * cov_samples[2, j] - z_1[j,2] * cov_samples[1, j]
+            w_1[1] = z_1[j,2] * cov_samples[0, j] - z_1[j,0] * cov_samples[2, j]
+            w_1[2] = z_1[j,0] * cov_samples[1, j] - z_1[j,1] * cov_samples[0, j]
+            
+            w_2[0] = z_2[j,1] * cov_samples[2, j] - z_2[j,2] * cov_samples[1, j]
+            w_2[1] = z_2[j,2] * cov_samples[0, j] - z_2[j,0] * cov_samples[2, j]
+            w_2[2] = z_2[j,0] * cov_samples[1, j] - z_2[j,1] * cov_samples[0, j]
+
+            phase = 0.
+            for d in range(field_dim):
+                phase += cov_samples[d, j] * pos[d, i]
+            for d in range(vec_dim):
+                summed_modes[d, i] +=  w_1[d] * cos(phase) + w_2[d] * sin(phase) #1 * cos(phase) + 1 * sin(phase)
+
+    return np.asarray(summed_modes)
+
+def summate_generic_vector_field(
+    const int vec_dim,
+    const double[:, :] cov_samples,
+    const double[:] z_1,
+    const double[:] z_2,
+    const double[:, :] pos
+    ):
+    cdef int i, j, d
+    cdef double phase
+    cdef double k_2
+    cdef int field_dim = pos.shape[0]
+
+    cdef double[:] e1 = np.zeros(vec_dim, dtype=float)
+    e1[0] = 1.
+    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((vec_dim, X_len), dtype=float)
+
+    for i in range(X_len):
+        for j in range(N):
+            k_2 = abs_square(cov_samples[:vec_dim, j])
+            phase = 0.
+            for d in range(field_dim):
+                phase += cov_samples[d, j] * pos[d, i]
+            for d in range(vec_dim):
+                #proj[d] = e1[d] - cov_samples[d, j] * cov_samples[0, j] / k_2 #!!! don't want incompressibility projector  here
+                summed_modes[d, i] += (z_1[j] * cos(phase) + z_2[j] * sin(phase))
+
+    return np.asarray(summed_modes)