Abstract The velocity divergence power spectrum is a key ingredient in modelling redshift-space distortion effects on quasi-linear and non-linear scales. We present an improved model for the z=0 velocity divergence auto and cross power spectrum which was originally suggested by Jennings et al. Using numerical simulations we measure the velocity fields using a Delaunay tessellation and obtain an accurate prediction of the velocity divergence power spectrum on scales k < 1 h Mpc−1. We use this to update the model which is now accurate to 2 per cent for both Pθθ and Pθδ at z= 0 on scales k < 0.65 h Mpc−1 and k < 0.35 h Mpc−1, respectively. We find that the formula for the redshift dependence of the velocity divergence power spectra proposed by Jennings et al. recovers the measured z > 0 P(k) to markedly greater accuracy with the new model. The non-linear Pθθ and Pθδ at z=1 are recovered accurately to better than 2 per cent on scales k < 0.2 h Mpc−1. Recently, it was shown that the velocity field shows larger differences between modified gravity cosmologies and Λ cold dark matter (ΛCDM) compared to the matter field. An accurate model for the velocity divergence power spectrum, such as the one presented here, is a valuable tool for analysing redshift-space distortion effects in future galaxy surveys and for constraining deviations from general relativity.