An inverse modeling approach - Wiley Online Library

WATER RESOURCES RESEARCH, VOL. 44, W08426, doi:10.1029/2007WR006635, 2008

Identifying discrete geologic structures that produce anomalous hydraulic response: An inverse modeling approach Michael J. Ronayne,1 Steven M. Gorelick,2 and Jef Caers3 Received 7 November 2007; revised 8 April 2008; accepted 30 May 2008; published 20 August 2008.

[1] Subsurface hydraulic response is controlled by discrete geologic structures in a

variety of aquifer settings. In this work, we analyze data from the well-studied alluvial fan system that underlies the Lawrence Livermore National Laboratory. There, subsurface flow and transport behavior is strongly influenced by discrete channel deposits, which are embedded within less permeable floodplain deposits. We model the distribution of these deposits using a multiple-point geostatistical approach that enables simulation of geologically plausible connected permeable bodies consistent with a training image. By coupling the geostatistical model with a dynamic flow model in a simulation inversion framework, we identify specific channel structures that are consistent with transient aquifer test data at six observation wells. We perform shortest-path analysis to characterize the geometry of the simulated channels and to locate critical high-permeability conduits that control the system’s response to pumping. In one area, all successful inverse solutions identify the same set of stacked channel deposits, a geologic feature that provides an explanation for anomalous response behavior observed during the aquifer test. Citation: Ronayne, M. J., S. M. Gorelick, and J. Caers (2008), Identifying discrete geologic structures that produce anomalous hydraulic response: An inverse modeling approach, Water Resour. Res., 44, W08426, doi:10.1029/2007WR006635.

1. Introduction [2] Subsurface structures exert a primary influence on hydraulic response dynamics in many geologic settings. Examples of discrete structures characterized by high or low permeability include faults, fractures, and depositional features such as coarse-grained channel fills. The spatial arrangement of these features often represents the first-order (i.e., most hydrologically significant) heterogeneity since their presence tends to produce large discontinuities in the permeability distribution [Fogg, 1986; Western et al., 2001; de Marsily et al., 2005]. Identifying and quantifying the impact of dominant structures at all relevant scales represents an important, continuing, research challenge [Proce et al., 2004; Schulz et al., 2006]. [3] In geologic environments where discrete high-permeability (k) features are present, the connectedness and geometry of these features is critical. Many recent studies in the groundwater literature have emphasized the fact that, when the most permeable subsurface material is highly connected, flow channeling is enhanced and solute transport behavior often belies predictions based on the classical advective-dispersive model [Zinn and Harvey, 2003; Zheng and Gorelick, 2003; G. Liu et al., 2004; Knudby and Carrera, 2005; Feyen and Caers, 2006]. 1 Department of Geological and Environmental Sciences, Stanford University, Stanford, California, USA. 2 Department of Environmental Earth System Science, Stanford University, Stanford, California, USA. 3 Department of Energy Resources Engineering, Stanford University, Stanford, California, USA.

[4] Although the availability of a connected path of highpermeability material is key, the geometry of that structural pathway is also important. Ronayne and Gorelick [2006] investigated fluid flow in synthetic systems spanned by high-k channel networks and demonstrated that the degree of flow channeling was very sensitive to the network geometry. In the study of porous composites, the geometry of permeable inclusions is often found to be an important influence on effective transport properties [Torquato, 2002]. In fractured rock settings, the orientation and arrangement of fractures are strong controls on flow channeling and the extent of fracture-matrix fluid exchange [Taylor et al., 1999; Matthaï and Belayneh, 2004]. [5] Percolation theory offers a general framework for studying connections and the geometry of connected material. Percolation models are commonly used to represent heterogeneous porous media [Sahimi, 1993; Hunt, 2001], and there are many useful results in the literature that relate fluid flow behavior to the geometrical characteristics of percolation clusters. One particularly important property is the shortest path; the length and geometry of the shortest path between sites on a cluster has been studied for various percolation models [Herrmann and Stanley, 1988; Dokholyan et al., 1998; Paul et al., 2002]. Knowledge of these properties is useful for understanding dynamics on permeable clusters. For example, it has been demonstrated that the minimum transport time for conservative solute migration between two locations on a cluster is strongly correlated with the shortest path length and fractal dimension [Lee et al., 1999; Andrade et al., 2000; Arau´jo et al., 2002]. These and other related theoretical results provide

Copyright 2008 by the American Geophysical Union. 0043-1397/08/2007WR006635

W08426

1 of 16

W08426

RONAYNE ET AL.: IDENTIFYING DISCRETE STRUCTURES

insight for field studies that aim to find a linkage between observations of hydraulic response and the properties of subsurface structures. [6] In this work, we consider data from the Lawrence Livermore National Laboratory (LLNL) field site in northern California. There, subsurface flow occurs within an alluvial sequence containing channel deposits. The permeability of the channel deposits is significantly higher than that of the adjacent overbank and floodplain material. Carle [1996] was among the first to explore how the structural heterogeneity at LLNL impacts hydraulic response behavior. Fogg et al. [2000] investigated the connectivity of the channels and proposed a ‘‘connected-network paradigm’’ to describe the heterogeneity, which may be broadly applicable to aquifers of fluvial origin. In this work, we build upon those previous studies and present a new quantitative method for identifying channel structures that control subsurface hydraulic response. This study is different from previous studies in that we use multiple-point geostatistics to develop high-resolution numerical models of the subsurface heterogeneity at LLNL. Also new is our application of a Bayesian inverse modeling technique, which allows us to identify specific channels that explain complex hydraulic response data measured during a site aquifer test. Previous investigators, most notably Fogg et al. [2000], have concluded that the aquifer test response behavior is a consequence of strong heterogeneity and a high degree of channel connectivity within the alluvial system. By applying an inverse model and performing shortest-path analysis, we go a step further with this work; namely, we locate specific high-k conduits, formed by interconnected channel deposits, that control the system’s response to pumping. Since our inverse solutions demonstrate good agreement with the measured hydraulic heads at several observation wells, and because we identify specific structures that explain observations of anomalous hydraulic response, the results of this study enhance the understanding of groundwater dynamics at LLNL. The inverse modeling technique that we employ, coupled with the shortest-path analysis, may be useful in a variety of geologic settings where the heterogeneity is dominated by discrete structures. [7] The remainder of this paper is organized as follows. In section 2, we review the hydrogeology of the LLNL site and describe the aquifer test data set. Section 3 covers the geostatistical simulation and inverse modeling techniques. In section 4, we describe our groundwater flow model, which we use to analyze the aquifer test observations. Results are presented in section 5 and, in section 6, we highlight our conclusions and discuss some avenues for future research.

2. Study Area and Site Data [8] Field data used in this study are from an alluvial fan deposit that underlies the Lawrence Livermore National Laboratory. Located in northern California (Figure 1a), LLNL is situated along the boundary between the Livermore Valley and the Diablo Mountains, part of the California Coast Ranges. In addition to ongoing hydrogeologic investigations to support remediation activities, the Livermore site has been the subject of numerous research studies. Of particular relevance to this work are previous efforts to model the subsurface heterogeneity using a discontinuous facies

W08426

framework [Carle, 1996; Carle and Fogg, 1996, 1997; Fogg et al., 2000]. In these studies, a transition probability Markov chain method was used to obtain stochastic realizations of the subsurface architecture. Recently, Lee et al. [2007] compared permeability fields generated with this method to fields simulated using a multivariate Gaussian model. They found that models constructed using transition probabilities were better able to explain the general response behavior observed during a site aquifer test. However, in the previous modeling efforts, no single realization has successfully reproduced the response data at a majority of observation wells. 2.1. Geology [9] Sediments that make up the LLNL alluvial fan system were delivered by streams originating in the Diablo Mountains. Consequently, the present fan system consists of discrete channel deposits embedded within finer-grained and less permeable floodplain deposits [Weissmann et al., 2002]. The floodplain material occupies a significantly higher volume fraction than the channel material. Previous conceptual models have envisaged a labyrinth geometry [Weber and van Geuns, 1990] to describe the distribution of channels. This interpretation is supported by the borehole data and is also consistent with general expectations for the depositional setting. [10] Noyes [1990], Carle [1996], and subsequent investigators [e.g., Fogg et al., 1998; Weissmann et al., 2002] demonstrated that textural descriptions and borehole geophysical data could be used to differentiate channel deposits from other finer-grained deposits within the alluvial system at LLNL. Following the general approach taken in those previous studies, we analyzed core descriptions and geophysical logs (single-point resistance and gamma logs) from 50 boreholes throughout the study area (extent shown in Figure 1b). In addition to estimating the fraction of each material type, the borehole data allowed us to determine realistic dimensions for the channel deposits. Channels comprise the coarsest subsurface material at LLNL. From textural descriptions, these deposits are identified as gravels, sandy gravels, and gravelly sands. On geophysical logs, channels are characterized by high resistivities and relatively low gamma signatures. Sandy silts, clayey silts, and silty sands are the dominant textures in the floodplain areas. These finer-grained deposits are distinguished on geophysical logs by moderate to low resistivities and gamma ray levels that are higher relative to the channels. Figure 2 shows example logs for two site boreholes. By jointly analyzing the geophysical data and textural descriptions from core samples, discrete channel intervals are identified with a high degree of confidence. Our analysis is consistent with the aforementioned previous studies, indicating that the volume fraction of the channel material is in the range of 15 –20%. For the subsurface area and sampling locations that we considered, the average channel thickness is 1.3 m; this is identical to the average thickness reported by Carle [1996]. [11] The boreholes shown in Figure 2 are separated by only 3.9 m and are situated along a line that is transverse to the main paleoflow directions. It is noteworthy that similar channel intervals are observed along both boreholes. We found this to be a common occurrence when analyzing data from nearby (20 m) that separate most boreholes. However, the successful correlation of channels for several closely spaced borehole pairs was important, not only to corroborate the conceptual notion that the channels are key depositional features that occupy significant volumes, but also to approximate a reasonable lower bound on the typical channel width. [12] One other important characteristic of the subsurface heterogeneity is the presence of laterally continuous paleosols. Given their high clay content and extremely low permeability, paleosol deposits can create extensive barriers to flow. Recent hydrogeologic analysis has led to the determination of several different hydrostratigraphic units (HSUs) within the alluvial sequence at LLNL [Mansoor et al., 2002; Hoffman et al., 2003]. Distinct HSUs have limited hydraulic communication, likely a result of bounding paleosol surfaces [Weissmann et al., 2002; Bennett et al., 2006]. In our model of the subsurface architecture, we include major paleosols (see section 4), which are key features affecting flow and transport behavior. [13] To summarize, the conceptual geologic model that we adopt in this study comprises three facies: (1) channel, (2) floodplain, and (3) paleosol deposits. This differs from previous efforts that rely on four depositional facies (channel, levee, debris flow, floodplain) to describe the subsurface heterogeneity at LLNL [Carle, 1996; Carle and Fogg, 1996, 1997; Fogg et al., 1998, 2000]. We contend that the simpler channel-floodplain dichotomy, along with the inclusion of major paleosol deposits, captures the most important aspects

of the heterogeneity controlling hydraulic response behavior. As will be discussed in section 4.2, all three of our modeled facies are characterized by substantially different hydraulic properties. We note that much of the evidence indicating the presence of laterally extensive paleosols is based on recent

Figure 2. Geophysical logs for two neighboring site boreholes. Solid and dashed lines are the resistivity and gamma logs, respectively. Shaded intervals represent coarse-grained textures (sands and gravels) identified from core samples.

3 of 16

W08426


W08426

2.3. Anomalous Response Curves [16] Spatially variable drawdowns that result from pumping are governed by a diffusion equation Ss ðxÞ

Figure 3. Composite drawdown plots for (a) observation wells screened in the upper unit, HSU-2; and (b) observation wells screened in the lower unit, HSU-3A. The x axis shows time (since the start of pumping) divided by distance squared; for each well, r is the distance from pumped location. Drawdown, s, is plotted on the y axis. data collection and analysis performed after publication of many of the previous studies. 2.2. Hydraulic Response Data [14] In this study, we analyze subsurface hydraulic response data from an aquifer test conducted at the LLNL site in April 1994. Groundwater was extracted from a single pumping well for a duration of 24 h. The initial pumping rate was 98.5 m3/d and, after 400 min, was reduced and maintained at approximately 69.1 m3/d. Pressure transducers were installed to record transient hydraulic heads at multiple observation wells. A more detailed account of the aquifer test is given by Carle [1996]. [15] Our analysis is focused on six observation wells (well locations are shown in Figure 1b). We corrected the observed hydraulic head data for atmospheric pressure fluctuations. Additionally, for two of the deep wells (W712 and W267), it was necessary to correct for a prevailing downward trend in the head data. A detailed discussion of the observation data set and the corrections that we apply is provided with the auxiliary material.1 Hereafter, we present only the corrected observation data. We will adopt the term drawdown to refer to the amount of head decline relative to the initial hydraulic head measurement at the start of pumping. 1 Auxiliary materials are available in the HTML. doi:10.1029/ 2007WR006635.

@hðxÞ ¼ r ð K ðxÞrhðx; t ÞÞ þ Qðx; t Þ @t

ð1Þ

where h is the hydraulic head, Ss is the specific storage, K is the hydraulic conductivity, K = (krg / m), and Q represents a source/sink term. Insight is often gained by considering the hydraulic diffusivity, D = K/Ss. In response to pumping, the propagation rate of a given magnitude of drawdown is proportional to the hydraulic diffusivity [Streltsova, 1988]. [17] Composite drawdown plots for the six observation wells are shown in Figure 3. The individual response curves are separated on the basis of each well’s screened interval. Three observation wells are screened in the upper hydrostratigraphic unit, HSU-2, and three are screened in the lower unit, HSU-3A. [18] It is useful to compare the field observations shown in Figure 3 to theoretical response curves that are obtained using well-studied idealized models. In the case of a medium with uniform D, a drawdown pulse moves through the system in a manner that is directionally invariant. For a confined aquifer system with uniform D and boundary conditions consistent with assumptions employed in the Theis [1935] solution, data points (s versus t/r2) from observation wells at various distances, r, would collapse onto the same curve. Inspection of Figure 3 reveals that this is clearly not the case for the LLNL aquifer test data. Departure from ideal behavior is especially pronounced for the HSU-2 observation wells. Of particular interest is the quick response and relatively large drawdowns recorded at W264, the most distant observation well. [19] Analyzed together, we note that the observed drawdown curves may be regarded as ‘‘anomalous’’ given that the response behavior differs so markedly from classical predictions [Cortis and Knudby, 2006]. In this work, we present a new spatially distributed model that relates this interesting hydraulic response to the locations, geometry, and connectivity of key channel features.

3. Geostatistical Methodology [20] Multiple-point geostatistics provides an effective framework for characterizing subsurface hydraulic property distributions that are dominated by discrete structures. The essence of the method is the application of higher-order statistics to capture patterns in heterogeneous fields. The ability to identify and reproduce complex multipoint patterns permits accurate modeling of features characterized by some irregular geometry. [21] Recent field applications of multiple-point geostatistics have focused on spatially variable petrophysical properties in petroleum reservoirs. For example, Strebelle et al. [2003] and Y. Liu et al. [2004] applied the technique to model sand body geometry in offshore clastic reservoirs. [22] The multiple-point method is well suited for modeling the subsurface channel distribution at LLNL given its emphasis on discrete structures and their geometry. Moreover, since the geometry of individual discrete objects affects the long-range connectivity properties for groups

4 of 16


W08426

of those objects [Garboczi et al., 1995; Baker et al., 2002], the multipoint approach offers some advantages (compared to two-point indicator geostatistical techniques) for realistically modeling the dynamic effects of connectivity. The ‘‘groups of objects’’ that are of interest in this study are connected channel deposits (or channel clusters), which will be discussed in section 5.3. [23] A general two-point correlation function (e.g., 0 R(a) 2 (u, u )) considers the probability of finding the same facies identifier, a, at locations u and u0. With the multipoint approach, higher-order correlation functions are assembled to characterize the spatial variability more completely. Information gathered from the subsurface is too sparse to directly infer these higher-order statistics. Instead, a so-called 3-D training image is created by interpreting the 3-D geologic structure from subsurface data. For example, if borehole data, geophysical records, and geologic interpretation indicate a fluvial depositional environment, then a training image model is built to describe general characteristics of the channel deposits. The training image only depicts the conceptual model and need not locally match the actual subsurface data. In addition to site-specific information gained from geologic sampling, other relevant data such as outcrop and analog data may be used to constrain the geometry and orientation of the channel features. [24] Several stochastic simulation methodologies have been developed to generate 3-D realizations that reproduce the style of heterogeneity depicted in a training image and, at the same time, are locally constrained to data [Strebelle, 2000; Zhang et al., 2006; Arpat and Caers, 2007]. In this paper we follow Strebelle’s approach, where multiple-point data events (essentially patterns) are extracted from the training image and anchored to subsurface data. The concept of the data event is critical for practical implementations of the technique [Strebelle, 2002; Okabe and Blunt, 2005]. [25] Our geostatistical simulation model for the LLNL site considers two dominant categories: (1) the channel facies and (2) the floodplain matrix facies. This discontinuous facies framework requires a discrete space random function (SRF) to describe the heterogeneity [Johnson and Dreiss, 1989]. Here, we use an indicator SRF defined as I ðu Þ ¼

1; 0;

if u 2 C otherwise:

m ¼ f I ðu1 Þ; . . . :; I ðuN Þg ¼ fI1 ; . . . :; IN g

where N is the total number of grid cells. The prior probability distribution, which describes the dependencies among those model parameters, can therefore be expressed as f ðmÞ ¼ Prf I ðu1 Þ ¼ iðu1 Þ; . . . :; I ðuN Þ ¼ iðuN Þg ¼ PrfI1 ¼ i1 ; . . . :; IN ¼ iN g

ð3Þ

ð4Þ

In our case, the full multivariate prior is not explicitly known; only conditional probability distributions are specified. Indeed, using the exact decomposition f ðmÞ ¼ PrfI1 ¼ 1g PrfI2 ¼ 1 j iðu1 Þg . . . : PrfIN ¼ 1 j iðu1 Þ; . . . :; iðuN 1 Þg

ð5Þ

only the univariate conditional distributions, describing the probability of each facies at a specific location, are required. In the approach that we use, these conditional distributions are modeled directly from the 3-D training image [Strebelle, 2000]. Sequential simulation [Johnson, 1987; Deutsch and Journel, 1998] is used to draw realizations. Realizations (or samples) drawn from this prior are termed unconditional realizations. [27] Figure 4 provides a simple schematic demonstrating how patterns in a training image are used to formulate the conditional probability distribution at a grid location. A conditional distribution is required at every cell in the domain where the facies is to be simulated. With the sequential simulation approach, a random path is defined that determines the order in which each grid cell is visited. As a result, the conditional distributions are constrained by more information as the simulation proceeds. Multiple realizations are required to explore the sensitivity of model results to differing random paths. Our inverse modeling procedure, described in section 3.1, blunts this sensitivity because the random path is routinely altered during the search for an optimal solution. [28] An important attribute of any practical geostatistical simulation method is the ability to condition to local site data. The following posterior probability distribution includes an observed data vector, d1, which contains n locations (borehole intervals) where the facies is known f ðmjd1 Þ ¼ PrfI1 ¼ i1 ; . . . :; IN ¼ iN jd1 g

ð2Þ

where u is the space coordinate and C refers to the set of all channels. If the function I(u) takes on a value of 1, then the channel facies is present at that location. Otherwise, the floodplain facies is present. In this study, we use multiplepoint statistics to describe the variability of the indicator SRF in (2). Notice that we do not include a separate paleosol facies in (2). Instead, we represent major paleosol deposits deterministically, which will be discussed in section 4. [26] For a particular simulation domain, the complete set of model parameters is given by [Caers and Hoffman, 2006]

W08426

ð6Þ

where d1 = {i(ua), a = 1, . . .., n}. Conditional realizations are generated by drawing samples from the posterior distribution in (6). Conditioning to direct observations of the indicator variable (observed data in d1) is performed using the same sequential simulation methodology [Deutsch and Journel, 1998]. This conditioning is straightforward since the direct data are treated the same as previously simulated values (Figure 4). 3.1. Integration of Dynamic Data [29] We also condition the geostatistical realizations to dynamic hydraulic response data. That is, we perform inverse modeling in an effort to find facies distributions that are consistent with transient hydraulic heads measured during a site aquifer test. These dynamic data, d2, which are

5 of 16


W08426

W08426

Figure 4. Two-dimensional example demonstrating the use of multiple-point data events to formulate a conditional probability distribution for location uj. The known values at u1, u2, u3, and u4 could be previously simulated locations, visited earlier in the sequential simulation procedure, or direct field observations (facies observed in site boreholes). For the specific configuration under evaluation (de), there are four replicates of this data event identified in the training image. In three of those four replicates, the central grid cell is occupied by a channel (i = 1). indirectly (nonlinearly) related to our indicator variable, can be incorporated within a Bayesian framework f ðmjd1 ; d2 Þ ¼

f ðd1 ; d2 jmÞf ðmÞ f ðd1 jmÞf ðd2 jmÞf ðmÞ ffi f ðd1 ; d2 Þ f ðd1 ; d2 Þ

ð7Þ

A forward operator F explains the functional relationship between the dynamic data and the model parameters; d2 = F(m) [Tarantola, 2005]. The last expression in (7) holds if we assume that d1 and d2 are conditionally independent. [30] Equation (7) is the full posterior distribution that includes both direct and indirect data types. The traditional sampling approach calls for explicit analytical specification of likelihood and prior. Except for the multi-Gaussian framework and linear data, such specification becomes analytically intractable or requires Markov chain Monte Carlo techniques [Mosegaard and Tarantola, 1995; Omre and Tjelmeland, 1996] that are prohibitively expensive when the forward model F is a high-resolution flow simulator. [31] Instead, we follow an alternative approach specified in detail by Caers and Hoffman [2006] that relies on a decomposition of the posterior in terms of two ‘‘preposterior’’ distributions. The preposteriors describe the dependence of an occurrence of I (the indicator variable of a grid cell) on each data type separately. We introduce some additional notation to further explain the method: Pr I uj ¼ 1jiðu1 Þ; . . . :; i uj 1 ; d1 ; d2 ¼ Pr Aj jBj ; d2 ð8Þ

where Aj = {I(uj) = 1} and Bj = {i(u1), . . .., i(uj 1), d1}. The distribution Pr{Aj j Bj, d2} in (8) is the local conditional probability distribution for the jth grid cell. It considers both data types, d1 and d2, as well as previously simulated facies values for the other (j 1) cells. To facilitate estimation of this distribution during sequential simulation, two preposteriors are defined, Pr{Aj j Bj} and Pr{Aj j d2}. The preposterior Pr{Aj j Bj} is modeled directly from the training image using the general approach illustrated in Figure 4. The preposterior for the dynamic

data, Pr{Aj j d2}, is estimated iteratively using the probability perturbation method [Caers, 2003; Caers and Hoffman, 2006]. [32] The probability perturbation method simplifies the inversion problem by focusing on one (global) or a few (regional) perturbation parameters. The following equation explains the effect of the perturbation parameter: h i Pr Aj jd2 ¼ ð1 rD Þ ið0Þ uj þ rD Pr Aj j ¼ 1; . . . :; N

ð9Þ

where rD is a perturbation parameter, i(0)(uj) is the simulated facies value corresponding to the previous best matching realization, and Pr{Aj} is the marginal probability. An rD value of 0 means that the existing facies should be maintained at location uj. An rD value of 1 corresponds to a maximum perturbation, such that Pr{Aj j d2} is set equal to the overall channel fraction. To find an optimal value of rD, flow simulations are required. After each simulation, an objective function is evaluated, and then the value of the perturbation parameter is updated in an effort to minimize the objective function. In our case, the objective function quantifies the mismatch between observed and modelsimulated drawdowns. We discuss this further in section 5.2. A flowchart describing the inverse procedure is provided with the auxiliary material. [33] To increase the efficiency of the inverse model, we optimize on two (regional) perturbation parameters simultaneously. We use one rD value to perturb the local probabilities in HSU-2 and a separate rD value for model cells in HSU-3A, the lower hydrostratigraphic unit. With the sequential simulation framework, the use of multiple perturbation parameters does not produce any artifacts [Hoffman and Caers, 2005]. That is, a sample from (7), without discontinuities between regions, can still be obtained. [34] We use the method suggested by Journel [2002] to combine the preposteriors and estimate Pr{Aj j Bj, d2}. In practice, Pr{Aj j Bj} is formulated during sequential simulation and combined with Pr{Aj j d2} to obtain Pr{Aj j Bj, d2}. Journel [2002] and Hoffman [2005] provide additional

6 of 16


W08426

Figure 5. Training image. details and further justification of this approach for similar applications.

4. Model Development [35] A key goal of this work is to provide a physical explanation for anomalous response behavior observed during the aquifer test. To accomplish this, we require a structural model of the subsurface heterogeneity and a dynamic model for simulating the response to pumping. The first part of this section addresses the heterogeneity model. Example realizations of the subsurface architecture are presented and discussed. In the second part of this section, we describe our groundwater flow modeling approach. 4.1. Subsurface Architecture [36] We use a training image-based multiple-point geostatistics approach to model the spatial distribution of the two dominant facies. The training image constrains the prior model. For the alluvial deposit investigated here, our training image contains specific information on the following properties: (1) the relative fractions of the channel and floodplain facies, (2) the thickness and (3) width of channel deposits, (4) the channel orientation (paleoflow direction), and (5) the sinuosity of channels. The first two items are estimated directly from available geologic data. As noted in section 2.1, the average channel thickness estimated from borehole samples is 1.3 m. Since the boreholes arbitrarily intersect portions of channel deposits that in reality are characterized by some curvilinear cross-sectional geometry, the sampled thickness is expected to be less than the maximum channel thickness. We use 2 m as a representative maximum thickness. Site data from closely spaced boreholes indicate that the typical channel width is at least 5 m. Drawing from an extensive review of data from the sedimentology literature, Gibling [2006] notes that many channel bodies preserved in alluvial fans are characterized by a ‘‘broad ribbon’’ geometry with width-to-thickness (W/T) ratios from 5 to 15. We assume W/T = 10, meaning that the maximum target width for an individual channel deposit is 20 m. On the basis of modern drainage patterns around LLNL, it appears that two channel orientations are likely, east-to-west and southeast-to-northwest [Hoffman and Dresen, 1990; Weissmann et al., 2002]. For the work presented here, we focus on the latter direction (i.e., we specify channels that trend to the northwest). The presence of channels in this direction was strongly indicated by the aquifer test response data that we studied. Perhaps the most uncertain of the properties listed above is the channel sinuosity. In our training image, we assume channel depos-

W08426

its with relatively low sinuosity. Given typical gradients that characterize a fan surface, highly sinuous channels are not anticipated, except possibly along the most distal reaches. [37] We used the Boolean technique of Deutsch and Tran [2002] to develop our 3-D training image, shown in Figure 5. Note that the training image is not conditioned to any local geologic samples or hydrologic observations; it represents a conceptual description that includes some quantitative information regarding the facies geometry. The size of the full training image is 500 m 500 m 16 m in the x, y, and z directions. In Figure 5, a subset of this volume is displayed (500 m 375 m 16 m). As demonstrated for the single layer (horizontal slice) shown in Figure 5, all channels in the training image are modeled as being laterally continuous. We permit some variability in the channel orientation, given the potential for local variations in paleoflow direction. More details on the training image generation, including input parameters for the Boolean simulator, are provided with the auxiliary material. [38] We generate 3-D realizations of the facies distribution using the algorithm of Strebelle [2000], with the 3-D training image as input. More specifically, multiple-point data events identified during a scan of the training image are used to inform simulated realizations of the facies distribution. In addition to data events (patterns) from the training image, other information used to guide the sequential simulation procedure includes direct conditioning data (observed facies from boreholes) and indirect (drawdown response) data. At specific locations and depth intervals sampled by boreholes, discrete facies (i.e., channel or floodplain material) were identified using the geologic analysis procedure described in section 2.1. These samples, which are used to condition the geostatistical realizations, constitute the d1 vector for our application. Sample locations are shown in Figure 1b. Along each sampled borehole, there are multiple depth intervals where the facies occurrence is prescribed. It should be noted that data were also available to support facies identification at the pumping well, W612, and at all but one of the observation wells. Accurate determination of the facies distribution is particularly important along the pumping well’s screened interval, since structural heterogeneity around the pumped location strongly controls early time response dynamics throughout the affected volume [Leven and Dietrich, 2006]. The data for W612 indicate that the well screen intersects a channel deposit. In fact, given the emphasis on remediation and the need to implement efficient pump-and-treat methods, most wells at LLNL are intentionally screened across high-k channel material. The combined physical/engineered alluvial system could therefore be described as a ‘‘plumbed’’ network of interconnected channel deposits contained within less permeable sediments [Fogg et al., 2000]. This suggests that some characteristics of the aquifer test response dynamics may be similar to those observed in fractured rock settings, where packers are used to isolate specific high-k fracture zones. [39] Two example conditional realizations from the geostatistical model are provided in Figure 6. We discretized the model space using a grid comprising uniform cells with dimensions Dx = Dy = 2.5 m and Dz = 0.4 m. The overall simulation domain contains 7,296,000 cells. This large, high-resolution grid was required in order to explicitly

7 of 16

W08426


W08426

Figure 6. (a and b) Two example realizations of the simulated facies distribution. Pumping well (circled) and observation well locations are shown in white. Only observation well W264 is screened across the displayed horizontal slice. The simulated volume shown here spans 525 m in the x direction, 400 m in the y direction, and 38 m in the z direction; y is the north-south direction. (c and d) Simulated versus observed drawdown at wells W264 and W714, respectively. In these drawdown plots, the solid black line gives the simulated response for the facies distribution shown in Figure 6a, which is an inverse solution. The other response curves correspond to noninverse realizations (i.e., samples from equation (6) that are conditioned only to direct data); the volume in Figure 6b depicts one such realization. represent individual channel deposits throughout the subsurface volume (>107 m3) that we consider. The volume selected for presentation in Figure 6 is a subset ( 45%) of the total simulated volume. [40] One noticeable characteristic of the realizations shown in Figure 6 is that individual channels are less (horizontally) continuous as compared to those included in the training image. This apparent difference in channel continuity is primarily a result of ergodic fluctuations [Deutsch and Journel, 1998] and has been described in more detail for multiple-point geostatistics by Caers and Zhang [2004]. For our work, practical considerations of

computational efficiency placed limits on the training image size; we used a training image that is slightly smaller than the simulation grid. Consequently, the simulation algorithm has difficulty reproducing the longest-range patterns because adequate sample statistics are not available to fully describe those patterns. However, as will be discussed in section 5.3, vertical connections among different channel deposits produce long-range connectivity of the high-k material. [41] In addition to the ergodicity issue, differences between the training image and output realizations are also a consequence of the stochastic, pixel-based multiple-point

8 of 16


W08426

Table 1. Hydraulic Property Values Hydrostratigraphic Unit HSU-2 HSU-3A HSU-2/HSU-3A

1

Facies

K (m/d)

Ss (m )

channel floodplain channel floodplain paleosol

90 0.3 72 0.24 6 10 5

5 10 6 1.8 10 4 4 10 6 1.4 10 4 1.8 10 4

algorithm that we use. Some variability in geometrical properties (beyond the variability contained in the training image) is introduced by the geostatistical simulator. This issue has been documented elsewhere [e.g., Caers and Zhang, 2004]. Most importantly, the output conditional realizations are realistic for the depositional environment that we consider, consistent with the local geologic data, and capture many of the key properties specified in the training image. [42] Also represented in Figures 6a and 6b is the contact between the two hydrostratigraphic units that comprise our modeled section. Previous work indicates that this contact is most likely formed by an extensive paleosol deposit. Using existing data and interpretations regarding its elevation at borehole locations, we contoured the paleosol surface. We incorporate this surface into our model by assigning the lowest k value to all grid cells that it intersects. However, at locations between the boreholes, we permit discrete erosional breaks in the paleosol surface, and we simulate those breaks stochastically. The aquifer test data support the notion of a discontinuous paleosol surface; observation points in HSU-3A respond relatively quickly to pumping in HSU-2. Therefore, at least in this particular area, there appears to be some hydraulic communication between the two HSUs, beyond the typical slow leakage response that would be observed with a competent and continuous low-k unit. In the model, we permit this communication by allowing breaks in the bounding paleosol. Our approach is simply to specify a break (i.e., assign high-k channel rather than low-k paleosol material) at any location where a simulated channel deposit coincides with the (contoured) paleosol surface. 4.2. Groundwater Flow Model [43] We apply the finite difference code MODFLOW2000 [Harbaugh et al., 2000] with a multigrid matrix solver [Wilson and Naff, 2004] to discretize and numerically solve (1) for our problem domain. The grid spacing used for the groundwater model is identical to that of the geostatistical simulation grid. Thus, no upscaling or downscaling step is required prior to running flow simulations. [44] The modeled vertical section ranges from approximately 28 to 32 m in thickness and spans the two hydrostratigraphic units where pumping-induced hydraulic head declines were observed. Top and bottom model boundaries follow low-permeability paleosol deposits and are specified as no flow (@h(x)/@nj0 = 0 for x on the boundary segment). Boundary conditions prescribed along the model sides are shown on the plan view map in Figure 1b. Inflow and outflow segments are treated using Cauchy boundary conditions; the computed flux into or out of the modeled

W08426

region depends on the difference between the boundary head (permitted to vary) and a fixed head located outside the domain. Also shown in Figure 1b are two no-flow boundary segments that coincide with regional flow lines. For the aquifer test conditions and hydraulic properties evaluated here, these external boundaries are sufficiently far from the pumping location and therefore have a negligible impact on simulated drawdowns. [45] For each realization of the subsurface architecture, we use the simulated facies distribution as a template and assign hydraulic property values accordingly. The values of hydraulic conductivity and specific storage for each material type are given in Table 1. The tabulated property values were obtained primarily through model calibration, although we used relevant field data to guide the calibration process. A discussion of these relevant data (field estimated K values), including some comparison to our calibrated conductivity values, is provided with the auxiliary material. We emphasize that the property values listed in Table 1 are well within the range of published values for each material texture [e.g., Freeze and Cherry, 1979]. [46] We apply our groundwater model to simulate the 24-h aquifer test described in section 2. Computed drawdowns are compared to measured drawdowns at the six observation well locations. Since we use a very large grid, deal with strong permeability contrasts, and run the transient flow model within a broader simulation inversion framework, computational efficiency is a key consideration. All simulation results presented in this paper were generated on a shared-memory Linux workstation. On this computing platform, both the geostatistical simulations and the forward flow simulations access over 1 GB of physical memory. The memory requirement of the geostatistical model (nearly 4 GB) is significantly higher than that of the flow simulation model. The run time to generate a single geostatistical realization of the facies distribution is approximately 20 min, while the flow simulation run time is typically in the range of 50 to 60 min. Therefore, each function evaluation (iteration) directed by the inverse model requires about 75 min. The inverse model performance is discussed in section 5.2.

5. Results [47] In this section, flow modeling results are presented for five different realizations of the subsurface architecture. Each realization was obtained using the inverse procedure described in section 3.1. 5.1. Flow Model Response [48] For all realizations, the dynamic model indicates that, in response to pumping, groundwater flow becomes focused within the high-k channel material. This is consistent with previous modeling results reported by Carle [1996] and Fogg et al. [2000]. Additionally, we find that the drawdown distribution is highly sensitive to the locations and geometry of specific channel deposits. This sensitivity enables the inversion approach that we use to identify individual channels. Because the hydraulic response behavior is highly sensitive, dynamic observations are expected to have good information content that can be exploited to resolve the channel structure.

9 of 16

W08426


Figure 7. Early time propagation of drawdown in vicinity of pumped well. (a) Facies distribution; channel deposits are shown in blue. (b) Simulated drawdown after 10 min of pumping. (c) Simulated drawdown after 60 min.

10 of 16

W08426

W08426


W08426

Figure 8. Comparison of simulated and measured drawdowns at observation wells. Open circles are the observed data. Lines represent the simulated response for a given realization; M1 through M5 refer to five distinct realizations. [49] Figure 7 shows the simulated drawdown distribution during the first hour of pumping for one particular realization. As revealed in Figure 7, early time drawdown closely tracks a large channel body that is intersected by the pumping well screen. The channel structure preferentially induces drawdowns to the northwest and southeast, an effect that facilitates the rapid response observed at a distant well, W264. From Figure 7, it is clear that the drawdown propagation rate within the channel material is higher than the rate for the floodplain matrix. This is due to the contrasting hydraulic diffusivities of each facies. As a result of this behavior, the response timing at a given observation point is a function not only of its distance, but also its connectedness (via the network of channel deposits) with the pumped area. 5.2. Performance of the Inverse Model [50] In Figure 8, model-simulated drawdown is compared to the measured drawdown at observation well locations. On each individual well plot, simulated results are presented

for all five inverse solutions. The model generally performs better at the HSU-2 observation wells. Difficulty in reproducing drawdown observations for the deeper wells may be partly related to uncertainties in the data. As discussed in the auxiliary material, within HSU-3A there was a trend of declining hydraulic heads prior to the start of the aquifer test. Although we used available measurements to correct for this trend at W712 and W267, the pretest variability is still a complicating factor that may introduce some (data) errors. [51] We use an overall sum of relative errors as our objective function to direct the inverse model Obj ¼

Nwell TO XX w¼1 t¼1

w ssim ðt Þ swobs ðt Þ swobs ðt Þ

ð10Þ

where swobs (t) is the observed drawdown for well w at time t, swsim (t) is the model-simulated drawdown for well w at time

11 of 16

W08426


Figure 9. Objective function value at each iteration for an example inverse modeling run. Circles identify iterations where there was an improvement over the previous best realization. The simulated facies distribution corresponding to the final iteration was retained as an inverse solution; this particular realization of the subsurface architecture produced an acceptable match to the drawdown measurements (20% error on average). t, TO is the total number of observation times (for a given well), and Nwell is the total number of wells. Application of a relative error measure is important since the smaller drawdowns (measured at early times) are often most informative about the subsurface variability. As demonstrated in Figure 7, dynamic head observations that immediately follow an applied stress can be quite sensitive to detailed heterogeneity and connectivity properties. In comparison, late-time measurements depend on larger-scale average properties [Meier et al., 1998]. [52] In Figure 9, we show the objective function (as the mean of the relative errors) at each iteration of the inverse model for an example run. The results indicate a relatively high number of nonproductive realizations. This reflects the discontinuous nature of the problem and the large permeability contrasts that are modeled. Because of the discontinuity (i.e., discrete facies framework), the application of a traditional gradient-based inverse technique is not feasible. Research is underway to explore new approaches for enhancing the efficiency of the probability perturbation method [e.g., Johansen et al., 2007]. Most important for the present work, however, is that the results in Figure 9 reveal a downward trend, showing clear progress toward an inverse solution. Through Monte Carlo sampling from (6), we attempted to find a realization that adequately reproduced the drawdown observations without performing the inversion. This effort was not successful, which indicates that the inverse modeling is necessary to obtain a good solution. Caers and Hoffman [2006] show that the sampling properties of the probability perturbation method (PPM) are similar to those of a classic rejection sampler; however, PPM is far more efficient. [53] Figure 6c further demonstrates the efficacy of the inverse model. Realizations sampled from (6), which are conditioned to the direct d1 data but not the dynamic data, consistently underpredict both the timing and magnitude of the measured drawdowns at well W264. It requires a specific high-k connection, identified by the inverse model, to reproduce the data at this location. This is true despite the fact that all realizations have the same bulk anisotropy behavior (resulting from similar channel fraction, size,

W08426

orientation, etc.). In section 5.3 below, we characterize the specific feature that explains the response behavior at this observation well. [54] During the search for a solution, the facies distribution is modified significantly. In comparing the simulated architecture for the first and last iterations, approximately 30% of the grid cells have different facies values. [55] To further evaluate the model performance at each observation well, we consider three different statistical measures, including the mean absolute error (MAE), the mean relative error (MRE), and the Nash-Sutcliffe coefficient of efficiency (NSE) [Nash and Sutcliffe, 1970] TO P

NSE ¼ 1

ðsobs ðt Þ ssim ðtÞÞ2

t¼1 TO P

ð11Þ ðsobs ðt Þ sobs Þ2

t¼1

Values of MAE, MRE, and NSE for each observation well are provided in Table 2. The tabulated results show that relative errors are generally below 25% except at well W267. As noted previously, the larger errors at this location may be due to data uncertainty introduced by an external forcing function. 5.3. Connectivity Properties and Shortest-Path Results [56] We now focus on the structural properties that give rise to the complex flow behavior discussed above. We perform cluster counting [Stauffer and Aharony, 1994; Pardo-Igu´zquiza and Dowd, 2003] to identify sets of connected channel cells for each realization of the facies distribution. Two locations are connected (i.e., on the same cluster) if there exists a continuous path of channel material between them. The probability of finding two connected channel locations decreases with distance and, for anisotropic models such as ours, is also a function of direction. For the conditional realizations that we generate, the cluster counting routinely finds a spanning cluster. A spanning (or percolating) cluster is a set of connected channel cells that spans the simulation domain in each major grid direction. Previous investigators at LLNL also observed a spanning cluster when channel deposits were generated using the transition probability Markov chain method [Fogg et al., 2000]. The transition probability approach, while capable of capturing important cross correlation among different facies, is based on a two-point spatial statistical measure. The multiple-point statistics approach allows for more

Table 2. Statistical Measures of Model Performance at All Observation Wellsa Observation Well

MAE (m)

MRE

NSE

W614 W714 W264 W616 W712 W267

0.0058 0.0028 0.0080 0.0190 0.0086 0.0149

0.22 0.16 0.13 0.25 0.22 0.38

0.63 0.91 0.84 0.57 0.78 0.24

a Reported values for each well are obtained by averaging over the five realizations shown in Figure 8.

12 of 16

W08426


Figure 10. Main channel body (from example realization) that connects the pumping well and distant observation well W264. explicit representation of channel geometry and connectivity. Given the fundamental differences in methodology, models generated with each technique are expected to have different percolation thresholds. This is a consequence of important differences in the simulated channel geometry. The multiple-point geostatistical model used here produces more of a curvilinear cross-sectional channel geometry and greater channel sinuosity compared to any previous studies. Both models, parameterized appropriately for LLNL, have a percolation threshold (pc) that is substantially lower than the pc value of 0.31 for standard (completely random) site percolation on a cubic lattice. For example, it appears that pc 0.18 (the site-specific channel fraction) for both models. In our case, given that we typically identify a very high fraction of channels (>80%) on the spanning cluster, the percolation threshold may be significantly less than 0.18. However, the effects of finite domain size come into play, making it difficult to identify a definitive pc value for any particular structural model [Harter, 2005]. [57] Like the transition probability – based model, our results indicate that vertical connections created by stacked or overlapping channel deposits are key in producing longrange connectivity (both vertical and horizontal) of the high-k material. A similar finding was reported by Proce et al. [2004], who investigated the three-dimensional connectivity of sand and gravel facies in glaciofluvial deposits. [58] Our inverse solutions consistently identify a set of vertically stacked channel deposits that provide a direct connection between the pumped interval and observation well W264. Figure 10 depicts this channel structure, which is part of a larger channel cluster, for the example realization presented in Figure 6a. Although in Figure 6a the two wells appear to intersect the same individual channel, their screened intervals are offset by over 3 m (after correcting for the regional stratigraphic dip). Therefore, the high-k connection between these locations requires vertically intersecting channel deposits. [59] We perform shortest-path analysis to further investigate the significance of the connectivity event involving wells W612 and W264. Using coordinates for each well screen, we define a 3-D vector h that describes the separation of these locations. khk = 206.7 m is the Euclidean distance between the two observation points. An important quantity is l, the length of the shortest path between points on the same channel cluster separated by h. The only constraint is that all grid cells along the path must be occupied by channel material. Clearly, l khk. We use Dijkstra’s algorithm to identify the shortest among all connected paths. For the W612$W264 connection in particular, we find that l assumes its lowest possible value (for the structural grid that we employ) in four of the five

W08426

calibrated models. In the fifth model, the shortest path is only slightly ( 200 m) observation well, where rapid hydraulic response was recorded. The structural feature that offers this direct path is a set of vertically stacked channel deposits. [64] The training image-based geostatistical approach used in this study requires information on the relative abundance and geometry of each modeled facies. For the depositional environment that we considered, key sitespecific data requirements included the width, thickness, and sinuosity of fluvial channel deposits. Although we relied on a large data set and a significant amount of prior geologic interpretation, there remains (as always) some uncertainty in the specified properties. Investigations like the one presented here, which attempt to understand the influence of structural heterogeneity, will benefit from new research developments on methods to characterize the geometry of subsurface structures. Regarding fluvial deposits in particular, modern channels offer important guidance. We echo the conclusions of Rubin et al. [2006] on the need to better understand the relationship between the geometry of surficial river channels and the geometry of channel deposits. [65] Aquifer tests have been and will likely remain an important tool in hydrogeologic investigations. Traditional interpretation of aquifer test data provides estimates of average hydraulic properties. However, as demonstrated in this work and elsewhere in the literature [Raghavan, 2004, and references therein], when combined with geologic data, a critical analysis of aquifer test results can yield additional insight into important subsurface details. A key research item moving forward will be to address the information content in anomalous response curves. For example, for a particular location and set of observations, what does such response behavior suggest about the presence and/or distribution of subsurface structures? Finally, we emphasize the importance of new test design and interpretation methods. Hydraulic tomography and similarly motivated techniques may significantly advance efforts to identify key structural features that control subsurface dynamics.

W08426

[66] Acknowledgments. Funding for this work was provided by a National Defense Science and Engineering graduate fellowship awarded to M. Ronayne. Additional support was provided by NSF grant EAR0537668. Any opinions, findings, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. We thank the Lawrence Livermore National Laboratory Environmental Restoration Department for providing the field data. Computing resources for this study were provided by the Stanford Center for Computational Earth and Environmental Science. We acknowledge the associate editor and three anonymous reviewers for constructive review comments that improved this paper.

References Andrade, J. S., Jr., S. V. Buldyrev, N. Dokholyan, S. Havlin, P. R. King, Y. Lee, G. Paul, and H. E. Stanley (2000), Flow between two sites on a percolation cluster, Phys. Rev. E, 62, 8270 – 8281, doi:10.1103/ PhysRevE.62.8270. Arau´jo, A. D., A. A. Moreira, H. A. Makse, H. E. Stanley, and J. S. Andrade Jr. (2002), Traveling length and minimal traveling time for flow through percolation networks with long-range spatial correlations, Phys. Rev. E, 66, 046304, doi:10.1103/PhysRevE.66.046304. Arpat, G. B., and J. Caers (2007), Conditional simulation with patterns, Math. Geol., 39, 177 – 203, doi:10.1007/s11004-006-9075-3. Baker, D. R., G. Paul, S. Sreenivasan, and H. E. Stanley (2002), Continuum percolation threshold for interpenetrating squares and cubes, Phys. Rev. E, 66, 046136, doi:10.1103/PhysRevE.66.046136. Bennett, G. L., G. S. Weissmann, G. S. Baker, and D. W. Hyndman (2006), Regional-scale assessment of a sequence-bounding paleosol on fluvial fans using ground-penetrating radar, eastern San Joaquin Valley, California, GSA Bull., 118, 724 – 732, doi:10.1130/B25774.1. Bunde, A., and S. Havlin (1996), Percolation I, in Fractals and Disordered Systems, 2nd ed., edited by A. Bunde and S. Havlin, pp. 59 – 113, Springer, New York. Caers, J. (2003), History matching under training-image based geological model constraints, SPE J., 8, 218 – 226, doi:10.2118/74716-PA. Caers, J., and T. Hoffman (2006), The probability perturbation method: A new look at Bayesian inverse modeling, Math. Geol., 38, 81 – 100, doi:10.1007/s11004-005-9005-9. Caers, J., and T. Zhang (2004), Multiple-point geostatistics: A quantitative vehicle for integrating geologic analogs into multiple reservoir models, AAPG Mem., 80, 383 – 394. Carle, S. F. (1996), A transition probability-based approach to geostatistical characterization of hydrostratigraphic architecture, Ph.D. dissertation, Univ. of Calif., Davis. Carle, S. F., and G. E. Fogg (1996), Transition probability-based indicator geostatistics, Math. Geol., 28, 453 – 476, doi:10.1007/BF02083656. Carle, S. F., and G. E. Fogg (1997), Modeling spatial variability with one and multidimensional continuous-lag Markov chains, Math. Geol., 29, 891 – 902, doi:10.1023/A:1022303706942. Cortis, A., and C. Knudby (2006), A continuous time random walk approach to transient flow in heterogeneous porous media, Water Resour. Res., 42, W10201, doi:10.1029/2006WR005227. de Marsily, G., F. Delay, J. Gonçalve`s, P. Renard, V. Teles, and S. Violette (2005), Dealing with spatial heterogeneity, Hydrogeol. J., 13, 161 – 183, doi:10.1007/s10040-004-0432-3. Deutsch, C. V., and A. G. Journel (1998), GSLIB: Geostatistical Software Library and User’s Guide, 2nd ed., 369 pp., Oxford Univ. Press, New York. Deutsch, C. V., and T. T. Tran (2002), FLUVSIM: A program for objectbased stochastic modeling of fluvial depositional systems, Comput. Geosci., 28, 525 – 535, doi:10.1016/S0098-3004(01)00075-9. Dokholyan, N. V., Y. Lee, S. V. Buldyrev, S. Havlin, P. R. King, and H. E. Stanley (1998), Scaling of the distribution of shortest paths in percolation, J. Stat. Phys., 93, 603 – 613, doi:10.1023/B:JOSS.0000033244. 13545.da. Feyen, L., and J. Caers (2006), Quantifying geological uncertainty for flow and transport modeling in multi-modal heterogeneous formations, Adv. Water Resour., 29, 912 – 929, doi:10.1016/j.advwatres.2005.08.002. Fogg, G. E. (1986), Groundwater flow and sand body interconnectedness in a thick, multiple-aquifer system, Water Resour. Res., 22, 679 – 694, doi:10.1029/WR022i005p00679. Fogg, G. E., C. D. Noyes, and S. F. Carle (1998), Geologically based model of heterogeneous hydraulic conductivity in an alluvial setting, Hydrogeol. J., 6, 131 – 143, doi:10.1007/s100400050139. Fogg, G. E., S. F. Carle, and C. Green (2000), A connected-network paradigm for the alluvial aquifer system, in Theory, Modeling, and Field Investigation in Hydrogeology: A Special Volume in Honor of Shlomo

14 of 16

W08426


P. Neuman’s 60th Birthday, edited by D. Zhang and C. L. Winter, Spec. Pap. Geol. Soc. Am., 348, 25 – 42. Freeze, R. A., and J. A. Cherry (1979), Groundwater, 604 pp., PrenticeHall, Englewood Cliffs, N. J. Garboczi, E. J., K. A. Snyder, J. F. Douglas, and M. F. Thorpe (1995), Geometrical percolation threshold of overlapping ellipsoids, Phys. Rev. E, 52, 819 – 828, doi:10.1103/PhysRevE.52.819. Gibling, M. R. (2006), Width and thickness of fluvial channel bodies and valley fills in the geological record: A literature compilation and classification, J. Sediment. Res., 76, 731 – 770, doi:10.2110/jsr.2006.060. Harbaugh, A. W., E. R. Banta, M. C. Hill, and M. G. McDonald (2000), MODFLOW-2000, the U.S. Geological Survey modular ground-water model—User guide to modularization concepts and the ground-water flow process, U.S. Geol. Surv. Open File Rep., 00-92, 121 pp. Harter, T. (2005), Finite-size scaling analysis of percolation in threedimensional correlated binary Markov chain random fields, Phys. Rev. E, 72, 026120, doi:10.1103/PhysRevE.72.026120. Herrmann, H. J., and H. E. Stanley (1988), The fractal dimension of the minimum path in two- and three-dimensional percolation, J. Phys. A Math. Gen., 21, L829 – L833, doi:10.1088/0305-4470/21/17/003. Hoffman, B. T. (2005), Geologically consistent history matching while perturbing facies, Ph.D. dissertation, Stanford Univ., Stanford, Calif. Hoffman, B. T., and J. Caers (2005), Regional probability perturbations for history matching, J. Pet. Sci. Eng., 46, 53 – 71, doi:10.1016/j.petrol. 2004.11.001. Hoffman, F., and M. D. Dresen (1990), A method to evaluate the vertical distribution of VOCs in ground water in a single borehole, Ground Water Monit. Rev., 10, 95 – 100, doi:10.1111/j.1745-6592.1990.tb00342.x. Hoffman, F., R. G. Blake, R. J. Gelinas, C. D. Noyes, Z. Demir, and P. F. McKereghan (2003), A conceptual model and remediation strategy for volatile organic compounds in groundwater in unconsolidated sediments: A Lawrence Livermore National Laboratory case study, Environ. Eng. Geosci., 9, 83 – 94, doi:10.2113/9.1.83. Hunt, A. G. (2001), Applications of percolation theory to porous media with distributed local conductances, Adv. Water Resour., 24, 279 – 307, doi:10.1016/S0309-1708(00)00058-0. Johansen, K., J. Caers, and S. Suzuki (2007), Hybridization of the probability perturbation method with gradient information, Comput. Geosci., 11, 319 – 331, doi:10.1007/s10596-007-9055-9. Johnson, M. E. (1987), Multivariate Statistical Simulation, 230 pp., John Wiley, New York. Johnson, N. M., and S. J. Dreiss (1989), Hydrostratigraphic interpretation using indicator geostatistics, Water Resour. Res., 25, 2501 – 2510, doi:10.1029/WR025i012p02501. Journel, A. G. (2002), Combining knowledge from diverse sources: An alternative to traditional data independence hypotheses, Math. Geol., 34, 573 – 596, doi:10.1023/A:1016047012594. Knudby, C., and J. Carrera (2005), On the relationship between indicators of geostatistical, flow, and transport connectivity, Adv. Water Resour., 28, 405 – 421, doi:10.1016/j.advwatres.2004.09.001. Knudby, C., and J. Carrera (2006), On the use of apparent hydraulic diffusivity as an indicator of connectivity, J. Hydrol., 329, 377 – 389, doi:10.1016/j.jhydrol.2006.02.026. Lee, S.-Y., S. F. Carle, and G. E. Fogg (2007), Geologic heterogeneity and a comparison of two geostatistical models: Sequential Gaussian and transition probability-based geostatistical simulation, Adv. Water Resour., 30, 1914 – 1932, doi:10.1016/j.advwatres.2007.03.005. Lee, Y., J. S. Andrade Jr., S. V. Buldyrev, N. V. Dokholyan, S. Havlin, P. R. King, G. Paul, and H. E. Stanley (1999), Traveling time and traveling length in critical percolation clusters, Phys. Rev. E, 60, 3425 – 3428, doi:10.1103/PhysRevE.60.3425. Leven, C., and P. Dietrich (2006), What information can we get from pumping tests?—Comparing pumping test configurations using sensitivity coefficients, J. Hydrol., 319, 199 – 215, doi:10.1016/j.jhydrol.2005.06. 030. Liu, G., C. Zheng, and S. M. Gorelick (2004), Limits of applicability of the advection-dispersion model in aquifers containing connected high-conductivity channels, Water Resour. Res., 40, W08308, doi:10.1029/ 2003WR002735. Liu, Y., A. Harding, W. Abriel, and S. Strebelle (2004), Multiple-point simulation integrating wells, three-dimensional seismic data, and geology, AAPG Bull., 88, 905 – 921, doi:10.1306/02170403078. Mansoor, K., M. Maley, Z. Demir, and F. Hoffman (2002), Implementation of deterministically-derived hydrostratigraphic units into a 3D finite element model at the Lawrence Livermore National Laboratory Superfund site, paper presented at Bridging the Gap Between Measurement and

W08426

Modeling in Heterogeneous Media: International Groundwater Symposium, Int. Assoc. of Hydraul. Eng. and Res., Berkeley, Calif. Matthaï, S. K., and M. Belayneh (2004), Fluid flow partitioning between fractures and a permeable rock matrix, Geophys. Res. Lett., 31, L07602, doi:10.1029/2003GL019027. Meier, P. M., J. Carrera, and X. Sańchez-Vila (1998), An evaluation of Jacob’s method for the interpretation of pumping tests in heterogeneous formations, Water Resour. Res., 34, 1011 – 1025, doi:10.1029/ 98WR00008. Mosegaard, K., and A. Tarantola (1995), Monte Carlo sampling of solutions to inverse problems, J. Geophys. Res., 100, 12,431 – 12,447, doi:10.1029/94JB03097. Nash, J. E., and J. V. Sutcliffe (1970), River flow forecasting through conceptual models, part I: A discussion of principles, J. Hydrol., 10, 282 – 290, doi:10.1016/0022-1694(70)90255-6. Noyes, C. D. (1990), Hydrostratigraphic analysis of the Pilot Remediation Test Area, LLNL, Livermore, California, M.S. thesis, Univ. of Calif., Davis. Okabe, H., and M. J. Blunt (2005), Pore space reconstruction using multiple-point statistics, J. Pet. Sci. Eng., 46, 121 – 137, doi:10.1016/ j.petrol.2004.08.002. Omre, H., and H. Tjelmeland (1996), Petroleum geostatistics, in Proceedings of the Fifth International Geostatistics Congress, vol. 1, edited by E. Y. Baafi and N. A. Schofield, pp. 41 – 52, Kluwer Acad., Dordrecht, Netherlands. Pardo-Igu´zquiza, E., and P. A. Dowd (2003), CONNEC3D: A computer program for connectivity analysis of 3D random set models, Comput. Geosci., 29, 775 – 785, doi:10.1016/S0098-3004(03)00028-1. Paul, G., S. Havlin, and H. E. Stanley (2002), Fractal behavior of the shortest path between two lines in percolation systems, Phys. Rev. E, 65, 066105, doi:10.1103/PhysRevE.65.066105. Proce, C. J., R. W. Ritzi, D. F. Dominic, and Z. Dai (2004), Modeling multiscale heterogeneity and aquifer interconnectivity, Ground Water, 42, 658 – 670, doi:10.1111/j.1745-6584.2004.tb02720.x. Raghavan, R. (2004), A review of applications to constrain pumping test responses to improve on geological description and uncertainty, Rev. Geophys., 42, RG4001, doi:10.1029/2003RG000142. Ronayne, M. J., and S. M. Gorelick (2006), Effective permeability of porous media containing branching channel networks, Phys. Rev. E, 73, 026305, doi:10.1103/PhysRevE.73.026305. Rubin, Y., I. A. Lunt, and J. S. Bridge (2006), Spatial variability in river sediments and its link with river channel geometry, Water Resour. Res., 42, W06D16, doi:10.1029/2005WR004853. Sahimi, M. (1993), Flow phenomena in rocks: From continuum models to fractals, percolation, cellular automata, and simulated annealing, Rev. Mod. Phys., 65, 1393 – 1534, doi:10.1103/RevModPhys.65.1393. Schulz, K., R. Seppelt, E. Zehe, H. J. Vogel, and S. Attinger (2006), Importance of spatial structures in advancing hydrologic sciences, Water Resour. Res., 42, W03S03, doi:10.1029/2005WR004301. Stauffer, D., and A. Aharony (1994), Introduction to Percolation Theory, 2nd ed., 181 pp., Taylor and Francis, Philadelphia, Pa. Strebelle, S. (2000), Sequential simulation drawing structures from training images, Ph.D. dissertation, Stanford Univ., Stanford, Calif. Strebelle, S. (2002), Conditional simulation of complex geological structures using multiple-point statistics, Math. Geol., 34, 1 – 21, doi:10.1023/ A:1014009426274. Strebelle, S., K. Payrazyan, and J. Caers (2003), Modeling of a deepwater turbidite reservoir conditional to seismic data using principal component analysis and multiple-point geostatistics, SPE J., 8, 227 – 235, doi:10.2118/85962-PA. Streltsova, T. D. (1988), Well Testing in Heterogeneous Formations, 413 pp., John Wiley, New York. Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation, 342 pp., Soc. for Ind. and Appl. Math., Philadelphia, Pa. Taylor, W. L., D. D. Pollard, and A. Aydin (1999), Fluid flow in discrete joint sets: Field observations and numerical simulations, J. Geophys. Res., 104, 28,983 – 29,006, doi:10.1029/1999JB900179. Theis, C. V. (1935), The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, Eos Trans. AGU, 2, 519. Torquato, S. (2002), Random Heterogeneous Materials: Microstructure and Macroscopic Properties, 701 pp., Springer, New York. Weber, K. J., and L. C. van Geuns (1990), Framework for constructing clastic reservoir simulation models, JPT J. Pet. Technol., 42, 1248 – 1297, doi:10.2118/19582-PA.

15 of 16

W08426


Weissmann, G. S., Z. Yong, G. E. Fogg, R. G. Blake, C. D. Noyes, and M. Maley (2002), Modeling alluvial fan aquifer heterogeneity at multiple scales through stratigraphic assessment, paper presented at Bridging the Gap Between Measurement and Modeling in Heterogeneous Media: International Groundwater Symposium, Int. Assoc. of Hydraul. Eng. and Res., Berkeley, Calif. Western, A. W., G. Blo¨schl, and R. B. Grayson (2001), Toward capturing hydrologically significant connectivity in spatial patterns, Water Resour. Res., 37, 83 – 97, doi:10.1029/2000WR900241. Wilson, J. D., and R. L. Naff (2004), The U.S. Geological Survey modular ground-water model—GMG linear equation solver package documentation, U.S. Geol. Surv. Open File Rep., 2004-1261, 47 pp. Zhang, T. F., P. Switzer, and A. Journel (2006), Filter-based classification of training image patterns for spatial simulation, Math. Geol., 38, 63 – 80, doi:10.1007/s11004-005-9004-x.

W08426

Zheng, C., and S. M. Gorelick (2003), Analysis of solute transport in flow fields influenced by preferential flow paths at the decimeter scale, Ground Water, 41, 142 – 155, doi:10.1111/j.1745-6584.2003.tb02578.x. Zinn, B., and C. F. Harvey (2003), When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate Gaussian hydraulic conductivity fields, Water Resour. Res., 39(3), 1051 – 1070, 1051, doi:10.1029/ 2001WR001146.

J. Caers, Department of Energy Resources Engineering, Stanford University, Stanford, CA 94305-2220, USA. S. M. Gorelick, Department of Environmental Earth System Science, Stanford University, Stanford, CA 94305-2115, USA. M. J. Ronayne, Department of Geological and Environmental Sciences, Stanford University, Stanford, CA 94305-2115, USA. (mronayne@stanford. edu)

16 of 16