The complexity of mass transfer processes often complicates solute transport simulations. We present a new approach for the implementation of the multirate mass transfer model into random walk particle tracking. This novel method allows for a spatially heterogeneous distribution of mass transfer coefficients as well as hydrodynamic parameters in three dimensions, and it is well suited for avoiding numerical dispersion and solving computationally demanding transport simulations. For this purpose the normalized zeroth spatial moments of the multirate transport equations are derived and used as phase transition probabilities. Performing a simple Bernoulli trial on the appropriate phase transition probabilities the particle distribution between the mobile domain and any immobile domain can be determined. The approach is compared satisfactorily to analytical and semianalytical solutions for one-dimensional, advective-dispersive transport with different types of mass transfer. Aspects of the numerical implementation of this approach are presented, and it is demonstrated that two restrictive criteria for the time step size have to be considered. By adjusting the time step size for each grid cell on the basis of the cell specific velocity field and mass transfer rate a correct simulation of solute transport is assured, while at the same time computational efficiency is preserved. Finally, an example is presented evaluating the effect of a heterogeneous intraparticle pore diffusion in a synthetic aquifer. The results demonstrate that for this specific case the heterogeneous distribution of mass transfer rates does not have a significant influence on mean solute transport behavior but that at low concentration ranges, differences between the different mass transfer models become visible.