Geological formations are ubiquitously heterogeneous, and the equations that govern flow and transport in such formations can be treated as stochastic partial differential equations. The Monte Carlo method is a straightforward approach for simulating flow in heterogeneous porous media; an alternative based on the moment-equation approach has been developed in the last two decades to reduce the high computational expense required by the Monte Carlo method. However, the computational cost of the moment-equation approach is still high. For example, to solve head covariance up to first order in terms of $\sigma_Y^2$, the variance of log hydraulic conductivity Y = ln Ks, it is required to solve sets of linear algebraic equations with N unknowns for 2N times (N being the number of grid nodes). The cost is even higher if higher-order approximations are needed. Zhang and Lu [J. Comput. Phys., 194 (2004), pp. 773--794] developed a new approach to evaluate high-order moments (fourth order for mean head in terms of $\sigma_Y$, and third order for head variances in terms of $\sigma_Y^2$) of flow quantities based on the combination of Karhunen--Loeve decomposition and perturbation methods. In this study, we systematically investigate the computational efficiency and solution accuracy of three approaches: Monte Carlo simulations, the conventional moment-equation (CME) approach, and the moment-equation approach based on Karhunen--Loeve decomposition (KLME). It is evident that the computational cost for the KLME approach is significantly lower than those required by the Monte Carlo and CME approaches. More importantly, while the computational costs (in terms of the number of times for solving linear algebraic equations with N unknowns) for the CME approach depend on the number of grid nodes, the cost for the KLME approach is independent of the number of grid nodes. This makes it possible to apply the KLME method to solve more realistic large-scale flow problems.