Stochastic analyses and Monte Carlo simulations were conducted for nonergodic transport of a nonreactive solute plume in three-dimensional heterogeneous and statistically anisotropic aquifers under uniform mean flow along the x axis. The hydraulic conductivity, K(x), is modeled as a random field which is assumed to be lognormally distributed with an anisotropic exponential covariance. The simulation model is validated with good comparison of the simulated and theoretical variogram model of the log K field, small mass balance errors, and a large number of Monte Carlo (MC) runs. The ensemble averages of the second spatial moments of a solute plume about its center of mass, Zii, and the plume centroid variances, Rii, (i = 1, 2, 3) of 1600 MC runs were simulated for three degrees of heterogeneity and a line source of three different lengths along the y or z axis. At least 800 MC runs are needed to stabilize the simulated moments even for the mildly heterogeneous aquifers of σY2 ≤ 0.5, and the ergodic condition is far from reaching for a line source with the initial dimension of several integral scales of log K field. The first-order theories predict well the longitudinal Z11 and R11 (thus the one-particle displacement variance, X11 = Z11 + R11), significantly overestimate the horizontal transverse Z22 and R22 (thus X22), and overestimate the vertical transverse Z33 but underestimate R33.