The integrodifferential equation for the spatial second moments X of the ensemble mean concentration in a heterogeneous aquifer is nonlinear due to statistical dependence of the particle displacement on X. This nonlinear equation is either linearized or quasi-linearized in previous studies to derive the linear and quasi-linear theories of time-dependent macrodispersion in aquifers. In this study a fully nonlinear analysis is carried out by solving the integrodifferential equation for X numerically and iteratively. The effects of the variance of log hydraulic conductivities &sgr;2Y, the local P¿clet number Pe, and the anisotropic ratio e are then investigated. Results show that in both statistically isotropic and anisotropic media, as compared with the linear theories, the effect of nonlinearity in X is to reduce the spatial longitudinal variance, X11, and enhance the transverse spreading of a solute plume except in isotropic media with &sgr;2Y≤1, where the linear theories may underestimate the longitudinal spreading of a solute. It is also shown that the effect of local dispersion on X11 can be neglected when Pe≥10 but on the transverse macrodispersion, this effect is significant for Pe as large as 100. Nevertheless, the effect of Pe on macrodispersion is secondary as compared with the effect of nonlinearity in X. Application of the nonlinear results shows good fits to the observed spatial variances of tracer concentration in the Borden experiment and excellent agreement with the simulated variances from recent Monte Carlo simulations. ¿ American Geophysical Union 1995