Sparse geologic dictionaries for subsurface flow model calibration: Part II. Robustness to uncertainty

Advances in Water Resources 39 (2012) 122–136

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Sparse geologic dictionaries for subsurface flow model calibration: Part II. Robustness to uncertainty Mohammadreza Mohammad Khaninezhad a, Behnam Jafarpour a,b,⇑, Lianlin Li a a b

Department of Petroleum Engineering, Texas A&M University, College Station, TX, United States Department of Chemical Engineering and Material Science, University of Southern California, United States

a r t i c l e

i n f o

Article history: Available online 25 November 2011 Keywords: Geologic dictionaries Sparse reconstruction Model calibration Geologic continuity Prior uncertainty lp-Norm inversion

a b s t r a c t We introduced sparse geologic dictionaries for formulating and solving subsurface flow model calibration problems in Part I. A key assumption was the availability of a reliable prior model library for construction of a relevant sparse geologic dictionary that was used to constrain the inversion solution. In practice, however, the spatial connectivity in rock physical property distributions has to be inferred from uncertain and incomplete sources of information, including qualitative geologic interpretations, formation type and outcrop maps, as well as scattered measurements (e.g., well core, log data, and seismic maps). Thus, the geologic continuity model that is used for generating prior models is likely to carry significant uncertainty. In Part II, we investigate the performance of the proposed method under uncertain and incorrect prior continuity models. We show that, unlike the conventional reduced-order parameterization methods such as truncated SVD, diverse geologic dictionaries are robust against uncertainty in the structural prior model. This robustness is attributed to the selection property of the sparse reconstruction algorithm and the connectivity exhibited by the dictionary elements. Under highly uncertain prior models, the reconstruction problem is reduced to identifying and combining relevant elements from a large and diverse geologic dictionary. Diversity of the dictionary further enhances solution sparsity since a large number of dictionary elements will have negligible contribution to the solution. Consequently, the inversion provides considerable flexibility to accommodate significant levels of variability (uncertainty) in prior geologic models. As an extreme case, we also evaluate the performance of the proposed method when even the formation type (e.g., Gaussian versus non-Gaussian) is uncertain. In addition, we examine the sensitivity of the proposed method to noise in the flow data, where we observe solution underestimation effects due to l1-norm approximation of l0-norm regularization. We address this issues by replacing l1-norm with an adaptive lp-norm, where p decreases with iterations from an initial value p = 1 to a value of p = 0+ at final iterations. Using several numerical examples, we illustrate these important properties of the proposed inversion approach and compare its performance with that of the truncated singular vector parameterization. Ó 2011 Published by Elsevier Ltd.

1. Introduction Data limitation in characterization of subsurface environments results in extensive interpretation and interpolation [14]. As a consequence, significant uncertainty is introduced into constructing static prior models. Calibration of prior models against dynamic data is a routine task to reduce the existing uncertainties in static prior models (e.g., [10,12,17–19,34,44,45]). While dynamic data provide important constraints for updating static models, the integral and scarce nature of these measurements lead to an information ⇑ Corresponding author at: Mork Family Department of Chemical Engineering and Material Science, Viterbi School of Engineering, University of Southern California, 925 Bloom Walk, HED 313, Los Angeles, CA 90089-1211, United States. Tel.: +1 213 740 2228. E-mail address: [email protected] (B. Jafarpour). 0309-1708/$ - see front matter Ó 2011 Published by Elsevier Ltd. doi:10.1016/j.advwatres.2011.10.005

content that does not support the identification of subsurface heterogeneity at an adequate resolution to predict the multi-scale nature of the flow and transport behavior [42,43]. Overcoming the challenges associated with subsurface characterization requires important breakthroughs in data acquisition technology for highresolution imaging as well as more sophisticated data fusion algorithms to effectively integrate various sources of data [42,43]. Physical modeling and computational simulation of geologic processes can also provide important constraints by generating formation specific physics-based prior models [16]. Such prior models can be extremely important in guiding stochastic simulation of complex (non-Gaussian) spatial connectivity patterns that are responsible for extreme flow behavior [16,27]. Parallel to advancing monitoring data collection systems, mathematical and computational approaches are needed to combine complex and highly uncertain prior models with diverse and disparate datasets to

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

provide a cohesive understanding of the dynamic flow behavior in complex geo-systems [12,42,43]. Many realistic geologic systems are complex and have local spatial connectivity structures that are not amenable to conventional second-order geostatistical descriptions of heterogeneity [8,9,13, 14,40,41,46]. Representing nonlinear, non-Gaussian geologic features that determine important preferential flow and transport patterns, such as meandering channels or thin shale lenses, requires more sophisticated higher-order statistical simulation algorithms such as multipoint geostatistics (MPS) [6,40]. The MPS methods perform the simulation at the grid level, which facilitates the integration of hard and soft static data. However, in addition to simulation of complex geologic patterns that are conditioned on various static data, MPS-based methods must be developed for integration of dynamic data to constrain the simulated models to reproduce the observed dynamic response in the field [7,23,38]. The importance of prior model cannot be neglected in solving dynamic inverse problems for subsurface characterization. Reliable prior models are critical for regularizing model calibration inverse problems, which are usually ill-posed because of the lack of sufficient constraining data. When complex geologic formations are considered, modeling the geologic processes that create the relevant patterns and structural connectivity can be used to generate the conceptual models needed to confidently constrain the plausibility of model calibration solutions. Just as in medical imaging where general knowledge about the organ of interest significantly constrains the type of solutions that can be deemed acceptable, even qualitative information about the type of geologic formations and depositional environment can also provide substantial constraining power in solving subsurface characterization inverse problems. However, an important difference in geosciences applications is that geologic models can be widely variable and subject to various sources and levels of uncertainty. Accounting for this uncertainty is of paramount importance in avoiding biased model calibration solutions that reproduce the existing data but fail to predict the future fluid flow and transport performance of the underlying aquifer/reservoir. For instance, simulation of patterns from a training image in multipoint statistical techniques constrains the generated results to high-order statistics that is encoded in the training image. When the training image is reliable, it provides considerable constraining power about the expected type of spatial variability and structural connectivity. However, conceptual geologic models such as a training image are generated from interpretation of incomplete and uncertain sources of information, including, but not limited to, the geologic history of the region, outcrop data, well log data and 3D seismic images. Consequently, one needs to be mindful of (and account for) the uncertainty that can exist in an adopted conceptual model. The success of model calibration process, in this case, becomes highly dependent on the credibility of the prior geologic model (i.e., training image) used [23,26]. Given the importance of accurate description of connectivity in subsurface flow properties and in particular for characterization of preferential flow paths/barriers, any error or bias in representing such features can lead to significant departures from the actual flow and displacement patterns in the field. Despite their uncertain nature, geologically-informed prior models of complex structures can still enhance model construction, inversion, and systematic uncertainty analysis. A practical approach to uncertainty quantification is to use several distinct and likely conceptual models that cannot be differentiated based on available (static) information. Additional flow and monitoring data can, subsequently, be used to further constrain plausible scenarios. This is the procedure that was adopted in [26], where multiple training images (representing alternative structural assumptions) were used during data integration to account for the uncertainty in the conceptual geologic model. The significance of each individual training image was identified

123

by the prediction performance of the samples generated from it and the dynamic flow data was used to distinguish between consistent and inconsistent training images. In Part II, we show how geologic dictionaries that were introduced in Part I [33] provide an effective framework for incorporation of prior structural uncertainty into the model calibration process. A key factor in our formulation is sparsity, which is the underlying property of recent developments in basis pursuit and compressed sensing paradigm [11,15]. Application of spatial and transform-domain sparsity (compactness) constraint to subsurface modeling and geophysical inverse problems has been considered in the past [2,21,22,28,35,36]. Here, using the sparsity-regularized inversion approach of Part I [33], we show that, interestingly, increasing the diversity in the prior model library implies that most of the elements in the dictionary will have little or no contributions to the reconstruction of the solution, a mechanism that trends to enhance solution sparsity. As opposed to reducedorder parameterization methods such as the truncated SVD, the sparse reconstruction formulation under consideration offers the flexibility of incorporating many diverse elements in the search for the solution, providing a considerably more powerful construction capability. We demonstrate that sparse geologic dictionaries that are learned from very uncertain prior models, in this case using the K-SVD algorithm [1], are quite effective in summarizing the content of the diverse features in the prior model realizations while preserving their distinct and complex spatial connectivity. As we will show in the next section, integration of dynamic data and the selective nature of the sparsity-regularized approach provide an effective and robust framework for identifying the consistent structural connectivity features in the diverse geologic dictionary. We also illustrate an important property of learning dictionaries with the K-SVD algorithm that presents a significant advantage over the truncated SVD approach; namely, continuity preservation. That is, even when highly uncertain prior model realizations are used to cover a wide range of possible features, the KSVD algorithm preserves the continuity in the original model realizations. We show that under the same circumstances, the linear SVD operation tends to disrupt the distinct continuity in the prior model realizations by decomposing them into spectral features that fail to exhibit the original continuity. We begin Part II by evaluating the performance of the K-SVD dictionary construction algorithm under a wide range of prior uncertainty and compare its reconstruction performance with that of the truncated SVD parameterization basis. We also consider the effect of observation error and examine the robustness of the proposed sparse reconstruction formulation in Part I [33] under noisecorrupted measurements. This leads to underestimated solutions related to the approximate nature of l1-norm regularization that is replacing the more complex l0-norm regularization. We present a simple modification of the algorithm by using the lp-norm, where p decreases from 1 to 0+ with the iterations. Finally, we discuss methods for incorporating the uncertainty about the formation type, which represents an extreme but in some cases realistic situation in dealing with uncertainty.

2. Prior uncertainty One approach to account for highly uncertain prior models is to use prior-independent or pre-constructed compressive bases (such as DCT and DWT) that are widely used in image compression [5,31]. Application of these compressive transforms to subsurface characterization has been studied in the literature (e.g., [3,4,20, 24,30,37]). This approach eliminates any reliance on specific prior information while incorporating implicit assumptions about the

124

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

structural attributes of the solution (e.g., smoothness) to reduce parameter dimensionality. Prior-independent parameterization methods have excellent model reduction properties; however, without high resolution and extensive coverage in the data they can have difficulty reproducing the connectivity in complex geologic formations. In general, even a qualitative prior knowledge about the underlying geology of a formation can provide distinctive information about the solution attributes. Hence, incorporation of such knowledge can play an important role in constraining the solution of model calibration inverse problems. In general, information about the depositional environment and geologic history of a region may be used to identify whether the existing features are amenable to second-order Gaussian modeling or higher order statistics must be incorporated in representing the existing connectivity. The knowledge about formation type and the conceptual geologic continuity model itself is of significant value in constraining inversion solutions even though the specific geometric attributes of the structural model may be unknown. For instance, if variogram-based modeling is determined to be appropriate for a particular formation, the variogram model type, its direction of anisotropy and its major correlation lengths can be considered as uncertain parameters. In practice, these parameters are estimated from statistical analysis of low-resolution seismic and scattered hard data using variogram estimation and modeling techniques [14]. In subsurface modeling, once a variogram (or any other structural) model is identified it is common to treat it with certainty. In some cases, the inversion solution is constrained to reproduce the variogram model. Given data scarcity in subsurface systems and the imperfect modeling assumptions that go into the construction of the continuity model, it may be more realistic to treat variogram model parameters with uncertainty. That is, instead of conditioning the solution to honor a prespecified variogram model, one could use the dynamic flow data as additional source of information to identify the parameters of a geologic continuity model. This view toward solving subsurface inverse problems is taken in [25], where the performance of ensemble Kalman filter is evaluated under variogram uncertainty. In [25], it is shown that the EnKF performance was strongly influenced by the assumptions about the prior variogram model parameters. This conclusion is likely to hold for most ill-posed inverse problems where the solution is strongly constrained by the prior model. Taking the parameters of the geologic continuity model as random variables, we generate several model realizations whose spatial distributions represent a wide range of structural variability. We use the dynamic data to characterize the type and form of geologic continuity in addition to identifying the distribution of flow properties. We evaluate the performance of the learned K-SVD and regular SVD algorithms under significant uncertainty in the geometric attributes and orientation of the resulting prior library samples. A desired property of the dictionary is compression of the main information in the library of prior models while preserving the wide range of continuity structures. We consider four sets of experiments as follows: (i) fluvial channel examples with unknown channel orientations; (ii) Gaussian heterogeneous fields with random anisotropy direction and correlation length scales; (iii) top layer of the SPE10 model for which the variogram is not known; and (iv) mixed channel and Gaussian priors to represent uncertainty in the formation type. We also evaluate the sensitivity of the method to measurement noise and modify the algorithm to be more robust against it. 2.1. Robustness against significant prior uncertainty We first examine the robustness of the K-SVD dictionary against prior uncertainty in a fluvial setting with channel features that are

Fig. 1. Training image used to generate prior library of facies models with the snesim MPS algorithm. The training image is adapted from [32].

not amenable to Gaussian descriptions. We use the snesim multipoint geostatistical algorithm to simulate prior fluvial facies realizations from the training image in Fig. 1. Figs. 2a–c show 20 samples from three different prior realizations (left column) and the corresponding K-SVD dictionaries with K = 2025 (second column) and K = 400 (third column) elements. In these examples, the sparsity of the K-SVD dictionary is S = 20. In each case, the rightmost column displays the leading 20 SVD basis elements for the same library. The prior channel facies realizations in Figs. 2a and b have predominantly South–North and East–West orientations, respectively. Since the true model in this example has predominantly South–North directionality, we have labeled the prior models in Figs. 2a and b as correct and incorrect, respectively. As can be verified from Fig. 2, the K-SVD elements tend to preserve the shape of the original facies samples. For K = 400, the dictionary compresses the prior samples such that with a considerably reduced size different S = 20 elements of the dictionary can still provide an accurate sparse description for all 2025 elements in the library. When reducing the dictionary size, the algorithm tends to remove the fine details and retain the main large-scale variabilities in the prior models, which is advantageous for preserving the continuity in channel systems. Fig. 2c displays similar K-SVD dictionary and SVD basis elements for prior samples with random channel orientations in the [0° 180°] range. This range was divided into nine uniform intervals and 400 samples from each interval were generated using the Latin Hypercube Sampling (LHS). Hence, a total of 3600 prior permeability samples were generated and used to construct the dictionary. A comparison between the K-SVD dictionary and SVD basis reveals an important distinction. While the SVD basis components lose the original continuity and orientation of prior channel features, the K-SVD dictionary elements preserve the prior channel connectivity. This connectivity-preserving property of the K-SVD dictionary presents an important advantage over the SVD basis for representing complex (nonGaussian) and uncertain geologic features. Preserving the prior continuity, while accommodating significant structural uncertainty, is critical in the subsurface model calibration context. Similar results and conclusions were obtained for heterogeneous models (not shown). A particularly interesting case is when even the type of prior continuity model (or formation type) is not known a priori. We postpone the discussion for this case to a later section of the paper. Next we discuss the model calibration performance of the two methods.

125

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

Fig. 2. Twenty sample log-permeability maps (first column) for fluvial channels with the corresponding K-SVD dictionaries with K = 2025 (second column) and K = 400 (third columns) elements and the leading SVD basis elements (fourth column) for channels with (a) vertical or correct, (b) horizontal or wrong, and (c) random or uncertain directionality. The reference log-permeability is vertical and shown in Fig. 3 (top).

2.2. Model calibration 2.2.1. Example 1: non-Gaussian fluvial channels We evaluate the calibration performance of the K-SVD dictionary using two sets of waterflooding experiments in two-phase (oil/water) two-dimensional oil reservoirs. In these synthetic examples a 450  450  10 m reservoir (hereafter, Reservoir A) is discretized into a 45  45  1 uniform grid system with 10  10  10 m uniform block sizes (see Table 1 for general reservoir description). The configuration in Reservoir A and the true log-permeability for this example are shown in Fig. 3a. We use the iteratively reweighted least squares solution approach (discussed in Part I) under different K-SVD dictionary assumptions. The reconstruction results with the correct, incorrect and mixed (random) K-SVD dictionary are shown in Fig. 3b–d, respectively. The significant K-SVD dictionary elements for different prior model

Table 1 General simulation information for Reservoirs A and B. Parameter

Reservoir A

Reservoir B

Simulation parameters Phases Simulation time Grid systems Cell dimensions Rock porosity Initial oil saturation Injection volume Number of injectors Number of producers

Two-phase (o/w) 1 year 45 by 45 by 1 10 by 10 by 10 0.20 1.00 1PV 1 8

Two-phase (o/w) 1 year 60 by 220 by 1 10 by 10 by 10 0.20 1.00 1PV 13 2

Assimilation information Observation intervals Obs. at injection wells Obs. at production wells

12 days Pressure Pressure & saturation

15 days Pressure Pressure & saturation

126

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

Fig. 3. Log-permeability reconstruction results for Experiment 1: (a) Reference log-permeability with a 9-spot well configuration; (b) reconstructed log-permeability using the K-SVD dictionary and the correct prior library; (c) reconstructed log-permeability map using the K-SVD dictionary and the incorrect prior library; (d) reconstructed logpermeability model using the K-SVD dictionary and the random prior library; (e) reconstructed log-permeability using 20 leading SVD basis elements; (f1)–(f3) reconstructed log-permeability with 30, 40, and 50 leading SVD basis elements. In all cases the K-SVD dictionary was generated with K = 400 and S = 20. In (b)–(d), the plots on the right show the selected dictionary components while the plot to the right in (e) shows the leading S = 20 SVD basis elements.

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

127

Fig. 4. Log-permeability reconstruction results for Experiment 2: (a) Reference log-permeability with a 9-spot well configuration; (b) reconstructed log-permeability using the K-SVD dictionary and the correct prior library; (c) reconstructed log-permeability model using the K-SVD dictionary and the random prior library; (d) reconstructed logpermeability using 20 leading SVD basis elements; (e1)–(e3) reconstructed log-permeability with 30, 40, and 50 leading SVD basis elements. In all cases the K-SVD dictionary was generated with K = 400 and S = 20. In (b)–(c), the plots on the right show the selected dictionary components while the plot to the right in (d) shows the leading S = 20 SVD basis elements in.

128

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

assumptions are shown on the right columns of Fig. 3b–d. In all cases, we have assumed the dictionary size and sparsity level to be K = 400 and S = 20, respectively. From these results, it can be observed that: (i) under a correct prior library, the solution can reproduce the main channel structure in the true model; (ii) when an incorrect structural prior model is used the reconstruction results degrade significantly; (iii) the mixed dictionary provides a solution similar to that obtained with the correct dictionary, which is largely attributed to the selective nature of the sparsity regularization. As shown in Fig. 3d, the final solution is obtained by selecting and combing relevant components from the highly diverse dictionary. Interestingly, the inconsistent structures that are included in the mixed library do not corrupt the quality of final solution because they do not contribute to reconstruction of the solution. Fig. 3e shows the reconstruction results using the S = 20 SVD leading basis functions. The solution in this case is not only far from the reference model, it also does not reproduce any channel structures. The latter is expected since the SVD basis elements are strongly affected by the presence of inconsistent elements (see Fig. 2). Shown in Fig. 3f are the results of inversion with increasing numbers of SVD leading modes (S = 30, 40, 50). Increasing the number of basis elements in this case did not lead to any noticeable improvement. 2.2.2. Example 2: heterogeneous Gaussian model In Example 2, we use a heterogeneous Gaussian permeability distribution as the reference model (Fig. 4a). This reference model is generated using the sgsim algorithm and the variogram parameters described in Table 2 (Reservoir B). For this example, we present the solutions with: (i) a prior library that was generated from the same variogram used to obtain the reference model (correct prior), (ii) a highly uncertain prior library generated by randomizing the variogram ranges and direction of anisotropy (see the third column in Table 2 for details). The 3600 realizations that were generated for this example were not conditioned on hard data for increased variability. Fig. 4b shows the reconstruction results with the correct variogram model. The significant dictionary elements are shown to the right while the final reconstructed permeability is displayed on the left. As expected, the solution captures the main features in the reference model. Fig. 4c shows the reconstructed log-permeability (left) and selected dictionary elements for the case with random variogram models. Also in this case, the estimated log-permeability is very similar to the reference model. Fig. 4d shows the reconstruction results for S = 20 leading elements of the SVD basis with random variogram parameters. While the high and low trends in the solution are approximately identified, the solution is not nearly as good as that obtained with the KSVD dictionary. Lastly, Fig. 4e shows the reconstruction results with increased SVD leading elements, which also in this case do not seem to improve the solution. The data mismatch errors for each case were consistent with the results shown for log-permeability estimation (not shown for brevity). We next consider the effect of measurement noise.

Table 2 Correct and uncertain variogram parameters. Parameter name heterogeneous model

Reservoir A: correct case variogram parameters

Reservoir A: mixed case variogram parameters

Amax AR = Amin/Amax Amin Azimuth Dip Rake

60 0.3 1

[10, 200] [0.2, 1] 1 ½ p2 ; p2  0 0

3. Measurement uncertainty Another important consideration for an inversion method to be viable is its robustness to measurement noise. Inversion of noisy measurements must be implemented with care to account for the presence of noise and avoid noise-fitting. In particular, under over-parameterized models where model resolution is not reconciled with measurement resolution, the solution can become quite sensitive to data noise. We explore this and other possible effects of measurement noise on the behavior of the proposed sparsityregularized inversion algorithm in this section. 3.1. Effect of measurement noise Under a noise covariance matrix Cy and a dimensionless data misfit (misfit normalized by magnitude of observed values), the IRLS multiplicative objective function can be expressed as (see [29] and Part I [33])

J ¼ ðy  gðv ÞÞT C1 y ðy  gðv ÞÞSPðv ; e1 ; e2 Þ;

. 2 where the sparsity term is SPðv nþ1 ; e1 ; e2 Þ ¼ i v nþ1 þ e1 i  ð12pÞ ðv ni Þ2 þ e2 . Setting the first order derivative of this objective function equal to zero leads to the following minimization iterations

vnþ1 ¼

 1 h i n JTy C1 JTy C1 ðy  gðv ÞÞT C1 y Jy þ bW y y ðy  gðv ÞÞ þ Jy v ; ð2Þ @gðv Þ . @v

where Jy is the Jacobian matrix, i.e.  12 with weights Wki;i ¼ ðv ðk1Þ Þ2i þ ek

4

0 0

As in Part I, we use the

e1 ¼ e2 ¼ e ¼ ky  gðv Þk22 . In

our reconstruction example, we assume that the measured data are of similar quality and use C1 y ¼ Iy , where Iy is the identity matrix with the same dimension as the measurements. Now, we study the effect of observation noise in an experiment similar to Example 2. For each measurement, we assume additive zero-mean Gaussian observation noise with a standard deviation equivalent to 5% of its variation range in time. The noise is added to the data that is obtained by running a forward flow simulation with the synthetic reference log-permeability model. We also considered the effect of noise under different dictionary sizes. The reconstruction results corresponding to geologic dictionaries with the correct variogram model and K = 2025, 400, 200, 100 are shown in Figs. 5a and b. Fig. 5a shows the solutions without any data noise while Fig. 5b displays the results when 5% noise is added to the observations. In general, presence of noise appears to have resulted in underestimation of the log-permeability maps. In addition, as expected, the heterogeneity of the solution is reduced with decreasing dictionary size. Similar results are also presented in Figs. 5c and d with the random prior variogram model parameters. Fig. 5c shows the estimation results in the absence of data noise while Fig. 5d displays the solutions for data with 5% noise. Underestimating the log-permeability features is also evident in this case. When observation noise is present the mathematical form of the sparsity-regularized inversion can be expressed as:

min kv k1

p

ð1Þ P 

s:t: ky  gðv Þk2Cy 6 r2 ;

ð3Þ Þk2Cy

where kv k1 is the l1-norm of the coefficient vector, ky  gðv is the weighted data mismatch norm, and r2 is the energy of the observation noise. In this formulation, any reduction in data misfit norm below the noise variance is not expected to improve the solution. Nonetheless, minimization of kv k1 may proceed to further

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

129

Fig. 5. Reconstructed log-permeabilities with reduced size K-SVD dictionaries and 5% observation noise: (a1)–(a4) reconstructed log-permeabilities using the correct K-SVD dictionary with K = 2025, 400, 200 and 100 and assuming noise-free observations; (b1)–(b4) reconstructed log-permeabilities with the correct K-SVD dictionary for K = 2025, 400, 200 and 100 and assuming 5% additive Gaussian measurements noise; (c1)–(c4) reconstructed log-permeabilities with the K-SVD dictionary derived from the randomized prior with K = 2025, 400, 200 and 100 and assuming noise-free observations; (d1)–(d4) reconstructed log-permeabilities with the K-SVD dictionary derived from the randomized prior with K = 2025, 400, 200 and 100 and assuming 5% additive Gaussian measurement noise.

reduce the misfit term ky  gðv Þk2Cy . This effect can be mitigated by replacing kv k1 with the original definition of sparsity, i.e. kv k0 . However, as discussed in Part I, kv k0 is not well-behaved and introduces optimization complications. The underestimation effect is partly attributed to the sensitivity of the l1-norm to coefficient magnitudes (unlike the l0-norm that is only affected by the1number of non-zero coefficients). We note that P p p the lp ðv Þ ¼ norm can be minimized through two mechai jv i j nisms: (i) reducing the number of nonzero elements and (ii) decreasing the magnitude of the active significant coefficients. The latter is not an intended outcome of the sparse reconstruction and can result in solution underestimation. The underestimation effect is also expected to be sensitive to data noise. Under noisy data, among the solutions that reproduce the observed data within the noise variance, underestimated solutions are selected because of their lower l1 -norm cost function.

A two-step implementation procedure was introduced in [22] where in the first step they identified the sparsity structure (i.e., non-zero components of the vector v) through l1-norm regularization; and in the second step, only updated the magnitude of the identified non-zero components from the first step to minimize the data misfit term. We note that coefficient underestimation effect is diminished as p approaches zero. Therefore, to mitigate the underestimation effect, we consider minimizing the lp-norm (0 < p 6 1) where p is adaptively decreased during the iterations. That is, instead of l1-norm, we include lp-norm in the objective function and decrease p from an initial value of p = 1 to p = 0+ at the final iterations. Effectively, as p approaches zero the objective function becomes insensitive to the magnitude of the coefficients v as desired. Fig. 6a shows the reconstruction results for the noisy data by minimizing the l1-norm regularized cost function. The solution

130

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

Fig. 6. Reconstruction results with l1 and lp norm penalty functions: (a1)–(a3) reconstructed log-permeability, sparse solution coefficient vector, and sorted significant dictionary elements using l1-norm minimization with IRLS algorithm; (b1)–(b3) reconstructed log-permeability, updated sparse solution coefficient vector, and sorted dictionary elements using a simple least-squares algorithm where only the active elements in (a2) were updated; (c1)–(c3) reconstructed log-permeability, sparse solution coefficient vector, and sorted significant dictionary elements using lp-norm minimization with decreasing p as a function of iterations of the IRLS algorithm (p is decreased with iterations following the function shown in (c2)). In all cases, K = 2025 and S = 20 were used.

for v at the final iteration is shown in the middle column. The sparsity of the solution is clearly observed in this plot where only very few active components are present. In Fig. 6b, using the two-step algorithm in [22], we updated only the active components from the solution in Fig. 6a (middle column) to minimize the data misfit norm. This was aimed at removing possible solution underestimation effects and improving the data match. As can be verified in Fig. 6b, the under-estimation effect is clearly improved by assigning larger values to significant components. Fig. 6c shows the results for our proposed approach where the lp-norm with decreasing p value (with iterations) is used in the objective function. For this example, we have used the following function for reducing the p value:

1 u ½p ðwðiÞ þ 1Þ þ pl ð1  wðiÞÞ with 2    1 Nit  i wðiÞ ¼ erf p 2 1 ; erf ðpÞ Nit  1

pðiÞ ¼

ð4Þ

where i denotes the current iteration, Nit the maximum number of iterations, and pl and pu, the lower and upper limits of p. The Rx 2 notation erf represents the error function erf ðxÞ ¼ p2ffiffipffi 0 et dt. We note that w(i) has the range 1 6 wðiÞ 6 1 where the lower and upper limits correspond to pl and pu. In this paper, we set Nit = 30, pl = p0+ = 0.05 and pu = 1. The middle panel in Fig. 6c shows the behavior of p(i) as a function of iteration for the shown example. The left column of the figure shows the estimated log-permeability while the right column displays twenty of the active elements in the reconstruction. In comparison with Fig. 6a (left), the underestimation effect is significantly improved in Fig. 6c (left). Note that since the minimization objective functions are different in each case, the solution in Fig. 6c is different from those in Figs. 6a and b. The middle panel in Fig. 6c also shows that the final solution is very sparse and many of the coefficients are indeed very close to zero (due to near l0-norm regularization at later iterations). Fig. 7a plots the behavior of kv k0 (left) and kv k1 (right) functions during the minimization iterations (for monitoring purpose

131

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

Fig. 7. The behavior of l0 and l1 norm functions at different iterations of IRLS minimization when (a) l1-norm and (b) l0-norm are minimized.

simulation). Fig. 8a shows the well configuration for this example. The production schedule is summarized in Table 1 (third column). Figs. 8b and c show the reference log-permeability model and its approximation using the best (in RMSE sense) S = 20 dictionary elements, respectively. Note that in generating the approximation in Fig. 8c, it is assumed that the log-permeability is known and, hence, the best achievable solution is obtained. Fig. 9a1 shows the log-permeability reconstruction results using l1 -norm minimization and noise-free dynamic flow measurements collected every 15 days for a total of one year. The solution in this case identifies the main large-scale trends in the reference log-permeability field. Note that the log-permeability solution tends to capture the permeability trend at the vicinity of the observation points (well locations). The results at locations away from the wells are not as accurate. The reconstruction results after adding 2% noise to the data is shown in Fig. 9a2. For this example, because the range of variation in each data was quite large, a 2% noise level translated into a significant noise level. The reconstructed model in this case is slightly inferior, but still identify the major high and low permeability trends. As expected, addition of noise resulted in slightly underestimated solution. Fig. 9a3 shows similar results for noisy data, but when lp-norm minimization (with variable p) is used. The solution is very similar to that in Fig. 9a2, except that the underestimation effect is ameliorated. Fig. 9b shows the selected dictionary elements for each case in Fig. 9a. A closer examination shows that the selected elements for the noisy data (Fig. 9b2 and b3) are almost the same, but reordered. Some of the elements selected for reconstructing the solutions in all the three cases are the same. These elements contain important features that contribute significantly to identification of low and high trends in the reference model.

only) when l1-norm regularization is used with 5% data noise. It can be observed from Fig. 7a that while kv k1 increases during the initial iterations, kv k0 continuously decreases. The monotonic decrease in kv k0 with simultaneous increase in kv k1 during the early iterations implies that while the number of active (non-zero) coefficients is decreasing, the magnitude of remaining active coefficients is increasing. At later iterations, however, kv k1 begins to decrease while kv k0 remains constant. Since, at this point in the iterations, the significant coefficients are stabilized, this behavior suggests that the number of active coefficients remain unchanged while the magnitudes of these coefficients decrease, which explains the underestimation effect. Similar results are plotted for thekv kp function (with decreasing p) in Fig. 7b. Also in this case kv k0 norm is monitored to understand the behavior of the algorithm. The kv k0 function monotonically decreases with the iterations, similar to the behavior observed in Fig. 7a. On the other hand, the kv k1 function initially increases and then starts to decrease between iterations 8 and 15. However, since p decreases with the number of iterations, as p becomes smaller kv k1 begins to increase without an accompanying increase in kv k0 . Hence, it can be concluded that the increase in the kv k1 function is primarily due to the increase in the magnitude of the significant solution coefficients.

3.1.1. Example 3: SPE10 model We now apply the proposed approach to the top layer of the SPE10 model [39], which represents a Brent sequence. The top layer of this model has 60  220 = 13,200 gridblocks. The prior model realizations in this case are constructed using the sgsim algorithm with highly uncertain variogram model parameters to account for the uncertainty in anisotropy direction (see Table 3). We have used 3600 realizations to construct K = 500 K-SVD dictionary elements with a sparsity level of S = 50 (10%). In our example, we consider 13 water injectors and 2 oil producers (two-phase

3.1.2. Example 4: uncertainty in the formation type In some cases there may be uncertainty not only in the parameters of the continuity model, but also in the formation type. In this

Table 3 Uncertainty introduced for variogram model for SPE10 top layer. Parameter name

Amax

AR

Amin

Reservoir C: mixed case variogram parameters

[10, 200]

[0.2, 1]

1

Azimuth p p

2;2

Dip

Rake

0

0

132

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

(a) Field setup Inj 1

Inj 2

Inj 4

Producer 2

Inj 6

Inj 7

(c) K-SVD approx.

(b) True Lperm 5

5

55

55

110

110

165

165

Inj 3

Inj 5

Inj 8

Producer 1 Inj 10

Inj 9

Inj 11

Inj 12

Inj 13

220

220 5

5

60

30

60

30

Fig. 8. Well configurations in Experiments 3 with the top later of SPE10 model (a), the reference log-permeability taken from the top layer of SPE10 model (b), and approximated/compressed representation of the reference log-permeability model with the K-SVD dictionary assuming complete knowledge about the reference model using the OMP algorithm (c). The K-SVD dictionary is generated from 3600 realizations of the log-permeability that are simulated with randomized variogram parameters reported in Table 3 (dictionary size K = 500 and sparsity S = 50 are used).

(a1) Estimated Lperm

(a2) Estimated Lperm

(a3) Estimated Lperm

5

5

5

55

55

55

110

110

110

165

165

165

220

220

220 5

60

30

5

(b1) Active elements

5

60

30

(b2) Active elements

(b3) Active elements

0

0

0

220

220

220

440

440 0

60

120

180

60

30

240

300

440 0

60

120

180

240

300

0

60

120

180

240

300

Fig. 9. Log-permeability reconstruction and leading structures in Experiment 3; reconstructed log-permeabilities using the K-SVD dictionary and (a1) l1-norm minimization with no observation noise, (a2) l1-norm minimization with 2% observation noise, and (a3) lp-norm minimization with decreasing p values with iterations and assuming 2% observation noise; (b1)–(b3) show the ordered significant dictionary elements. In all cases K = 500 and S = 50 was used and the prior library and dictionary were identical to those in Fig. 8.

section, using both Gaussian and fluvial channel prior models, we investigate the performance of the proposed method under

uncertainty about the type of continuity model. We consider two approaches: first, constructing a geologically trained dictionary

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

from a library that has both channel and Gaussian formations; and second, generating one geologic dictionary for each formation type separately and solving the inverse problem with a combined dictionary. We call the latter a ‘‘hybrid dictionary’’ to distinguish it from a single dictionary. For this example, we use a setting similar to Examples 1 and 2. 3.1.3. Dictionary for the mixed library For this section we mix the library elements in Examples 1 and 2 to generate a highly uncertain geologic library. We apply the K-SVD algorithm to compress the prior library and construct a S = 20 sparse dictionary with K = 400 elements. The library elements and the related K-SVD dictionary and SVD basis are shown in Figs. 10a–c, respectively. The K-SVD dictionary elements tend to have both channel shape and heterogeneous features to be

133

able to approximate every model in the library. As expected, the SVD operation tends to combine the library elements and provide dominant features that do not resemble any of the formations types. The data integration results for the K-SVD dictionary and SVD basis are shown in Figs. 11a and b, respectively. The results are presented for both the heterogeneous and channel reference models. The superiority of the reconstruction results with the K-SVD dictionary is evident in both of these examples. 3.1.4. Hybrid dictionary Mixing prior library models without any consideration to their fundamental characteristic differences may not be the most effective way of constructing a geologic dictionary. An alternative approach is to generate two separate K-SVD dictionaries, one for

Fig. 10. K-SVD dictionary and SVD basis derived from a mixed prior library representing two different formation types: (a) prior library realizations; (b) corresponding 20 KSVD dictionary elements with K = 400 and S = 20; (c) S = 20 leading SVD elements.

Fig. 11. Reconstructed log-permeabilities and the corresponding significant dictionary elements: (a1) and (a2) reconstructed log-permeability and 20 selected K-SVD dictionary elements for the heterogeneous log-permeability model in Experiment 2; (b1) and (b2) reconstructed log-permeability and 20 selected K-SVD dictionary elements for the channel example in Experiment 1; (c1) and (c2) reconstructed log-permeability and 20 (out of 40) leading SVD basis elements for the Gaussian heterogeneous permeability model in Experiment 2; (d1) and (d2) reconstructed log-permeability and 20 (out of 40) leading SVD basis elements for the channel permeability model in Experiment 1.

134

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

channel example and one for Gaussian example and combine them into one single dictionary to perform the inversion. This approach would ensure that the geologic features are not blended, which is preferable especially when the two formations are completely distinct. In our example, we generated K=200 K-SVD dictionaries for each prior library in Examples 1 and 2 (see Fig. 12) and combined them to form a K = 400 dimensional dictionary. In this case, the solution can be expanded as

u20251 ¼ ½uG 20251 þ ½uCh 20251 ¼ ½UG 2025200  ½v G 2001 þ ½UCh 2025200  ½v Ch 2001 ;

ð5Þ

where uG is the contribution of the Gaussian dictionary and uCh is the contribution of the channel dictionary. The matrices UG and UCh respectively represent Gaussian and channel elements of the dictionary while vG and vCh denote the coefficient vectors corresponding to the Gaussian and channel dictionaries, respectively.

Fig. 13a and b show the reconstruction solutions with the resulting hybrid dictionaries. Shown in Fig. 13 are the selected elements from the combined dictionary and the individual contribution of the Gaussian and channel geologic dictionaries, calculated as ½uG 20251 ¼ ½UG 2025200  ½v G 2001 and ½uCh 20251 ¼ ½UCh 2025200  ½v Ch 2001 . Note that even though the formations types are distinct it is possible to have elements from one dictionary to contribute to reconstruction of a model from a different formation type, as can be observed in Fig. 13. One may also apply the reconstruction with each individual dictionary, i.e. by performing two separate model calibrations, and decide which prior library provides a better solution. This becomes rather cumbersome with increasing number of prior models. The discussion of this section can be easily extended to multiple formation types and multiple prior libraries. In that case, the hybrid dictionary combined with the selective nature of the sparsity-regularization provides an attractive calibration framework, an approach that is currently being investigated.

Fig. 12. Hybrid geologic dictionary obtained by combining the dictionaries that are separately obtained from the channel and Gaussian prior libraries; (a) K-SVD dictionary derived from the channel prior library with K = 400 and S = 20; (b) K-SVD dictionary derived form the Gaussian prior library with K = 400 and S = 20.

Fig. 13. Log permeability reconstructions with the hybrid dictionary: (a1) and (b1) 20 selected leading hybrid dictionary elements to reconstruct the Gaussian and channel log-permeability models, respectively; (a2) and (b2) overall contribution of the Gaussian log-permeability dictionary to the reconstructed Gaussian and channel logpermeability models, respectively; (a3) and (b3) overall contribution of the channel log-permeability dictionary to the reconstructed Gaussian and channel log-permeability models, respectively; (a4) and (b4) reconstructed Gaussian and channel log-permeability models, respectively.

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

4. Conclusions In Part II of the paper, we discussed the application of geologically learned sparse dictionaries (introduced in Part I) for subsurface flow model calibration under significant prior uncertainty and noisy measurements. We showed that the reconstruction results are far less sensitive to prior uncertainty than alternative prior-based parameterization techniques such as truncated SVD method. In fact, the diversity of the prior library can be a useful property because it can further enhance the sparsity of the solution since, in that case, many of the existing dictionary elements are likely to have little or no contributions to the solution. Unlike the SVD basis, the K-SVD dictionary is able to preserve the original continuity in the prior models despite significant diversity in the prior information. The selective property of the sparsity regularization exploits the redundancy in the dictionary rather elegantly by choosing a small number of relevant elements from a diverse set to reconstruct the solution. We also examined the performance of the algorithm under noise-corrupted measurements and showed that minimization of an l1-norm sparsity regularized objective function can lead to underestimated solutions. The original l0-norm sparsity regularization, on the other hand, does not exhibit this issue. Since l0-norm is not a mathematically well-behaved function and causes singularity during optimization, we proposed the use of an adaptive lp-norm sparsity regularization where during minimization iterations p decreases from 1 to a very small but non-zero value (0.05 in this case). We showed that this solution approach mitigates the underestimation problems as it practically reduces to an l0-norm at final iterations and shows less sensitivity to coefficient magnitude. Important developments in the sparsity-regularized model calibration include an ensemble-based implementation to quantify the solution and prediction uncertainties as well as non-uniqueness issues, more efficient algorithms for generating physics-based learned geologic dictionaries, methodologies for reducing the computational cost of the model calibration process such as iteratively thresholding algorithm, and application of the method to large scale realistic problems. The importance of sparsity is recently recognized in signal processing and is currently being widely studied as a useful and frequently encountered property in modeling and analysis of many systems in various disciplines. The theoretical breakthrough in sparse reconstruction was achieved by the compressed sensing theory of [15]. The approach has since received increasing recognition and is rapidly percolating into several branches of science and engineering. In this paper, we demonstrated a sample application of it in effectively integrating uncertain prior models with complex structural connectivity into the model calibration process. Sparsity is a rich concept and prevalent in geosciences that has a potential to transform several aspects of subsurface modeling, including characterization, numerical simulation, and reduced-order modeling. References [1] Aharon M, Elad M, Bruckstein A. K-SVD: an algorithm for designing over complete dictionaries for sparse representation. IEEE Trans Signal Process 2006;54(11):4311–22. [2] Ajo-Franklin JB, Minsley BJ, Daley TM. Applying compactness constraints to differential seismic travel time tomography. Geophysics 2007;72:R67–75. [3] Bhark E, Jafarpour B, Datta-Gupta A. An adaptively scaled frequency domain parameterization for flow data integration. J Petrol Sci Eng 2011;75(3– 4):289–303. [4] Bhark E, Jafarpour B, Datta-Gupta A. A generalized grid-connectivity based parameterization for production data integration (in revision). Water Resour Res; 2010b. [5] Britanak PC, Yip P, Rao K. Discrete cosine transform: general properties, fast algorithms, and integer approximation. Academic Press; 2006. [6] Caers J, Zhang T. Multiple-point geostatistics: a quantitative vehicle for integrating geologic analogs into multiple reservoir models. In: Grammer

[7]

[8] [9]

[10] [11] [12]

[13] [14] [15] [16]

[17] [18] [19] [20] [21]

[22]

[23]

[24]

[25]

[26]

[27]

[28] [29] [30]

[31] [32]

[33]

[34]

[35] [36] [37] [38]

[39]

135

GM. editor. AAPG memoir integration of outcrop and modern analogs in reservoir modeling, vol. 80, 2004. p. 383–94. Caers J, Hoffman BT. The probability perturbation method: a new look at Bayesian inverse modeling. Math Geol 2006;38(1):81–100. doi:10.1007/ s11004-005-9005-. Carle SF, Fogg GE. Modeling spatial variability with one and multi-dimensional continuous Markov chains. Math Geol 1997;29(7):891–918. Carle SF, LaBolle EM, Weissman GS, van Brocklin D, Fogg GE. Conditional simulation of hydrofacies architecture: A transition probability/Markov chain approach, in hydrogeologic models of sedimentary aquifers. In: Fraser GS, Davis JM, editors. Concepts in hydrogeology (society for sedimentary geology). Springer; 1998. p. 147–70. Carrera J, Alcolea A, Medina A, Hidalgo J, Slooten LJ. Inverse problem in hydrogeology. Hydrol J 2005;13(1):206–22. Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM J Sci Comput 1998;20:33–61. de Marsily G, Delhomme J, Delay F, Buoro A. Fourty years of inverse problems in hydrogeology. CR Acad Sci SerII A Earth Planet Sci 1999;vol. 329(2): 73–87. de Marsily G, Delay F, Goncalves J, Renard PH, Teles V, Violette S. Dealing with spatial heterogeneity. Hydrogeol J 2005;13(1):161–83. Deutsch CV, Journel AG. GSLIB – geostatistical software library and user’s guide. New York, Oxford: Oxford University Press; 1998. Donoho DL. Compressed sensing. IEEE Trans Inform Theory 2006;52(4):1289–306. Fogg GE, Carle SF, Green C. A connected-network paradigm for the alluvial aquifer system. In: Zhang D. editor. Theory, modeling and field investigation in hydrogeology: a special volume in honor of Shlomo P. Neuman’s 60th birthday, Boulder, Colorado, Geological Society of America Special Paper, vol. 348; 2000. p. 25–42. Gavalas GR, Shah PC, Seinfeld JH. Reservoir history matching by Bayesian estimation. Soc Petrol Eng J 1976;16(6):337–50. Hill MC, Tiedeman CR. Effective groundwater model calibration. New York, NY: John Wiley & Sons; 2007. Jacquard P, Jain C. Permeability distribution from field pressure data. Soc Petrol Eng J 1965:281–94. Jafarpour B. Wavelet reconstruction of geologic facies from nonlinear dynamic flow measurements. IEEE Trans Geosci Remote Sensing 2010;49(5):1520–35. Jafarpour B, Goyal VK, McLaughlin DB, Freeman WT. Transform-domain sparsity regularization for inverse problems in geosciences. Geophysics 2009;74(5):R69–83. Jafarpour B, Goyal VK, McLaughlin DB, Freeman WT. Compressed history matching: exploiting transform-domain sparsity for regularization of nonlinear dynamic data integration problems. Math Geosci 2010;vol.42(no.1): 1–27. Jafarpour B, Khodabakhshi M. A probability conditioning method (PCM) for nonlinear flow data integration into multipoint statistical facies simulation. Math Geosci 2010;vol. 43(2):133–64. Jafarpour B, McLaughlin DB. Reservoir characterization with discrete cosine transform. Part-1: Parameterization. Part-2: History matching. Soc Petrol Eng J 2009;14(1):182–201. Jafarpour B, Tarrahi M. Assessing the performance of the ensemble Kalman filter for subsurface flow data integration under variogram uncertainty. Water Resour Res 2011;47:W05537. doi:10.1029/2010WR009090. Khodabakhsi M, Jafarpour B. Flow-calibrated simulation of complex geologic facies from multiple uncertain training images. Water Resour Res 2011 [in revision]. Koltermann CE, Gorelick SM. Heterogeneity in sedimentary deposits: a review of structure-imitating, process-imitating, and descriptive approaches. Water Resour Res 1996;32(9):2617–58. 10.1029/96WR00025. Last B, Kubik K. Compact Gravity Invers: Geophys 1983;48:713–21. Li L, Jafarpour B. Effective solution of nonlinear subsurface flow inverse problems in sparse bases. Inverse Problems 2010;26(10):105016. Lu PB, Horne RN. A multiresolution approach to reservoir parameter estimation using wavelet analysis. In: SPE annual technical conference and exhibition, Dallas, TX, 1-4 October; 2000. Mallat S. A wavelet tour of signal processing. 2nd ed. Academic Press; 1999. Mirowski P, Tetzlaff D, Davies R, McCormick D, Williams N, Signer C. Stationarity scores on training images for multipoint geostatistics. Comp Geos 2008;41(4):447–74. M-Khaninezhad MR, Jafarpour B. Sparse geologic dictionaries for subsurface flow model calibration: Part I. Inversion Formulation:. Adv Water Resour 2011 [in revision]. Oliver DS, Reynolds AC, Liu N. Inverse theory for petroleum reservoir characterization and history matching. Cambridge, UK: Cambridge University Press; 2008. Parker R. Geophysical inverse theory. Princeton, NJ: Princeton Univ. Press; 1994. Portniaguine O, Zhdanov M. Focusing Geophys Invers Images: Geophys 1999;64:874–87. Sahni I, Horne RN. Multiresolution wavelet analysis for improved reservoir description. SPEREE 2005;8(1):53–69. Sarma P, Durlofsky LJ, Aziz K. Kernel principal component analysis for efficient, differentiable parameterization of multipoint geostatistics. Math Geosci 2008;40:3–32. SPE Comparative Solution Project. .

136

M.M. Khaninezhad et al. / Advances in Water Resources 39 (2012) 122–136

[40] Strebelle S. Conditional simulation of complex geological structures using multiple-point statistics.. Math Geol 2002;34(1):1–22. [41] Western AW, Bloschl G, Grayson RB. Toward capturing hydrologically significant connectivity in spatial patterns. Water Resour Res 2001;37(1):83–97. [42] Yeh TCJ, Lee CH, Hsu KC, Illman WA, Barrash W, Cai X, et al. A view toward the future of subsurface characterization: CAT scanning groundwater basins. Water Resour Res 2008;44:W03301. doi:10.1029/2007WR00637. [43] Yeh TCJ, Lee CH, Hsu KC, Tan KC. Fusion of active and passive hydrologic and geophysical tomographic surveys: the future of subsurface characterization. In: Hyndman DW. et al., Subsurface hydrology: data integration for properties

and processes, geophys. monogr. ser., vol. 171, AGU, Washington, DC; 2007. p. 109. [44] Yeh WWG. Review of parameter estimation procedures in groundwater hydrology: the inverse problem. Water Resour Res 1986;22:95–108. [45] Yeh WW-G, Yoon YS. Aquifer parameter identification with optimum dimension in parameterization. Water Resour Res 1981;17(3):664–72. [46] Zinn B, Harvey CF. 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 2003;39(3):1051. doi:10.1029/2001WR00114.