Attribution: You must attribute the work in the manner specified by the author or licensor (but no in any way that suggests that they endorse you or your use of the work).Noncommercial ‒ you may not use this work for commercial purpose.No Derivative works ‒ You may not alter, transform, or build upon this work.Sharing ‒ Individual scientists are hereby granted permission, without fees or further requests to GSA, to use a single figure, a single table, and/or a brief paragraph of text in other subsequent works and to make unlimited photocopies of items in this journal for noncommercial use in classrooms to further education and science.

Abstract

This article presents a genesis method for characterizing heterogeneous media representing alluvial deposits. This method simulates the main steps of the medium genesis for meandering, braided, and incising streams and generates facies, which are then translated into hydraulic conductivities to simulate flow and transport. In order to compare this “genetic” model with other methods commonly used to characterize heterogeneous media, a basic sequential Gaussian indicator method was applied to the same site: a 5200-m-long reach of the Aube River floodplain (France). Ten different geostatistical realizations were generated. An equivalent homogeneous representation was also included. Flow and transport simulations in the different heterogeneous numerical media were conducted with Visual MODFLOW. The results were analyzed and compared in terms of permeability fields, plume spreading, and equivalent longitudinal dispersion. Emphasis is on the ability of the genetic model to represent continuous channels that can serve either as conduits or as barriers to flow, which, we think, is unique.

INTRODUCTION

How to represent heterogeneous media in groundwater flow and transport simulations has been the subject of intensive research for many years; see, e.g., a review in Marsily et al. (2005). One of the difficulties is caused by the concept of connectivity of high- (or low-) permeability lenses or strata that may act as conduits (or barriers) to flow. In the stochastic approach, Zinn and Harvey (2003) showed that classical statistical tools such as the probability distribution and the covariance function of the permeability field are insufficient to correctly characterize the connectivity of high- or low-permeability facies and therefore to correctly represent flow and transport. As will be discussed later, this conclusion has been challenged by other authors (e.g., Proce et al., 2004), who have proposed new geostatistical concepts to better represent connectivity.

This article focuses on a special type of aquifer, alluvial deposits. To better describe the connectivity of facies, one can use a “genetic” approach, which mimics the deposition mechanisms that created the alluvium. Thus, this approach is able to generate patterns of “channels” as a result of the incising and filling-up of river beds. The channels may act as preferential paths or barriers to the groundwater flow depending on the type of deposits (coarse or fine material). In the case of a meandering system, simulated channels have a sinusoidal shape, which represents a paleo–stream bed with perfect connectivity unless interrupted by another incision-infilling sequence. Trying to describe their morphology and connectivity within a simple geostatistical framework by a covariance function, which only considers correlations as a function of distance, is likely not to be very effective: a covariance can describe ellipsoidal morphologies, but not sinusoidal ones. In general, geostatistics applied to alluvial media concludes that the permeability correlation distance is very short, on the order of tens of meters, because the method cannot infer the presence of these long-range continuous sinusoidal structures even if anisotropy is assumed while processing the data.

Our comparison is limited to these two approaches: a new “genetic” model, capable of generating the channel structures, and a basic geostatistical model, both applied to the same data from an alluvial aquifer. An equivalent homogeneous medium is also considered, for reference. In the conclusion, new geostatistical methods will also be discussed.

The comparison between the genesis approach and the other techniques applied to the Aube River floodplain was partly presented in Teles et al. (2004). The comparison focused only on three-dimensional structures, head, and velocity fields. The transport results were only briefly investigated. Here, we analyze in detail the transport simulation results.

To begin with, we briefly present the two methods: the genesis model and the sequential Gaussian indicator method. Next, we discuss their application to an alluvial aquifer in the Aube Valley (France) and present the main conclusions on the generated geological structures obtained by both methods; groundwater flow simulations and results are then summarized. Results of solute transport simulations in the numerical aquifers created by the different models are analyzed and compared.

HETEROGENEITY CHARACTERIZATION METHODS USED IN THIS STUDY

The Genesis Model

The idea underlying genesis methods is that subsurface heterogeneity is primarily due to the variability of geological processes during the creation of the geological strata and thus that the creation mechanisms should be considered in order to characterize the resulting heterogeneity.

At geological time and space scales, it is not practically feasible to solve hydrodynamic and sediment transport equations that can be used for engineering purposes at smaller scales. A few genesis models exist: Koltermann and Gorelick (1992, 1996) simulated the sedimentation in the San Francisco bay over 600,000 yr, using a model based on Tetzlaff and Harbaugh's (1989) approach with simplifications of flow, sediment load transport, and deposition equations; Paola et al. (1992) on the one hand, and Grandjeon (1996), Grandjeon and Joseph (1999), and Quiquerez et al. (2000) on the other, used a diffusion equation to model two-dimensional (2-D) alluvial and three-dimensional (3-D) marine sedimentary basin creation, respectively; Webb and Anderson (1996) coupled a random network generation with hydraulic equations to generate the geometry of channels and bars of braided river systems; Howard and Knutson (1984), Howard (1992), Sun et al. (2001a, 2001b), and Lancaster and Bras (2002) developed models of meandering channel evolution, but they did not include a stratigraphic description of the alluvial floodplain deposits. Ritzi et al. (1995) developed hydrofacies distribution in valley aquifers, but used a statistical approach based on genesis concepts by defining proportions and transient probabilities. Fogg et al. (1998), Weissmann and Fogg (1999), and Weissmann et al. (1999, 2002) also used the Markov chain approach to represent alluvial aquifers.

A few authors have studied the implications on flow and transport in modeled aquifers with genesis approaches, and, in particular, the importance of connectivity in solute transport prediction. Webb and Anderson (1996) concluded that their genesis model of braided deposits led to enhanced preferential transport but did not compare their simulations with geostatistical ones. Scheibe and Freyberg (1995) developed a geometric model of alluvial aquifers that simulates various bars shapes. Later on, Scheibe and Cole (1994) and Scheibe and Murray (1998) compared transport in this modeled alluvial aquifer to that in several alternative geostatistical representations. They concluded that connectedness of extreme-valued parameters are critical features for transport in heterogeneous porous media. Western et al. (2001) also pointed out with several examples the importance of connectivity in hydrologic processes. From the geostatistic community, Journel and Deutsch (1993) compared flow and transport behavior in geostatistically simulated aquifers to that in a digitized outcrop image.

The genesis model in this study (Teles, 1999; Teles et al., 2001, 2004) was developed to mimic the genesis of alluvial media and to obtain a three-dimensional heterogeneity characterization. It uses simple geometric rules to model several types of fluvial systems: braided, meandering, and channel incision. The model does not simulate the hydrodynamics of sedimentation. Instead, the main sedimentary processes are modeled, either at the stream scale or at the plain scale, by geometric sedimentary rules chosen from the literature on fluvial geomorphology or as assumptions based on observations. The model is fully described in Teles et al. (2001, 2004), to which the reader is referred for details; only a short description is provided here.

The computational domain consists of horizontal square meshes, on or from which “sediment” entities may be deposited or eroded during a simulation. Sedimentary structures such as bars are defined by several of these “entities.” Altogether, the “sediment” entities form the modeled alluvial medium. The model generates sedimentary structures through time, based on an inferred sequence of genesis periods. Genesis periods represent time lags in which the modeled system exhibits a given behavior such as braided net aggradation, channel incision, general silting of the domain, etc. This method needs an expert interpretation of actual observations, in terms of geological sequences of deposits, their horizontal extent, and associated sedimentary processes. As natural records might be truncated, since deposits may have disappeared by erosion, some climatic periods might not be represented by the genesis periods if they are based only on lithofacies observations. In such complex systems and without dating, it is sometimes impossible to determine whether two lithofacies result from contemporaneous processes. By definition, then, a climatic period might correspond to several genesis periods. Adding more knowledge (e.g., by dating) would be beneficial to the genesis approach because it could lead to features that are not recognized on the studied cross section but that exist elsewhere in the floodplain. The discrepancy between genesis and climatic periods is, however, not really important because the objective is to reproduce the structures in the cross sections and not necessarily follow their exact genesis sequence.

Figure 1 shows a flow chart of a genesis-period simulation. The major active geomorphological process is defined for each genesis period according to actual observations. It determines both the kind of model (meanders, braided bars…) to consider and the spatial scale and action area of the process. Depending on the modeled system, a time step is arbitrarily chosen. Braided-system simulation is divided into several time steps because the entanglement of several bars is the key to a correct simulation of the geometry of such deposits. In contrast, when the deposit is homogeneous over the whole domain, or in a channel, or during incision, there is no need to divide the simulation into several time steps, and only one time step is considered.

For all types of environment, the model can either erode or deposit. If the period is an accumulation period, the type(s) of entities to deposit are chosen according to lithological information specific to the modeled area.

Braided systems are simulated as random processes. Considering various classifications (Miall, 1985, 1992; Collinson, 1986; Fujita, 1989; Boggs, 1995), three schematic bar shapes were chosen to represent the braided bars in the model (longitudinal bars in the stream direction, transverse bars that grow across the stream direction, and lateral bars attached to the banks.) These bars are the basic elements in the braided-system modeling. The model randomly simulates these bars within an active zone. Once the bar shape is defined, sedimentary entities are deposited on the modeled floodplain in such a way that deposited entities cover the bar shape. When the floodplain is rising or being incised, “sediment” entities are taken off around the bar to simulate the incision of surrounding channels. A braided-genesis period is simulated by a succession of several bar generations covering the whole active domain. When the bar lithology is homogeneous, it leads to a homogeneous formation with geometric elongated features in the valley axis.

Considering the spatial extent of floodplains, meandering systems are only crudely represented as a succession of sinusoidal channels from upstream to downstream within the active zone of the alluvial plain. Period, amplitude, and length of the channels are drawn at random according to the plain size. If observed data are available, the model can be constrained to simulate channels passing at some given locations.

The simulation of a genesis period is validated quantitatively against available data, such as the thickness of the sediments deposited during a given period, or the elevation of sediments, terraces, and paleochannels. This genesis model does not simulate physically based processes but rather reproduces the sedimentary geometries. The sequence of genesis periods is important because a given period simulation may affect the next ones by defining the topography. The model can thus reproduce the interlocking pattern of different alluvial units.

The genesis model uses the geologist's knowledge to its maximum. Thus, it is much more deterministic than the geostatistical method, especially concerning lithofacies, which are prescribed. That is the reason why a single simulation was run with this approach. Its variability was considered as negligible compared with that of the geostatistical realizations.

The Sequential Indicator Simulations

A classical (simple) geostatistical approach to characterize heterogeneities, the sequential indicator method (Gomez-Hernandez and Srivastava,1990;JournelandGomez-Hernandez, 1993; Gomez-Hernandez and Cassiraga, 1994; Marsily et al., 1998, 2005), was used to simulate the same area with the same initial observed data set. This indicator approach, described in Journel (1989) and Deutsch and Journel (1992), divides a variable into a set of classes and computes the probability of the parameter falling into each class. This method allows the simulation of multimodal or highly skewed distributions (Marsily et al., 1998). The approach used here is an extension of this method applied to discrete variables. It was presented in detail in Teles et al. (2004) and is only summarized here; it uses a code that was specifically developed by one of us (F. Delay). Since sequential simulations are intended for Gaussian fields, a normal score transform (Deutsch and Journel, 1992) must be performed on the data to obtain a standard Gaussian distribution with zero mean and unit variance. For each layer, a geostatistical analysis is made to infer variograms for each lithofacies. The simulations give the probability of occurrence of each lithofacies throughout the layer. They are carried out with an anisotropy factor of ten for the variogram in the main downstream direction on a grid with a resolution of 100 m × 100 m. Ten realizations are generated for each level and each available variogram. Simulations are back-transformed to the initial distribution with an inverse normal score transform. The layers are constructed one by one according to the resulting probability of the lithofacies at each node of the grid and the number of lithologies supposed to appear at that level. The three-dimensional medium is then reconstructed by combining all the simulated layers.

The use of this relatively simple geostatistical tool may be seen as “not giving its chance” to the geostatistical technique when compared with the genesis model. Although parts of this criticism may be valid, and will be discussed again later, we wish to emphasize that, given the very limited source of data, it was not possible to infer real variograms in the longitudinal direction, nor to develop a full stochastic simulation using, for example, multi-Gaussian truncated approaches (Matheron et al., 1987; Galli et al., 1994; Chiles and Delfiner, 1999; Armstrong et al., 2003), or a multiple point geostatistical approach (Strebelle, 2002) or Boolean models. Additional reconnaissance and inference of the structure, shape, and geometry of the facies in the medium would have been required to improve the geostatistical model, as will be emphasized in the conclusion. We wanted to use a very basic geostatistical method in the simplest way. Thus, we chose the sequential indicator simulator because it is the most currently used method in the oil industry to simulate heterogeneous media.

NUMERICAL REPRESENTATION OF THE AUBE RIVER

Site and Data Set Description

The alluvial aquifer considered here is a 5.4-km-long by 1.4-km-wide area of the Aube River, France (Fig. 2). A cross section was constructed from a series of 44 auger-hole logs and a geomorphologic description along a road crossing the floodplain (P. Antoine and J.F. Pastre, 1999, personal commun.). Twenty-three lithofacies were observed. At the bottom of the valley, on the top of the chalk bedrock, Quaternary gravelly fluvioglacial sediments were deposited probably in a braided periglacial environment during the last glaciation period of the Pleistocene. Several paleochannels incised the gravel during the Holocene (10,000 yr B.P. to present). These channels have different widths, depths, and in-fillings: sand, silt, or peat sediments. Finer loamy deposits were also identified on top of coarse ridges, as well as silted-up channels. These deposits were probably deposited during widespread floods that affected the entire floodplain.

The actual data, 44 auger-logs along a line crossing the floodplain, present a difficulty for the geostatistical inference, since they cannot provide spatial correlations in the floodplain horizontal axis direction, which is crucial for extrapolating the observations. In order to provide this information, three other cross sections were arbitrarily built from the actual one (Teles et al., 2004). Vertically, sediments are described by ten different layers representing 4-m-thick alluvia. The four transects of 44 logs with ten layers sums up to a total of 1760 conditioning points, which are considered as “true” data to be exactly reproduced in the stochastic geostatistical simulations. In other words, we tentatively and indirectly extended the correlation length of data in the main floodplain direction. We will see later that this is not enough to make geostatistics mimic channel structures.

Simulated Structures

For both simulation methods, lithofacies are translated into hydraulic conductivities. These were not measured at the site, but orders of magnitude were inferred from similar alluvial deposits. 01Table 1 shows the horizontal and vertical hydraulic conductivities assigned to each lithofacies. They span five orders of magnitude. The two lithofacies simulation methods have different vertical and horizontal resolutions. Hydraulic conductivities are homogenized by upscaling to the groundwater-flow simulation resolution. Our comparison of hydraulic conductivities was made after homogenization.

The simulated sedimentary structures are very different (Fig. 3, plan view, top level, and Fig. 4, cross sections). At the surface, layers 1–5, which represent the first meter of alluvia, are similar for all heterogeneous representations. They are almost homogeneous with a hydraulic conductivity of 5.10−8 to 10−7 m s−1. Layers 6 and 7 of the genesis simulation show a large continuous channel (Kh ∼ 3.10−4 m s−1, where Kh is the horizontal hydraulic conductivity) embedded in a low conductivity medium. For layers 8–13 of the genesis model, a higher-conductivity matrix (Kh ∼ 5.10−4 m s−1) corresponding to the gravel deposit surrounds several entangled channels with hydraulic conductivities ranging from 2.10−4 to 4.10−7 m s−1. To represent these layers (6–13), the sequential Gaussian simulations generate elongated structures in the flow direction but with no apparent continuity significant at the scale of the floodplain, at least in these two-dimensional slices (three-dimensional connectivity may exist in some cases, and can be characterized, see, e.g., Proce et al., 2004, but would occur here only randomly, as each slice of the geostatistical model is generated independently of the others). Below, down to layer 18, narrow channels are present in the genesis simulation, the geostatistical realizations are broadly homogeneous. The sequential Gaussian simulator was unable to find any correlation between the channel lithofacies that appear only punctually at some auger holes. Note that if we had artificially increased the correlation length of the variograms, the results would hardly have improved. There are not enough conditioning data points, even with the three fictitious transects added to the data set (see previous discussion), to allow the basic geostatistical method to generate continuous channels. The last two layers are homogeneous in both approaches. The vertical continuity of “genetic” channels is visible on the cross sections of Figure 4, whereas the geostatistical realization exhibits scattered channels as each layer is simulated independently from the others.

GROUNDWATER FLOW SIMULATIONS AND RESULTS

Hydraulic Conductivities and Flow Conditions

In order to assess the value of the genesis approach in characterizing heterogeneity as compared to the geostatistical approach, groundwater-flow and transport simulations were conducted on three different representations of the alluvial aquifer: the genesis model reconstruction, the sequential Gaussian simulations, and an equivalent homogeneous representation. The homogeneous representation was estimated by upscaling the hydraulic conductivities of the genesis model using the simplified renormalization method proposed by Le Loc'h (1987) and also described in Renard and Marsily (1996). The horizontal and vertical upscaled hydraulic conductivities were Kh = 5.10−4 m s−1 and Kv = 10−6 m s−1, respectively. Spatial resolutions for these simulations were set to horizontal 100 × 100 m square meshes and 20-cm-thick layers.

Groundwater-flow simulations were made with the Visual MODFLOW commercial software (McDonald and Harbaugh, 1988). The aquifer was considered to be confined in order to always have a fully saturated domain including all the layers. Constant-head boundary conditions were applied to the top layer at each end of the domain, upstream and downstream. Noflow conditions were defined on all other sides. The horizontal head gradient imposed on the top layer was 5 m over 5200 m (9.62 × 10−4).

Summary of Flow Simulation Results

Head and discharge results, described in Teles et al. (2004), show that the geostatistical aquifers are more permeable overall than the genesis simulations. The homogeneous equivalent aquifer computed by upscaling the “genetic” aquifer should theoretically have a hydraulic conductivity yielding a flux equivalent to that of the genesis model for the same boundary conditions. With the same prescribed head difference, the head map is similar, but surprisingly, the homogeneous case has a higher discharge, similar to that of the geostatistical simulations. This suggests that the heads are strongly constrained by the domain size and by its boundaries and that the renormalization method might not be very accurate for such heterogeneous structures with high permeability contrasts in channel-like features poorly represented by a covariance function. Nevertheless, one reason why the geostatistical and homogeneous models have almost similar discharges is that they represent a semicontinuous medium of equivalent bulk conductivity. On the contrary, the strong channeling (i.e., semidiscontinuous medium) stemming from the genesis model has a much weaker discharge because flow is effective only through a limited portion of the aquifer, i.e., the high-conductivity gravel at the bottom of the medium.

SIMULATED TRANSPORT RESULTS

Transport simulations were done on each aquifer simulation, for the same steady-state flow conditions as described already. An initial concentration of 106 mg L−1 of a conservative solute was prescribed at the upstream end, in each cell and each layer of the simulated aquifers, as a quasi-Dirac delta function (one time step of approximately one week followed by a natural evolution without any additional injection for ∼100 yr). Density effects were neglected since the concentration unit is arbitrary; it is only the ratio between the input concentration and the measured one that matters. Neither local dispersion nor sorption was considered, a perfect tracer was assumed. Including local dispersion in the model could have an effect on large-scale transport by enhancing exchanges between layers; we did not do it in this comparison, because it would also reduce the differences between the two models. The MT3D code was used to simulate the advective transport alone. Transport parameters were homogeneous in the whole aquifer and for all simulations. The effective porosity was set at 0.15. It is clear that the porosity may be different for each facies. This effect could easily be included, if porosity measurements were available, but it was not considered in the present simulations for simplicity, since the effects are likely to be of second-order.

Each layer had a different behavior as regards contaminant transport (Fig. 3). At the surface, hydraulic conductivities were low overall. The concentration fields exhibited patterns that were strongly linked to the propagation of the solute in the deeper, more conductive layers. There was almost no horizontal propagation of the solute in the first surface layers. The solute flowed down into the deeper, more conductive layers. The concentration visible at the surface depended on occasional more conductive “windows” in the surface layers through which the solute could propagate upward. At the bottom, homogeneous high-conductivity layers had a similar pattern from one representation to the other. Figure 5 shows cross sections of hydraulic conductivities and associated concentration fields. In the “genetic” medium, the contaminant flowed first in the deeper layers (left panels) and later arrived in the less permeable top of the simulated aquifer. In the geostatistical realizations, the contaminant plume had a greater thickness and appeared more scattered than in the genesis simulation. The end-tail of the plume flowed through the same channels throughout time. In the genesis model, the lithological change between the bottom and the top of the alluvial aquifer was due to a change in deposition processes and to a finer lithology toward the end of the simulation. This led to two well-differentiated zones. Thus, the “genetic” flow was constrained effectively through the deeper portion of the aquifer, as suggested by the significantly lower discharge (Teles et al., 2004). The simple geostatistical realizations were unable to simulate this vertical geological sequence, although more sophisticated geostatistics might be able to do it. Simple geostatistical realizations generate noncontinuous low-conductivity patterns, imbedded in a high-conductivity matrix, and the contaminant can therefore flow around these patterns and spread more easily. For both approaches, the represented channels and zones of low permeability affected the concentration field, but in a different way. For example, the impervious channel on the left (marked A in Fig. 5) of the “genetic” cross section clearly acted as a hydraulic barrier to the contaminant, and this did not occur in the geostatistical simulation.

These transport simulations show the importance of the generated structures for the spread of pollution. The added value of the genesis method is to generate (by construction) continuous features that have geometry and location estimated according to actual observations. These patterns can greatly influence the transport behavior of the aquifer, creating either preferential paths or hydraulic barriers.

In an attempt to quantitatively characterize the transport, breakthrough curves were extracted from these transport simulations. The idea here was to compare the simulation results as a hydrogeologist would when studying an aquifer in the field. We assumed that the heterogeneity was unknown and that transport was simulated in an equivalent homogeneous medium, where the heterogeneity was described only by the equivalent dispersion coefficient that best represented the observed spread of the plume. Velocity was thus assumed constant over the whole domain (both vertically and horizontally), and velocity variations due to the alluvial structure (channels, etc.) were assumed to be represented by the equivalent dispersion coefficient. The spreading and dispersion was only studied in the horizontal direction. Therefore, we calculated average concentrations on a vertical plane crossing the valley, at different distances from the source and at different times. These concentrations are the arithmetic averages in the meshes of the plane. For the two heterogeneous cases, one could consider averages weighted by the Darcy velocity in the x (downgradient) direction to obtain the solute flux. But we assumed here that the medium was homogeneous, and that the velocity variations were unknown and only represented by the equivalent dispersion coefficient. Thus, these averages were those that would be measured in situ, if many local concentration measurements were taken on each vertical plane, without any knowledge of the local velocity distribution in space, which is generally unknown.

Breakthrough curves in time (Fig. 6) show the average concentration for each given 2-D transect at different distances from the source. For the genetic, homogeneous and three geostatistical representations, breakthrough curves were computed for every possible vertical transect, i.e., at 100 m resolution. Curves for transects close to the source were characterized by a dominant vertical redistribution of the contaminant and nonhomogeneous downward flow. Subsequently, these curves were eliminated for simplicity. Figure 6 shows the curves for transects farther downstream from the zone where the main movement was vertical. The homogeneous case exhibits classical breakthrough curves that widen as the distance from the source increases, and which have a maximum that decreases as the contaminant flows downstream. This behavior is due to the longitudinal dispersion of the plume with time, induced by the boundary conditions imposed only on the surface layer, thus creating an initially downward flow, and different travel lengths and velocities along the flow paths; numerical dispersion may also play a small role, since pure advective transport was simulated. Breakthrough curves for all the heterogeneous cases exhibit long tails with non-negligible concentrations. These late arrivals of solute are due to less-permeable heterogeneities that delay the migration of the contaminant. The widening of the curve and the decrease of the maximum concentration are faster for the “genetic” curves. The “geostatistical” curves are more difficult to interpret. They show nonhomogeneous maximum decrease and curve widening that may be related to the more heterogeneous and nonconnective appearance of the geostatistical simulations.

For each curve, the mean transit time and apparent velocity were computed as if the medium were homogeneous (i.e., no information was available on the heterogeneity, and the breakthrough curves were interpreted for an equivalent homogeneous medium). They are defined as follows:  
formula
and  
formula
where tm is the mean transit time, x is the distance from the source, C(x,τ) is the concentration at distance x and time τ, τ is the time integration variable, and 〈u(x)〉 is the mean apparent velocity.
For each transect, variances were defined by:  
formula
Apparent dispersion coefficients (Da) in the time domain were estimated from these curves, for each distance, with the following definition:  
formula

Although the variances and apparent dispersion coefficients calculated in this way would only correctly represent a truly Fickian transport, which is not an accurate description of the complex mechanisms taking place in these heterogeneous fields, we still used them for comparison between the two models, as a way to characterize the heterogeneity, as one would do in the field by a simple interpretation of the observations with the classical model.

Alternatively, the effective dispersion coefficient (De) was defined as the derivative of the variance:  
formula

These two definitions should converge for statistically homogeneous media. Variances, apparent and effective dispersions are shown in Figure 7. They show very different behaviors for each type of model. For the homogeneous case, variances and dispersions quickly reach a constant value. For the genetic medium, computed variances decrease with distance from the source. The same trend is visible for the geostatistical realizations, although there is a higher variability from one point to the other. As part of the contaminant is trapped upstream by low-permeability structures (long tails on Fig. 6), variances are computed in truncated breakthrough curves or decreasing contaminant volume in the downstream direction. This leads to a decrease in the variance with distance from the source. This decrease of variance is clearly an artifact due to boundary conditions, which strongly influence the transport behavior, particularly in the downstream part. Similarly, for the homogeneous representation, the apparent dispersion coefficient is almost stable around 1.10−4 m2 s−1; the effective dispersion has small negative values, no lower than −1.2 × 10−3 m2 s−1. For the genetic medium, the apparent dispersion increases to a maximum of around 2.5 × 10−4 m2 s−1 then decreases steadily to 1 × 10−4 m2 s−1. The genetic effective dispersion is almost constant at the beginning then starts to oscillate over an order of magnitude around 10−3 m2 s−1 (between 2000 and 5000 m on the x axis, i.e., for the downstream part). Geostatistical realizations have highly variable apparent dispersion coefficients, larger by one order of magnitude. Their effective dispersion also shows a higher variability. They reach values that are one order of magnitude larger (around 5 × 10−2 m2 s−1). This high variability is not significant, we believe, and is due to insufficient calculation points; only the general trend is of interest. The decrease of the apparent dispersion, computed as a linear function of the variance, is due to the decrease in the variance because of truncated breakthrough curves. The effective dispersion is computed from the variance derivative and thus is less sensitive to the decrease in variance but does show higher and more variable values for the last part of the curve (between 2000 and 5000 m on the x axis). For both heterogeneous approaches, the effective dispersion coefficient is higher than the apparent dispersion by one order of magnitude

Dispersion from spreading in space can also be computed, at different times. Variance in the x-axis direction (σx2), or contaminant longitudinal spreading, was computed at every available time step for each generated medium. The effective dispersion in space (Ds) is then defined as the variation with time of the variance in the x-axis direction:  
formula

All heterogeneous representations exhibit increasing longitudinal variances as time increases and reach a similar threshold value (Fig. 8). Variances calculated for the geostatistical realizations increase faster than the one for the genetic representation. As a consequence, longitudinal dispersion coefficients of geostatistical realizations reach their maximum (around 2 × 10−3 to 2.5 × 10−3 m2 s−1) faster than the genetic coefficients, the maximum of which stays under 2 × 10−3 m2 s−1(Fig. 8). These results are consistent with the fact that the sequential Gaussian simulations have a more heterogeneous and nonconnective appearance, which leads to a higher local variance of the velocity field (Teles et al., 2004) and, thus, a higher value of the dispersion coefficient. As expected, the equivalent homogeneous representation exhibits variances and dispersion smaller by one order of magnitude. The decrease of both variance and dispersion for late times is due to a loss of mass at the downgradient outlet of the system.

For the homogeneous and geostatistical cases, results from both definitions (the apparent dispersion from time distributions and the spatial spreading dispersion) are similar in orders of magnitude but not in behavior. For the heterogeneous media, effective dispersion coefficients are higher than the apparent ones. Theoretically, definitions (equations 4 and 5) of dispersion are consistent with the Gaussian assumption for the permeability and velocity fields (Dagan, 1989, 2002; Gelhar, 1993), which may explain these similarities. This assumption is not verified for the simulated media, which may be the reason why the different definitions do not give similar values. All methods show higher dispersion coefficients for the geostatistical realizations, which have more scattered noncontinuous heterogeneous structures.

These results suggest that the different dispersion definitions are not equivalent for heterogeneous non-Gaussian fields. In fact, dispersion coefficients defined in the time domain include a strong assumption on the velocity field: a single mean velocity U(x)> is used in the definition of the equivalent homogeneous medium. The dispersion definition in space does not contain this assumption and computes the dispersion coefficient directly from the observed spreading due to velocity variations. It should be preferred in practical field applications.

DISCUSSION

Results of groundwater flow and transport simulations conducted in this study show that flow velocities and contaminant plumes are significantly different in the different representations, depending on which three-dimensional structure is generated, although groundwater simulations are strongly influenced by the prescribed boundary conditions. In this example, both the hydraulic conductivity distribution and connectivity are different in geostatistical realizations and the genesis simulation. Zinn and Harvey (2003) showed that, even for nearly identical log-normal distributions of hydraulic conductivity, connectivity of high-conductivity or low-conductivity structures rather than connected structures of intermediate values (in the multivariate log-Gaussian case) can lead to flow and transport behaviors that are not predicted by the stochastic theory. In the case of a Gaussian stationary geostatistical field, defined by a mean and a covariance, the temporal and spatial definitions of the apparent dispersion should coincide, if the domain is infinite and if the travel distance of the solute is long enough to reach ergodic behavior (in general, approximately ten times the correlation length; Gelhar, 1993). In the geostatistically generated medium, these two definitions do not coincide, but are different by a factor of only 2.5 (approximately). This (small) discrepancy is due to the effect of conditioning the medium and to the finite size of the domain; the maximum travel distance (5200 m) is indeed just about ten times the correlation length (500 m for the largest indicator variogram range; Teles et al., 2004), but only for late times. For the genetic medium, the two definitions differ by a factor of ∼20. This larger discrepancy can be explained as follows: the migration of solutes in a continuous channel in the main direction of flow imbedded in a low-permeability medium can be seen as equivalent to the migration in a perfectly stratified medium, with no or little exchange between strata; it is well known that in such a case, the dispersion increases linearly with time, and not with the square root of time; no ergodic limit can be reached, and the spatial and temporal definitions of dispersion will never converge. The same is true if the continuous channels are of low permeability and act as barriers to transport. Therefore, the very strong continuity of the lithologic facies created by the genetic method explains the difference between the spatial and temporal definitions of apparent dispersion. More simply, suppose that the major part of the solute moves within a single fast channel, but that a minor part remains close to the source, trapped in a low-permeability region. At a given distance from the source, the breakthrough curve only sees what happens in the fast channel and displays the relatively small apparent dispersion within that fast channel. On the contrary, the spatial spreading of the plume is then very large, because it also considers the minor part of the plume that remains close to the source. Thus, a dispersion calculated with time distributions is small, whereas the correlation length of Lagrangian velocity fields predicts a large dispersion calculated on spatial particle spreading (see definition of dispersion in the Lagrangian context in, e.g., Dagan, 1989).

These results of the comparison between apparent dispersion coefficients underline the fact that field methods commonly used in hydrogeology rely heavily on strong assumptions on the nature of subsurface three-dimensional structures. Most of the time, field hydrogeologists do not know the aquifer heterogeneity, neither do they have access to the apparent dispersion coefficient defined in space by monitoring the contaminant spatial spreading, which would require a very extensive observation network. The common method for estimating dispersion coefficients is therefore to make tracer tests and to determine the apparent dispersion coefficient at fixed distances by the widening with time of the breakthrough curves. This definition of the dispersion coefficient embeds the assumption that a single mean velocity value can be used in modeling heterogeneous fields. For Gaussian fields, this assumption may be valid if the travel distance is large enough so that the ergodic limit has been reached; the dispersion coefficient estimated in the time domain can then be used in space. But this approach, as shown here, will give wrong values of the actual dispersion in space in highly heterogeneous media with strong channeling or permeability differences. One could argue that, when channel flow occurs, as shown for instance with the genesis model, a classical convection-dispersion equation can no longer predict solute transport, and that it is unnecessary to identify an equivalent apparent dispersion coefficient (or a tensor). This would be true when accurate facies simulations are available, with well-constrained scenarios (well-known boundary conditions, injection conditions, and acting transport mechanisms). But in practice, an overall equivalent description of the system is much more useful and relevant in view of the uncertainty attached to facies-constraining data. Upscaled models that can be conditioned and inverted even on scarce data are, we believe, still useful. Extracting a “time-domain” apparent dispersion coefficient from available field data and transforming it into a true “spatial-domain” apparent dispersion is not straightforward. To our best knowledge, neither a sound theoretical basis nor empirical recipes exist for this, at least for highly contrasted velocity fields. The transport simulations reported here emphasize the need for further work on this issue, both for porous media and for fractured rocks.

CONCLUSION

Geostatistical methods consider subsurface formations as mathematical realizations of random fields, the spatial correlation parameters of which are statistically inferred from observations. Very often, a multi-Gaussian distribution of the random field is assumed, but non-Gaussian distributions can also be considered. One of their advantages is that they can be fully conditioned to any kind of real information, i.e., they will reproduce observations at each observation point. They are very useful and objective tools, especially because they require few pragmatic or subjective choices from the modeler. This approach has been found to be valid for many geological formations, but, as we have shown here, and for the type of geostatistical tool that we used, they are poorly suited to describe continuous “channel-like” structures, as found, for example, in alluvial deposits.

Genesis models, on the contrary, have the opposite philosophy. They rely strongly on pragmatic choices concerning formation processes and their succession. They are able to consider general characteristics that cannot be inferred from sparse local data but are associated with formation processes that can be studied and observed elsewhere. This is the case, for instance, of the nature of sedimentation processes, such as meandering, incising, and braiding, which are valid for a given region and a given time period, and not just locally or for the type of deposited sediments (gravel, sand, silt, etc.). They can more realistically reproduce sediment patterns, which geostatistical methods can hardly represent, such as continuous channels acting as flow conduits or barriers. One drawback of the genesis method is, however, that it is presently very difficult to condition them to actual local observations. Multipoints geostatistics (e.g., Strebelle, 2002; Strebelle et al., 2006) trained on 3-D images created by a genetic method may, however, be one option to condition the genetic methods and to reconcile the two approaches.

In the particular case of alluvial formations, such as the Aube valley, the actual data set composed of 44 auger-holes distributed over a cross section of the floodplain favors the genesis approach over the basic sequential Gaussian indicator. It makes it easier to determine the succession of periods (with their major acting geomorphological processes, their sediment types, etc.), but, conversely, this data set does not provide the crucial information on all three directions, which is required by the geostatistical method. Better-distributed data would have improved the strength of the geostatistical approach. However, even with the same number of observation points evenly distributed over the domain, the geostatistical approach would not have been able to generate features such as channels or take into account lithofacies scarcely represented here, such as peat. With a more even distribution of observed data, the genesis method probably would have needed additional information such as dating by 14C of organic debris, commonly carried out in geomorphic studies, in order to better infer the succession of genesis periods.

Today, upscaling, renormalization, and other techniques to determine whether or not a numerical reservoir correctly mimics a macroscopic behavior addressed by field data are based on flow simulations. One of the main problems raised by this paper is that, for flow problems, almost any upscaling technique works, with small discrepancies. For transport problems, this is quite another story. This paper shows how a standard tool that succeeds and passes tests for flow, may fail for transport.

We have shown here that a correct representation of transport processes in an alluvial environment can hardly be obtained with simple geostatistical tools, and that genesis methods can better handle such situations. Work remains to be done to compare the genesis approach with other geostatistical techniques, such as the Markov chains, the multiple-points, and the Gaussian threshold, and also to develop methods to condition the genesis method to actual local observations. However, improving methods to accurately represent and condition, for example, alluvial reservoirs, as done by the genesis model, produces other difficulties. We have shown, for instance, that transport within channeled flow fields over relatively large distances does not lead to a macroscopic behavior as regards solute dispersion. The paradigm in transport studies is that, on the one hand, accuracy is sought to better understand the physics of transport processes at the fine scale and lower the risk of neglecting an important mechanism, and, on the other hand, these physics must be simplified at the large scale for pragmatic purposes. Resolving these apparently contradictory issues is the task of modern approaches to reservoir modeling.

We received very strong criticisms from colleagues for our view that the genetic approach is better than indicator simulation. Our conclusion that geostatistical tools are in general “poorly suited to describe continuous channel-like structures” was considered not true. We supposedly did not apply indicator simulation properly and used an outdated approach in indicator simulation, whereas newer, more appropriate forms of indicator simulation could have simulated channels in fluvial and alluvial deposits. Such simulations correctly show a high degree of channel interconnectivity. We were told that to effectively simulate unit types with indicator geostatistics, all of the following steps should have been followed: (1) study how the proportions of unit types vary over space from the borehole logs, and create zones in which it makes sense to assign proportions to the unit types, as different from other zones; (2) characterize the geometries of the unit types, and infer lateral correlation scales based upon the expected geometry, not from borehole statistics; inferences from borehole data can be spurious, and two-point bivariate correlation models can be improved from an understanding of the geology; (3) simulate the unit types in the zones separately, and merge the simulated zones together; an indicator simulation method exists which, through transformation of coordinates, can simulate curvilinear structures; this has been used for folded strata, and can probably be applied to curved fluvial channels; and (4) verify that the simulations have the proper proportions, length statistics, and transitioning in each zone; this must be checked before a realization is accepted; the important point is that general knowledge of the proportions and geometry of the unit types can and should be used to choose the indicator correlation structure, even when boreholes are not sufficient for defining lateral correlation; references to the work of Johnson (1995), Carle (1996, 1998), Weissmann and Fogg (1999), Weissmann et al. (1999), Ritzi (2000), Ritzi et al. (2000), Proce et al. (2004), Dai et al. (2005) and Harter (2005) were provided to support these comments.

We consider this criticism as a valuable contribution to a discussion of this paper, which we wish to publish with it, but not as a direct conclusion of our work, nor something we fully agree with since we have not tested it. We maintain that our conclusions are valid for the simple and classical geostatistical method that we have used. By comparing a simple genetic approach to a simple geostatistical method, we wanted to advocate that we absolutely need to integrate geological knowledge of the natural systems and processes in order to be able to correctly describe them. We believe most, if not all, of our colleagues would agree with that, even if the approach we used differs from that suggested in these comments.

*Corresponding author gdm@ccr.jussieu.fr

We thank the reviewers and the editor of Geosphere for their comments and help, and the editor for an interesting discussion and conclusion on the balance between criticism, editorial decisions, and freedom to express one's own views.