  We use large volume N-body simulations to predict the clustering of dark
matter in redshift space in f(R) modified gravity cosmologies. This is the
first time that the nonlinear matter and velocity fields have been resolved to
such a high level of accuracy over a broad range of scales in this class of
models. We find significant deviations from the clustering signal in standard
gravity, with an enhanced boost in power on large scales and stronger damping
on small scales in the f(R) models compared to GR at redshifts z<1. We measure
the velocity divergence (P_\theta \theta) and matter (P_\delta \delta) power
spectra and find a large deviation in the ratios \sqrt{P_\theta \theta/P_\delta
\delta} and P_\delta \theta/P_\delta\delta, between the f(R) models and GR for
0.03<k/(h/Mpc)<0.5. In linear theory these ratios equal the growth rate of
structure on large scales. Our results show that the simulated ratios agree
with the growth rate for each cosmology (which is scale dependent in the case
of modified gravity) only for extremely large scales, k<0.06h/Mpc at z=0. The
velocity power spectrum is substantially different in the f(R) models compared
to GR, suggesting that this observable is a sensitive probe of modified
gravity. We demonstrate how to extract the matter and velocity power spectra
from the 2D redshift space power spectrum, P(k,\mu), and can recover the
nonlinear matter power spectrum to within a few percent for k<0.1h/Mpc.
However, the model fails to describe the shape of the 2D power spectrum
demonstrating that an improved model is necessary in order to reconstruct the
velocity power spectrum accurately. The same model can match the monopole
moment to within 3% for GR and 10% for the f(R) cosmology at k<0.2 h/Mpc at
z=1. Our results suggest that the extraction of the velocity power spectrum
from future galaxy surveys is a promising method to constrain deviations from
GR.
