The field-scale transport of reactive and nonreactive solutes by groundwater in a statistically anisotropic aquifer is examined by means of high-resolution, three-dimensional numerical solutions of the steady state flow and transient advection-dispersion equations. The presence of physical and chemical heterogeneities in the aquifer media is modeled with the use of a geostatistical description of the hydraulic conductivity and chemical distribution coefficient. The geostatistical parameters describing the spatial variations of the log-transformed hydraulic conductivity fields, ln<K(x)>, are chosen to resemble those of the Borden aquifer. For sorbing solutes the spatial variation in the log-transformed distribution coefficient fields, ln<Kd(x)>, is generated so as to have a geometric mean and variance similar to that estimated by Durant (1986) for the organic chemical tetrachlorethane with Borden sand. It is further assumed that the ln<K(x)> and ln<Kd(x)> fields are inversely correlated to each other and that they possess the same spatial correlation structure.

Five realizations of media which incorporate these characteristics are generated by means of the Fourier spectral technique of Robin et al. (1993). It is shown that joint K(x) and Kd(x) variability can impart a large-scale ''pseudokinetic'' behavior, in that the ensemble mean bulk retardation factor can increase with time and plume displacement distance, even though sorption is modeled as being linear and instantaneous at the scale of any single heterogeneity and the flow field is assumed to be steady state. From realization to realization, however, the apparent bulk retardation factor can either increase dramatically at early time in a manner similar to that observed by Roberts et al. (1986) during the Borden tracer test or it can decrease. At large time the ensemble mean velocity of the centroid of the reactive plume is close to a value given by the mean fluid velocity divided by the arithmetic mean of the locally variable retardation factors. Local-scale transport nonidealities, such as intraparticle diffusion and relatively rapid kinetic sorption, are shown to have minimal influence on the plume centroid velocity and macrodispersivity, relative to the effect of aquifer heterogeneity. The results of the numerical simulations further demonstrate that the first-order stochastic analyses of Gelhar and Axness (1983), Dagan (1988), and Naff (1990) tend to overestimate the actual field-scale longitudinal spreading of a nonreactive solute. This result is believed to occur because these analyses include the artificial effect of plume centroid dispersion about the ensemble mean position, in addition to the actual spreading of each plume realization about its center of mass.

When the effect of plume centroid dispersion is added to the numerical simulation results, a reasonable agreement with the solution of Naff (1990), which takes into consideration local-scale dispersion, is obtained for the longitudinal macrodispersivity. In the transverse direction, the ensemble mean spreading agrees reasonably well with the solution of Dagan (1988). It is also shown that the longitudinal macrodispersivity of a reactive solute can be enhanced relative to that of a nonreactive one. Its actual value, however, is overpredicted by the first-order solution of Garabedian (1987) for the reasons given above. Also explored, in a preliminary fashion, are some issues related to prediction uncertainty for both reactive and nonreactive solutes migrating through heterogeneous aquifers.