Conditional Simulation of Hydrofacies Architecture: A Transition Probability Approach

Abstract

A practical and now-popular application of geostatistics in hydrogeology is to use conditional simulation to assess impacts of subsurface heterogeneity on flow and transport processes. A transition probability approach is amenable to integration of geological concepts into analysis and modeling of spatial variability and conditional simulation of hydrofacies architecture. Hydrofacies may be defined from geological perspectives (e.g., a fluvial system with channel, splay, levee, or flood plain deposits) or textural descriptors (e.g., mud, silt, sand, gravel,…). As an alternative to variogram or covariance-based geostatistical modeling approaches, categorical geostatistical methods can be implemented with the geologically interpretable transition probability. The Markov chain model, as applied in a spatial context, is amenable to integration of geological ways of thinking as an alternative to empirical reliance on exhaustive or “training” data sets to implement geostatistical conditional simulation algorithms. Applications of conditional simulation to hydrogeological problems will benefit from using geological approaches to separate out the more deterministic and stochastic components of subsurface characterization. The overarching goal is to reduce uncertainty in flow and transport analysis and interpretation of hydrogeological data.

1. Introduction

1.1 Stochastic aspects of hydrofacies architecture

Subsurface characterization – the foundation of many practical hydrogeological studies – has both deterministic and statistical or “stochastic” aspects. A familiar stochastic problem in hydrogeology is “borehole correlation,” wherein the spatial distribution of the hydrogeological categories of subsurface materials of interest (e.g., “sand and mud” or “aquifers and aquitards”) is predicted to occur between boreholes. We refer to hydrogeological categories relevant to subsurface hydrological interpretation as hydrofacies and the three-dimensional (3-D) spatial distribution of the hydrofacies as hydrofacies architecture.

As illustrated in Figure 1, the “reality” of the borehole correlation problem is that rarely if ever are enough data available to characterize the true hydrofacies architecture. Two challenges develop:

1. In the practice of manually drawing cross-sections by geologic interpretation, there can be a tendency to overestimate the lateral extent of sand (or aquifer) to the point that the sands (or aquifers) are unrealistically layered and laterally connected.

2. In many hydrogeological applications, a realistic 3-D model is needed – one that will be useful yet not burdensome to produce for use in groundwater flow and transport analysis.

Figure 1.

Illustration of the borehole correlation problem for sand and mud (clay and silt undifferentiated) hydrofacies.

Beginning in the mid 1990s, “indicator” geostatistical approaches [1, 2] became recognized as an automated means for addressing the stochastic aspect of characterization of hydrofacies architecture within alluvial or fluvial depositional systems [3, 4, 5, 6]. However, we found that the indicator approach was difficult to apply in a practical 3-D hydrogeological setting because of its data-intensive and empirical approach to the geostatistical model development. At that time, we were looking for a method that would enable development of the geostatistical model by either empirical or conceptual approaches. In particular, we wanted to be able to relate the parameters of the geostatistical models to observable attributes pertinent to development of realistic 3-D representations of hydrofacies architecture.

1.2 Translating geological concepts into a stochastic model

A 2- or 3-D picture or “training image” of a heterogeneity pattern may be used to conceptualize, formulate, and parameterize a stochastic model. Here, we will focus on how to translate a picture of hydrofacies architecture (whether on paper or in one’s mind) into parameters of stochastic model. Figure 2 is a 3-D illustration of a fluvial environment which could serve as a realistic geological model of hydrofacies architecture. The illustration conveys the spatial attributes of the four hydrofacies – splay sands, channel sands, peat, and interchannel silt/clay – including their proportions, length scales, and juxtapositional tendencies (how one hydrofacies tends or tends not to occur adjacent to another). We will emphasize two main points here and later:

• hydrofacies definitions are best defined in a geological context, which does not necessarily directly correspond to texture (e.g., particle size descriptors), and

• proportions, mean lengths, and juxtapositional tendencies are fundamental observable attributes associated with spatial characteristics of hydrofacies architecture that can be directly integrated into the parameterization of a geostatistical model.

Figure 2.

A 3-D illustration of a mixed-load fluvial depositional system (modified from W.E. Galloway, personal communication; see also [7]).

Subsurface heterogeneity exists at multiple scales including smaller scales within the hydrofacies described above. Figure 3 provides an example of three scales of heterogeneity relevant to channel and splay sands hydrofacies. In this case, the hydrofacies correspond to the “megascopic” scale of heterogeneity. Depending on the hydrogeological problem, the multiple scales of heterogeneity may or may not need be addressed in either a stochastic or deterministic manner. In many (if not most) site-scale hydrogeological problems, the megascopic scale heterogeneity is recognized as a major uncertainty pertinent to analysis and prediction of flow and transport observations at wells because so many processes are strongly influenced by geometry and connectivity of the major aquifer and non-aquifer hydrofacies [8, 9, 10, 11].

Figure 3.

Heterogeneity at three scales relevant to fluvial hydrofacies complexes (modified from W.E. Galloway, personal communication; see also [7]).

2. Transition probability approach to conditional simulation

2.1 Transition probability

With an intention of putting more “geo” into geostatistics, we will outline the transition probability approach to geostatistical modeling in comparison to the preceding variogram-based approach. In the 1990s, approaches oriented toward mining geostatistics [12] were being extended to stochastic simulation [1, 2]. The mining-oriented approach relies on data-intensive empirical analysis of spatial variability using (cross-) variograms. Because hydrogeological systems are inherently 3-D and sparsely sampled in the lateral directions, these empirical approaches were not readily translatable to hydrogeological applications.

Spatial statistics such as the indicator (cross-)variogram and covariance are difficult to relate to the geologically interpretable observable attributes as applied to characterization of hydrofacies architecture:

• proportions – the volume fraction of a hydrofacies,

• mean length – the mean thickness or lateral extent of a hydrofacies, and

• juxtapositional tendencies – how frequent one hydrofacies tends to occur adjacent to another hydrofacies, including directional tendencies (e.g., upward-fining sequences; levees laterally adjacent to channels [13, 14, 15]).

Alternatively, the transition probability applied in a spatial context has a simple and interpretable definition in terms of a conditional probability

$$t_{jk}\left(h\right)=\mathrm{Pr}\left[k\text{occurs}atx+h\text{|}j\text{occurs}atx\right]$$E1

where j and k are discrete categories (e.g., hydrofacies), x is a spatial location, and h is a spatial separation vector or “lag” in geostatistical terms. The transition probability directly addresses a geologically interpretable question:

$$\text{If}j\text{occurs}atx,\text{what is the probability that}k\text{occurs}atx+h?$$

The transition probability statistic can also be applied in an embedded context, which pertains to discrete occurrences of the category in analysis of non-random juxtapositional tendencies such as fining upward cycles [16, 17, 18]. Of practical importance, the parameters of the embedded transition probability analysis are relatable to parameters of a spatial Markov chain model [13, 14, 15].

In geostatistical applications to categorical variables such as hydrofacies, mathematical formulations, analyses of spatial variation, and methods for stochastic simulation can be formulated with the transition probability instead of the indicator (cross-) variogram or covariance [13, 14, 15, 19, 20]. A spatial Markov chain model has been extended to 3-D applications and found to be a simple but effective stochastic model for hydrofacies architecture in alluvial and fluvial settings [13, 15, 20, 21, 22].

2.2 Conditional simulation

Conditional simulation refers to a stochastic algorithm designed to generate spatial distributions or realizations that honor both (1) known or “hard” data (e.g., hydrofacies identified at boreholes) and (2) the geostatistical models of spatial variability [1, 2]. The “tsim” code implements a transition probability approach to conditional simulation of facies architecture [19, 20]. In comparison, the “sisim” code, from which the tsim was largely developed, uses indicator variables and variogram-based modeling approaches [1, 2]. Here, we emphasize differences in transition probability and variogram approaches rather than the codes or numerical algorithms used to perform conditional simulation. Subsequent to tsim, other algorithms/codes have been developed for conditional simulation using Markov chain models [23, 24], lending further credence to a transition probability approach.

The original intended use of the transition probability approach and tsim was for conditional simulation of the hydrofacies architecture of an alluvial fan aquifer system underlying the Lawrence Livermore National Laboratory (LLNL) [13, 14, 15, 19, 20]. Hydrogeologists at LLNL defined four relevant hydrofacies – channel, levee/bar, debris flow, and floodplain – based on geological interpretation of hundreds of borehole logs and outcrops in roadcuts and excavation pits [11, 25]. The hydrofacies could be recognized in resistivity logs and textural assemblages (macroscopic heterogeneity in Figure 3) apparent in geological descriptions of core. There is not a direct correspondence between textural categories and hydrofacies in this application. Here, the transition probability approach is applied in the geometric context of megascopic-scale hydrofacies similar to the fluvial system in Figures 2 and 3, but adapted to the alluvial fan system of interest at LLNL, which includes a debris flow hydrofacies.

The remainder of this chapter begins with relatively simple example applications of the transition probability approach and finishes with more practical 3-D conditional simulation applications. The discussion emphasizes the conceptual side to geostatistical model development. In the 25+ years of our experience in applying and teaching geostatistical modeling, the concepts we discuss are essential to implementation of a transition probability approach to conditional simulation.

3. Example applications and approaches to stochastic simulation

3.1 Empirical approach using an exhaustive dataset

An exhaustive dataset provides an opportunity to verify the efficacy of a model or method. An exhaustive dataset included in [1] was used for verification testing of variogram approaches to stochastic simulation of lithofacies categories, the results of which showed inaccurate reproduction of the measured spatial variability [26]. For comparison, we will apply a transition probability approach to the same dataset.

The dataset shown in Figure 4(A) consists of four categories on a 50 $\times$ 50 grid with proportions of 0.4, 0.3, 0.2, and 0.1. It is readily apparent that the heterogeneity pattern is isotropic. Assuming the dataset represents “reality”, the objective here is to verify that a transition probability approach using a Markov chain model can characterize this pattern of spatial variability such that the “realizations” reproduce a pattern of heterogeneity similar to the “reality.” Visual inspection of stochastic simulation results shown in Figure 4(B) and (C) indicates that that the transition probability approach meets this objective and is extensible to 3-D as needed for flow and transport applications in hydrogeology.

Figure 4.

From empirical application of the transition probability approach: (A) exhaustive dataset and (B) stochastic simulation results for four 2-D realizations, and (C) a 3-D realization.

3.2 Modeling spatial variability

Development of the model of spatial variability is often the most challenging part of producing a stochastic simulation. In mining-oriented variogram approaches [1, 2, 7, 20], a two-step procedure is emphasized in modeling of spatial variability:

1. acquire data in a closely-spaced pattern to resolve spatial variability,

2. empirically fit a model of spatial variability to variogram measurements.

In hydrogeological applications, step 1 is often a roadblock to geostatistical analysis – there is rarely adequate data spacing to empirically characterize spatial variability of hydrofacies in the lateral directions. Hydrogeological problems may require more geological input, as outlined below:

1. use borehole data to analyze spatial variability in vertical direction,

2. use geological knowledge to assist in producing realistic parameters for the models of spatial variability, particularly for the lateral directions, and

3. consider using a Markov chain as a model of spatial variability because of its simplicity and relatability to fundamental observable attributes [13, 15, 20, 22].

The three-step approach outlined above will require a deeper understanding of the geostatistical model parameters beyond empirical fitting of mining-oriented parameters of “nugget, “range,” and “sill” [1, 2, 12].

Figure 5 shows the transition probability measurements and Markov chain model derived from the exhaustive data set in Figure 4(A) assuming isotropy. The measurements and models are organized in a matrix, with each entry representing the transition from the row category to the column category. The Markov chain fits what appear to be rather complex “S” and “hump” shaped curves that are not readily fitted by empirical curve-fitting procedures [26]. In the off-diagonal entries, the S-shaped structures are characteristic of two categories that tend not occur adjacent to each other, and the hump-shaped structures are characteristic of two categories that tend to occur adjacent to each other (a juxtapositional tendency). The cause of these transition probability structures is obviously attributable to the observable juxtapositional tendencies in the dataset:

Figure 5.

Matrix of measured and modeled transition probabilities for the dataset in Figure 4.

$$1\to 2\text{or}32\to 13\to 1\text{or}44\to 3.$$

A take-home message here is if a priori information exists on the juxtapositional tendencies, such as from geological knowledge of facies architecture, that insight could assist in the development of a realistic model of spatial variability.

The Markov chain, while appearing to be complex by its ability to produce variety of model curves – sixteen different transition probability model curves in this four-category application – is simply formulated as a matrix exponential function:

$$T\left(h\right)=\exp \left(Rh\right)$$E2

where h is the lag and R is a matrix of transition rates [13, 15, 20]. Importantly, while the original indicator simulation methods [1, 2, 26] would have required cumbersome indicator variogram analysis (construction of sample variograms and fitting of models) for each of the 16 plots, here all 16 are modeled with a single, matrix Eq. (1).

At the time we were developing transition probability approach to conditional simulation of hydrofacies architecture, a variogram modeling approach applied to the four facies of Figure 4 was shown to require individual models for ten (cross-) variograms [26], which may seem an advantage over modeling sixteen transition probabilities. However, for discrete categorical variables such as hydrofacies, the laws of probability will reduce the number of transition probability model curves required to model T(h) – from sixteen to nine for the four-category case – by imposing the row and column summing constraints shown on Figure 5 [13, 14, 15, 20]. Probability constraints are conveniently applied by choosing a “background” category, for which dark gray (category 2) is an obvious choice in Figure 4(A). Matrix entries in Figure 5 involving category 2 are shaded yellow to indicate constraint by the laws of probability. By assuming symmetry, the required number of transition probability model curves reduces to six for this application [15].

Mathematically, the 4 $\times$ 4 Markov chain R matrix of Eq. (1) for a four-category system with symmetric juxtapositioning will require specification of a total of nine parameters, six of which are transition rates and three of which are proportions. The proportions determine the “sill” on each column entry. The transition rates for the three unshaded diagonal entries are determined by the mean length ${\overline{L}}_k$ as $r_{kk}=-1/{\overline{L}}_k$, which is the slope of t_kk(h) at h = 0 [13, 14, 15, 20].

In summary, the 16 transition probability models shown in Figure 5 are formulated by the Markov chain using only nine parameters, all of which are relatable to observable attributes of proportions, mean length, and juxtapositional tendencies [13, 14, 15, 20]. This relatability to fundamental observable attributes enables a priori knowledge or geological insight to assist in development of a model of spatial variability.

3.3 Sand lens example

In some applications, hydrofacies may be defined from textural categories (e.g., clay, mud, silt, sand, muddy sand, gravel) with or without a depositional context. For example, intensive geostatistical analysis and modeling has been applied to the spatial variability of sand embedded within silt [27]. Figure 6(A) shows a “training image” of sand occurrences shaded in black. Different geostatistical methods were applied in [27] including a transition probability approach with two textural categories – sand and silt. However, this application of the transition probability approach did not distinguish what appear to be two sand types occurring as relatively small and large sand lenses that obviously differ in their observable attributes. As a result, the stochastic simulation results lacked an apparent bi-modal distribution of small and large sand lens occurrences.

Figure 6.

Example for hydrofacies defined by two sand lens categories.

In this application to textural data, we would recommend distinction of two categories of sand lens (SL): SL-large and SL-small. Figure 6(B) shows proportions and mean lengths estimated from the training image and tables provide in [27]. An evident juxtapositional tendency is that the sand lenses tend to transition to silt (white), which readily functions as the background category. We used estimates of sand lens proportions and mean length and a near-100% tendency of the sand lenses to transition to silt to formulate a three-category Markov chain model. Figure 6(C) shows a stochastic simulation result with the SL-small hydrofacies highlighted in red. By distinguishing the two size-related categories of sand lens, the transition probability approach could produce a bi-modal distribution of sand lens sizes consistent with the patterns of sand lens occurrence in the training image. This example emphasizes the importance of integrating the interpretative aspects of facies concepts into the stochastic model development.

4. Conditional simulation including curvilinear features

A conditional simulation algorithm may need to be extended to honor curvilinear features such as variable angles of dip or major direction of deposition or anisotropy. Figure 7 shows an example of stochastic simulation results from tsim including curvilinear dip angles [13]. The realization shown in Figure 7(A) is unconditional to borehole data and assumes zero or no dip at this site having a ground surface sloped from right to left. The realizations shown in Figure 7(B) and (C) are conditional to hydrofacies data at the left and right ends of the cross-section (highlighted in green) and dip angles. The constant dip case of Figure 7(B) assumes the dip is parallel to the ground surface. The variable dip case of Figure 7(C) used seismic data interpretations to determine the dip angles.

Figure 7.

Example of conditional simulation including dip angles.

To include curvilinear features, tsim uses a method similar to “local anisotropy kriging” [28], where the variation in the local major directions of anisotropy (or local stratigraphic coordinates) are prescribed a priori [13, 29]. One may analyze variable major directions of anisotropy as a separate stochastic process, an approach which has been implemented using tsim to simulate dipping or meandering characteristics of hydrofacies in a fluvial depositional system [20, 21].

Of practical relevance to tsim, a “checkerboarding” effect can result depending on the discretization or the magnitude of variation in the major directions of anisotropy. Undoubtedly, the numerical algorithms used by tsim have room for improvement! Nonetheless, these examples illustrate that there are working mechanisms to include curvilinear features in a transition probability approach to conditional simulation.

5. 3-D hydrogeological example applications

5.1 LLNL alluvial fan

The LLNL alluvial fan example application of geostatistical modeling [13, 14, 15, 19, 20] was originally motivated by need for a joint geologic and stochastic approach to the hydrogeological characterization of a contaminated aquifer system [3, 4, 5]. LLNL geologists recognized that the vertical hydrofacies sequences exhibited fining-upward cycles characteristic of fluvial deposition, which are manifested by non-symmetric juxtapositional tendencies [11, 13, 15, 17, 18]. A fundamental drawback to the variogram approach is the implicit assumption of symmetry in the juxtapositional tendencies inherent to the indicator (cross-) variogram spatial statistic [13, 14, 26].

The Markov chain was found to be a simple but effective model to fit observations of both the embedded and spatial aspects of measured vertical transition probabilities of the LLNL alluvial fan hydrofacies [13, 15, 20]:

• non-symmetric series of vertical juxtapositional tendencies (or fining upward cycle) of channel $\to$ levee $\to$ floodplain hydrofacies,

• variable thickness and lateral extent of the different hydrofacies, and

• expected lateral juxtapositional tendencies for a fluvial depositional system based on Walther’s Law concepts [17, 30].

Markov chain models were developed for the lateral directions using a combination of borehole data, outcrops, and geological insight, with the resulting 3-D Markov chain serving as the stochastic model of spatial variability for the hydrofacies architecture [15]. The example conditional simulation for the LLNL alluvial fan application shown in Figure 8 is conditioned to over two-hundred boreholes [13, 20]. Subsequent uses of tsim-generated conditional simulations of the LLNL alluvial fan hydrofacies architecture included several studies on the effects of heterogeneity on groundwater flow and transport [8, 10, 31, 32, 33, 34, 35].

Figure 8.

3-D conditional simulation of hydrofacies architecture of LLNL alluvial fan.

Consideration of the radial and dipping aspects of alluvial fan deposition was another important aspect of the geostatistical modeling effort applied to hydrogeologic characterization of the LLNL alluvial fan. Accordingly, the conditional simulations of the LLNL alluvial fan included curvilinear features of variable dip and radial variation in major direction of deposition, as evident in Figure 8. Curvilinear features were included in other early applications of tsim [13, 20, 21, 22] before the advent of multiple-point geostatistical methods for conditional simulation [36, 37].

5.2 Verification

Verification includes determination if a model is functioning as intended [38], which for geostatistical simulation may include visual inspection as in Section 3.1 or comparison of modeled and simulated spatial variability [13, 26, 39]. In Figure 9(A), conditional simulation results conditioned to data along four vertical strings are compared for the variogram and transition probability approaches using sisim and tsim, respectively. The variogram approach produces less spatial continuity and less-ordered juxtapositioning of the hydrofacies as compared to transition probability approach. Figure 9(B) compares modeled and simulated transition probabilities to a Markov chain fitted to extensive vertical transition probability data at LLNL [13, 14, 15, 20]. The transition probability approach produced simulation results that closely match the Markov chain, including non-random asymmetric juxtapositional tendences (e.g., channel $\to$ bar/levee). The variogram approach was unable to match the Markov chain model, notably for proportions of hydrofacies and asymmetry in juxtapositional tendencies.

Figure 9.

Comparison of results for variogram and transition approaches for (A) stochastic simulation and (B) simulated transition probabilities.

5.3 South Lake Tahoe contaminant transport example

A geostatistical model or simulation is usually not the end-product of hydrogeological studies; applications to groundwater problems are oriented toward gaining understanding of flow and transport processes and insight to interpretation of measurements at wells [3, 4, 5, 6, 8, 10, 13, 21, 31, 32, 33, 34, 35, 39, 40, 41, 42, 43]. In the late 1990s, the gasoline additive Methyl tert-butyl ether (MTBE) was reaching drinking water wells in the South Lake Tahoe area of California at levels approaching the regulatory standard of 40 ug/L. Understanding how the MTBE plumes behave was highly relevant to design and implementation of effective remediation strategies and for protecting freshwater supplies [44].

A potential application of geostatistics is to use 3-D conditional simulation to assess effects of subsurface heterogeneity on transport of a mobile contaminant such as MTBE from source to well. For this example, a flow and transport model was developed to investigate possible effects of a point-source release of MTBE in a geological setting representative of the South Lake Tahoe area [40]. The geological information available indicated that this area is underlain by co-mingling glacial and fluvial sediments, with glacial sediments prevalent in the upstream direction. Till, outwash, and lacustrine hydrofacies were identified within the glacial sediments, and fluvial hydrofacies were largely identified by textural information.

Figure 10(A) shows a 3-D geostatistical realization developed for this South Lake Tahoe application using the transition probability approach at two scales:

1. comingling depositional systems – glacial and fluvial, and

2. hydrofacies within each depositional system.

Figure 10.

3-D geostatistical realizations for the South Lake Tahoe example: (A) hydrofacies and (B) hydraulic conductivity.

Hydraulic properties were spatially distributed by Gaussian random fields superposed on the hydrofacies. Figure 10(B) shows a hydraulic conductivity realization that is conditioned by the hydrofacies realization in Figure 10(A). The total number of grid cells is 7.2 million.

The geostatistical realizations of hydraulic properties were used to parameterize simulations of groundwater flow and transport. The steady-state groundwater flow field shown in Figure 11(A) was simulated using the ParFlow code [45], and the transport was simulated using a streamline-based particle-tracking algorithm [46] with added dispersive and diffusive transport [47]. Transport simulation results shown in Figure 11(B) assume a 1.5 kg/day release of a mobile contaminant for 300 days. The simulation results indicate the possibility of the contaminant reaching one or more drinking water wells within a few years depending on pumping schedules and the heterogeneous distribution of hydraulic properties, which imparts significant macro-dispersive effects. Within the stochastic framework, these simulation results are consistent with the observations.

Figure 11.

Flow and transport results for the South Lake Tahoe example: (A) Simulated steady-state flow field, and (B) simulated plume and contaminant breakthrough at wells.

What is not apparent in the simulated plume picture of Figure 11(B) is that the contaminant transport initially migrates rapidly along preferential pathways in the shallow subsurface and, later, pumping effects draw the contaminant downward through a complex network of interconnected pathways of relatively high hydraulic conductivity. The geostatistical components of this model were essential to understanding how the effects of geological heterogeneities on groundwater flow and transport processes may explain observations of elevated levels of contaminants in wells.

6. Summary and conclusions

The transition probability approach to conditional simulation of hydrofacies architecture was developed in the mid-1990s as an extension of existing variogram-based geostatistical approaches, which were then state-of-the-art [1, 2]. At that time, interest was growing in applying categorical simulation approaches to hydrogeological problems [3, 4, 5, 6]. The transition probability approach was conceived to facilitate integration of geologic interpretation and improve the consistency between observed and modeled spatial variability [13, 14, 15, 19, 20].

In geostatistical applications to hydrogeological problems, the conceptualization and/or identification of hydrofacies should be rooted in both geologic and geometric contexts, with attention to the relevant scales and stochastic aspects of subsurface characterization [8, 10, 11]. In some applications, hydrofacies may be defined by distinctive depositional or hydrogeological units with characteristic textural assemblages; in others textural categories may be used but may need to be either grouped or further subdivided in a geologically reasonable manner to account for the heterogeneity effects of interest.

In hydrogeological applications of conditional simulation, uncertainty can be reduced by separating out the deterministic and stochastic aspects of the problem [20, 21, 22, 48]. Deterministic aspects include conventional geologic and geophysical stratigraphic analysis to identify and map the major formations, depositional systems, or hydrogeological zones, including their bounding unconformities and geologic structures. Stochastic aspects include unknown spatial variability within deterministic units and may extend to structural variations or curvilinear features, such as for folded stratigraphy or meandering directions of deposition [13, 20].

Different geostatistical methods, and not necessarily a transition probability approach, should be considered depending on the hydrogeological problem of interest. The transition probability approach has been useful for gaining insight on the conceptual aspects of flow and transport processes, such as:

• Connectivity of the highest permeability hydrofacies can strongly influence flow and transport processes and, therefore, connectivity is itself a valuable metric for assessing the relative suitabilities of different models of heterogeneity [8, 10, 21, 31, 39, 40].

• The abundance and architecture of low-permeability hydrofacies are key factors in governing the behavior of contaminant breakthrough, timing and efficacy of pump-and-treat remediation, and diffusive fractionation of groundwater isotopes [32, 33, 34, 35, 43].

However, it should be recognized that in some applications, different geostatistical methods may produce effectively similar outcomes in assessment of the groundwater flow and/or transport processes of interest [41].

Arguments can be made to add further complexity (e.g. cross-bedding or surficial worm-like “channel” structures) into geostatistical models. Subsequently, more complex multiple-point geostatistical approaches have been developed [34, 35]. However, in some situations the benefit of adding further complexities to geostatistical models do not outweigh the costs, which may include acquisition of extensive “training” data sets [41, 49]. Moreover, the added effort and complexity may not be necessary for interpretation at the scales of observation (e.g., bulk aqueous concentration or hydraulic head in wells) or be geologically realistic considering that sedimentary features on the surface are not necessarily preserved in the subsurface [50].

There remains room for improvement in the practicality, flexibility, or complexity of geostatistical conditional simulation algorithms and codes. In hydrogeological applications to groundwater problems, the usefulness of a geostatistical model may be best measured by how it enhances understanding flow and transport processes and reduces uncertainty in interpretation or prediction of observational data, particularly at wells.