To study statistics of velocity fields in three-dimensional heterogeneous multi-Gaussian saturated hydraulic conductivity fields and the accuracy of their prediction, we performed high-resolution Monte Carlo (MC) analyses. The MC analyses included variances of the log hydraulic conductivity in the range of 0.5 ≤ σY2 ≤ 3.0 and anisotropy ratios in the range of 0.017 ≤ e ≤ 1. The statistics of the velocity fields from the MC analyses are compared with analytical solutions of the first- and second-order approximations of the stochastic flow equation. This paper shows that the second-order approximations fit significantly better to the univariate statistics of the Darcy velocity from the MC analyses. For isotropic cases the second-order approximations correspond fairly well to the univariate statistics of the velocity. For anisotropic cases the accordance is given only for the mean velocity and the variance of the transverse vertical component of the velocity. The MC analyses show that the spatial correlation of the velocity decreases more rapidly with increasing σY2. This was more pronounced for the anisotropic than for the isotropic case. The negative correlations, in absolute terms, of the transverse velocity components simultaneously decrease with increasing σY2. This is in contrast to the first-order approximation of the spatial correlations of the velocity. It is assumed that the discrepancies between approximate solutions of the stochastic flow equation and the results of the MC analyses are strongly dependent on the nonnormality of the probability density distributions of the velocity.