Sparse Randomized Maximum Likelihood (SpRML) for subsurface flow model calibration and uncertainty quantification

Advances in Water Resources 69 (2014) 23–37

Contents lists available at ScienceDirect

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

Sparse Randomized Maximum Likelihood (SpRML) for subsurface flow model calibration and uncertainty quantification Mohammadreza M. Khaninezhad, Behnam Jafarpour ⇑ Viterbi School of Engineering, University of Southern California, United States

a r t i c l e

i n f o

Article history: Received 23 February 2013 Received in revised form 17 February 2014 Accepted 18 February 2014 Available online 1 March 2014 Keywords: Sparsity Compressed sensing Laplace priors Bayesian inversion Randomized Maximum Likelihood Model calibration

a b s t r a c t Despite their apparent high dimensionality, spatially distributed hydraulic properties of geologic formations can often be compactly (sparsely) described in a properly designed basis. Hence, the estimation of high-dimensional subsurface flow properties from dynamic performance and monitoring data can be formulated and solved as a sparse reconstruction inverse problem. Recent advances in statistical signal processing, formalized under the compressed sensing paradigm, provide important guidelines on formulating and solving sparse inverse problems, primarily for linear models and using a deterministic framework. Given the uncertainty in describing subsurface physical properties, even after integration of the dynamic data, it is important to develop a practical sparse Bayesian inversion approach to enable uncertainty quantification. In this paper, we use sparse geologic dictionaries to compactly represent uncertain subsurface flow properties and develop a practical sparse Bayesian method for effective data integration and uncertainty quantification. The multi-Gaussian assumption that is widely used in classical probabilistic inverse theory is not appropriate for representing sparse prior models. Following the results presented by the compressed sensing paradigm, the Laplace (or double exponential) probability distribution is found to be more suitable for representing sparse parameters. However, combining Laplace priors with the frequently used Gaussian likelihood functions leads to neither a Laplace nor a Gaussian posterior distribution, which complicates the analytical characterization of the posterior. Here, we first express the form of the Maximum A-Posteriori (MAP) estimate for Laplace priors and then use the Monte-Carlo-based Randomize Maximum Likelihood (RML) method to generate approximate samples from the posterior distribution. The proposed Sparse RML (SpRML) approximate sampling approach can be used to assess the uncertainty in the calibrated model with a relatively modest computational complexity. We demonstrate the suitability and effectiveness of the SpRML formulation using a series of numerical experiments of two-phase flow systems in both Gaussian and non-Gaussian property distributions in petroleum reservoirs and successfully apply the method to an adapted version of the PUNQ-S3 benchmark reservoir model. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Predictive modeling of fluid flow and displacement processes in geologic formations, such as groundwater aquifers and hydrocarbon reservoirs, hinges upon the accuracy and representativeness of the underlying numerical model(s) used to generate them. Construction of high-resolution geologic models from prior information and dynamic response data is of particular interest to development of subsurface resources. High resolution representations of subsurface flow models are desirable for detailed forward modeling and prediction of the system response to alternative ⇑ Corresponding author. Address: 925 Bloom Walk, HED 313, Mork Family Department of Chemical Engineering and Material Science, University of Southern California (USC), Los Angeles, CA 90089, United States. Tel.: +1 (213) 740 2228. E-mail address: [email protected] (B. Jafarpour). http://dx.doi.org/10.1016/j.advwatres.2014.02.005 0309-1708/Ó 2014 Elsevier Ltd. All rights reserved.

development strategies as a basis for optimizing resource recovery performance. However, building predictive flow models is challenging for several reasons, including the uncertainty in the prior knowledge about rock properties, the complexity of the connectivity in geologic formations, the nonlinearity between model parameters and observed quantities, and the lack of sufficient spatial data to constrain the variability in uncertain rock property distributions. In addition, finely discretized subsurface flow models that are typically used for high resolution numerical simulation often lead to an over-parameterized description of the inverse problem. Since grid-based high-resolution representations of distributed model inputs result in a large number of unknown parameters relative to available data, the resulting imaging inverse problem becomes ill-posed (e.g., [8–10,13,20,24,30,41,44,48]). Solving ill-posed inverse problems typically involves the minimization of a suitably defined data mismatch objective function

24

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

subject to specified constraints, known as regularization terms, that reflect direct or indirect prior knowledge about the solution. In the last decade, significant progress has been made in conditioning grid-based subsurface models to flow and transport data both in hydrogeology and petroleum engineering [13,23,41]. Deterministic and stochastic inversion algorithms with varying levels of complexity and different regularization techniques have been developed and applied to subsurface flow inverse problems to characterize the hydraulic properties of groundwater aquifers and petroleum reservoirs. An intrinsic attribute of geologic formations that is of fundamental importance in regularizing and/or parameterizing their physical properties is geologic continuity. The continuity in geologic formations translates into large-scale spatial correlations that makes reservoir properties amenable to low-order representations. Two methods that are commonly used to mitigate the ill-posedness of the subsurface flow inverse problem are (i) reducing the number of unknown parameters, i.e. parameterization (e.g., [20,24,30,44]), and (ii) incorporating prior information in the form of regularization constraints (e.g., [35,41,47]). The existing parameterization/regularization methods make restrictive a priori assumptions about the form of the solution, which introduces bias and lacks robustness and flexibility. Correlated rock property maps generally lend themselves to sparse (compact) representations in a properly designed decorrelating basis. That is, the connectivity features in grid-based descriptions of heterogeneous rock property distributions can be more compactly described using a featurebased representation. A feature-based description gains its efficiency from compacting the high-dimensional connectivity patterns into a set of known representative basis elements. The basis elements can be generic such as those used in image compression or they can be learned from a prior training data. The main characteristic of the resulting basis elements is that they are highly effective in reproducing the most salient connectivity features in the solution. This view transforms the original gridbased estimation of subsurface flow properties to a feature selection problem. Many realistic geologic systems have complex architectures with spatial connectivity patterns that are not amenable to conventional second-order geostatistical descriptions (e.g., [7,14,21,50]). Representing nonlinear, non-Gaussian geologic features, such as meandering channels that generate important preferential flow and transport patterns requires more sophisticated higher-order statistical simulation techniques such as multipoint geostatistics [5,22,45]. The subsurface flow inverse problem in such cases can be significantly facilitated by incorporating prior knowledge such as a training image or a prior geological dataset. The resulting prior training dataset can be used to learn compact geologic basis functions for parameterization. Khaninezhad et al. [31] present learned ‘‘sparse geologic dictionaries’’ for compact and flexible representation of geologic variability in a given formation. A sparse geologic dictionary in this context refers to a diverse collection of geologic features that are learned from a prior training data with the property that only a small subset of these dictionary elements is needed to represent the patterns in the prior information. Such a dictionary can be used to formulate and solve the subsurface flow inverse problem using sparse reconstruction techniques [31]. Data scarcity and imperfect assumptions in describing subsurface flow and transport processes often lead to significant uncertainty in characterization and prediction of the behavior of these system. In particular, proper representation of uncertainty in spatially distributed rock properties can greatly impact the subsurface flow and transport behavior, and hence resource development strategies. Underestimation of uncertainty in model calibration

can result in biased solutions and future predictions that are difficult to detect and correct. On the other hand, overestimating the uncertainty can lead to under-constrained solutions where many distinct, but plausible, scenarios can honor the existing data but provide very different forecasts, thereby rendering the inverse modeling exercise inconclusive. In general, however, overestimation of uncertainty is preferable to its underestimation since the latter can misguide future development and result in irreversible losses whereas the former often underscores the need for additional constraining data [29]. In general, the main sources of uncertainty in solving an inverse problem are: (i) the uncertainty in the knowledge about the unknown model parameters (i.e., prior uncertainty), (ii) the uncertainty in the model that is used to relate the input parameters to the observed data (a.k.a., modeling error), and (iii) the uncertainty in measuring and representing the observed outputs (i.e., observation errors). In the Bayesian formulation of an inverse problem, measurement errors are typically represented through a likelihood function. A topic less explored in the literature is the representation of uncertainty in the modeling assumption, or the modeling error, partly because proper characterization of the modeling error is not trivial. In many cases, the errors introduced in deriving the governing flow equations is considered to be smaller than the uncertainty in describing the subsurface physical properties as model input parameters. Hence, in many cases the inversion is presented as the integration of a prior model PDF that reflects the uncertainty in rock properties (model parameters) with a likelihood function that represents the uncertainty in measurements and their predictions. The cost and difficulties associated with data acquisition coupled with the complex heterogeneity in the distribution of subsurface properties place significant uncertainty on the prior information. Typically, dynamic flow data provide aggregate information about rock properties and limit the resolution at which model parameters can be estimated. The limitations in the amount and resolution of the dynamic flow data lead to significant level of unresolved uncertainty in the model and its predictions even after data integration. The uncertainty in describing reservoir properties can also affect the type and level of model parameterization. For example, if reliable prior information is available, one would prefer to strongly constrain the parameterization and, hence, the solution of the inverse problem by using the existing prior information. A popular example of these types of parameterization is the Principle Component Analysis (PCA) [25,34], which constrains the parameterization to honor the second-order statistical moment of the parameter field. The PCA provides an effective approach for parameter reduction, which is quite useful for model calibration [34]. However, placing too much emphasis on the prior information can introduce bias into the solution. On the other hand, when prior information is unreliable or unavailable, one would have to depend on the available data to determine the main spatial connectivity in the field. In general, even qualitative prior knowledge/assumptions based on the physics of the problem may be useful for constraining the inverse modeling solution. For instance, Tikhonov type regularization techniques may be used to impart global spatial smoothness on the parameters, and low-frequency Fourier-type subspaces may be suitable for parameterization of correlated and smooth subsurface properties. A major limitation of such regularization measures is that they are not universally applicable and can introduce bias by favoring overly smooth solutions. In general, even qualitative knowledge about the formation type and the form or shape of the expected geologic patterns can be quite helpful in prescribing an effective and suitable parameterization or regularization method. The Bayesian approach [46] provides an elegant probabilistic framework for conditioning uncertain prior model parameters on

25

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

available flow data. However, complete closed-form characterization of the posterior distributions can only be obtained for very simple forms of probability density functions (such as jointly Gaussian data and model parameter distributions). In practice, however, real parameter fields rarely follow such simple distributions and approximate solution methods must be sought. Two classes of approximate solutions are commonly used in practice. The first approach estimates a few point statistics of the posterior distribution (mainly the first and second order statistical moments) while the second class of methods use Monte-Carlo simulation techniques to approximately represent the posterior distribution via a finite set of conditional realizations. By generating multiple conditional model realizations, ensemble approaches provide a systematic Monte-Carlo simulation-based method for quantifying the flow prediction uncertainties (e.g., [1,42]). The ensemble Kalman filter (EnKF) [1,17], Markov-Chain Monte-Carlo (MCMC) [46], and the randomized maximum likelihood (RML) [33,43] are examples of ensemble methods that have been applied to subsurface inverse modeling problems. In this paper, we use sparse geologic dictionaries that are learned from prior training datasets to sparsely represent prior geologic models. We formulate a sparsity-promoting inversion approach using the Bayesian framework with the Laplace distribution to describe the sparse representation of the prior models. We solve the resulting Bayesian inverse problem using the Sparse RML (SpRML) approach as practical method for uncertainty quantification. The minimization involved in the SpRML approach is carried out using the iteratively reweighted least-squares (IRLS) algorithm [6,36]. To examine the performance of this approach, we apply this framework to two-dimensional and three-dimensional model calibration problems for uncertainty quantification under sparsity assumptions. 2. Bayesian formulation 2.1. RML with multi-Gaussian priors The general Bayesian formulation of the inverse problems that arise in multiphase flow in porous media is expressed as [41,46]

pðdobs juÞpðuÞ pðujdobs Þ ¼ pðdobs Þ

ð1Þ

where u 2 RN is the unknown rock property (e.g. log-permeability) and dobs 2 RM is the vector of measured data (e.g., well pressure and fluid flow rates). In this paper, we only consider the integration of dynamic data; however, we point out that in many cases direct measurements of the properties at the well location are available and can be included in model calibration. The probability density functions (PDF) pðuÞ, pðdobs Þ, pðdobs juÞ and pðujdobs Þ represent the prior, observation, likelihood, and conditional (posterior) distributions. In the case of jointly Gaussian model and data distributions, the posterior PDF can be characterized as pðujdobs Þ / expðJðuÞÞ, where T

T 1 JðuÞ ¼ ðgðuÞ  dobs Þ C1 d ðgðuÞ  dobs Þ þ ðu  up Þ Cu ðu  up Þ:

ð2Þ

where gðuÞ 2 RM is the vector of simulated measurements, Cd is the data covariance and up and Cu are the mean and covariance, respectively, of the prior distribution of u. The simulated data are usually nonlinearly related to model parameters u; hence, JðuÞ is a nonlinear function of the unknown parameters u. For linear models, i.e., gðuÞ ¼ Gu, the function JðuÞ takes a quadratic form, and the posterior PDF pðujdobs Þ becomes multi-Gaussian. In that case, the analytical form of the posterior PDF is readily to provide a full characterization of the solution [41,46].

The complete characterization of the posterior PDF in the case of a nonlinear model is not trivial. Additionally, even under linear models, if the multi-Gaussianity assumption about the prior model is not justified, analytical characterization of the posterior PDF may not be possible. Alternative methods for approximate characterization of the posterior PDF in each of the above cases have been used in the past. The MCMC method is a computationally demanding but rigorous approach for sampling from the posterior distribution [46]. The EnKF is a Monte-Carlo implementation of the Kalman filter that has been introduced as a computationally efficient and practical data integration approach for nonlinear models. The original Kalman filter is an optimal state estimator for linear statespace dynamical systems with jointly Gaussian statistics. For linear models with jointly Gaussian states and measurement PDFs, the EnKF can be shown to asymptotically (for large sample sizes) converge to the Kalman filter. Despite its theoretical limitations for handling nonlinear non-Gaussian problems, application of the EnKF to several large-scale nonlinear data assimilation problems has produced promising results [1,17]. The application of the EnKF to subsurface flow model calibration has been extensively studied in the past decade and important observations have been reported [1,12,18,40,49]. The RML is another method that also provides an approximate sampling approach for nonlinear and non-Gaussian problems [33,43]. For nonlinear models gðuÞ with jointly Gaussian distribution, nens samples from the posterior distribution in the RML approach are obtained by solving the following minimization problem [33,41,43] T

i

i

i min Jðui Þ ¼ ðgðui Þ  dobs Þ C1 d ðgðu Þ  dobs Þ ui

T

i i i ¼ 1; 2; . . . ; nens þ ðui  uiuc Þ C1 u ðu  uuc Þ

ð3Þ

where uiuc is an unconditional realization from the prior distribution i of u with uiuc  Nðup ; Cu Þ, and dobs is a random realization of the i observation vector taken from dobs  Nðdobs ; Cd Þ [33,43]. In this context, Nðup ; Cu Þ represents a multivariate Normal Gaussian distribution with mean up and covariance matrix Cu . For linear models, i.e., gðuÞ ¼ Gu, and jointly Gaussian PDFs, the closed-form solution of the above minimization problem can be expressed as [41]

 1   i ui ¼ uiuc þ Cu GT G Cu GT þ Cd Guiuc  dobs :

ð4Þ

However, in the linear case when the full conditional multiGaussian PDF is available, one can directly sample from it using other direct sampling methods. The RML minimization can be interpreted as an approximate conditional sampling in which unconditional realizations are slightly perturbed to reproduce the observed data within the measurement error given by Cd . Linearized error analysis can be used to show that the conditional samples generated with the RML approach reproduce the conditional mean and covariance matrix as long as the distance between the conditional mean and the simulated conditional samples is small enough so that the same linearization applies to both of them [33,43]. When this condition does not apply, the generated samples can only be considered as approximations and are not necessarily samples from the theoretical conditional distribution. In general, one could use an importance sampling criterion such as the Metropolis–Hasting approach to accept or reject the simulated samples. Despite the approximations involved, application of the RML method to subsurface flow model calibration problems in the past has been successful [19,41,43]. For nonlinear models, the above minimization problem must be solved as many times as the number of conditional realizations needed. However, since minimizing the above objective function involves several timeconsuming forward flow simulation runs, in practice the number of samples generated with this method is limited. The RML

26

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

formulation presented above is general and can be used with other types of prior distributions. In the next section, we consider the RML formulation under sparse priors. 2.2. Sparse RML with Laplace priors Given a set of expansion functions fui gi¼1;2;...;n , we express a property of interest u (e.g., gridblock log-permeability values) as a linear expansion of the form

u¼

n X

v i ui

ð5Þ

i¼1

where v ¼ ½v 1 ; v 2 ; v 3 ; . . . ; v n T are the expansion coefficients that represent u in fui gi¼1;2;...;n ; the U matrix is the linear transformation in sparse domain named the dictionary in sparse signal reconstruction literature. The vector v is said to be s-sparse if only s out of its n elements are non-zero. In that case, Eq. (1) can be reduced to an sterm expansion

u¼

s X

v i ui

ð6Þ

i¼1

The sparse reconstruction problem involves the estimation of the parameter u from its incomplete observations dobs ¼ gðuÞ. Since u is linearly related to v (i.e. u ¼ Uv ) we can write nonlinear and linear measurement equations in terms of v , as dobs ¼ gðv Þ and dobs ¼ Gv, respectively. The sparse reconstruction literature focuses on recovering the sparse vector v from its incomplete observations dobs ¼ Gv by minimizing the support (number of non-zero coefficients) of v (i.e., kv k0 ) subject to measurement constraints, i.e.

ðP0 Þ :

min kv k0

s:t: Gv ¼ dobs

ð7Þ

Recent developments in statistical signal processing that have resulted in the compressed sensing paradigm, specify that under certain, relatively mild, conditions one can find the sparse solution to (P0) by solving the following (P1) problem [11,15,16]

ðP1 Þ :

min kv k1 s:t: Gv ¼ dobs

ð8Þ

Reformulation of the (P1) problem to allow for observation errors leads to [16]

ðPe1 Þ :

min Jðv Þ ¼ kGv  dobs k22 þ kkv k1

ð9Þ

where k is the regularization parameter. The norms for t > 0 used in 1 P t t the above equations are defined as kxkt ¼ . i¼1 ðxi Þ The problem ðP e1 Þ provides the inspiration for formulating the Bayesian form of the sparse reconstruction problem. This can be achieved by writing the Bayesian formulation for sparse parameters v using the Laplace (double exponential) distribution [37]. That is, we consider the prior PDF of v to take the form

pðv Þ / expðkjv  v p jÞ

ð10Þ

where v p is the prior mean of v and k is the scale parameter. This formulation assumes that each component of the sparse vector v is uncorrelated with other elements in this vector. In addition, it is assumed that the parameter of the Laplace distribution for each element is identical and constant. With a Gaussian choice for the likelihood function, the posterior PDF from the Bayes’ rule reduces to pðv jdobs Þ / expðJðv ÞÞ with T

Jðv Þ ¼ ðgðv Þ  dobs Þ C1 d ðgðv Þ  dobs Þ þ kjv  v p j1 :

ð11Þ

For traditional sparse priors we have v p ¼ 0 (the expected value of the coefficients is zero, as implied by sparsity) and the resulting exponential argument in the posterior PDF can be simplified to T

Jðv Þ ¼ ðgðv Þ  dobs Þ C1 d ðgðv Þ  dobs Þ þ kjv j1 :

ð12Þ

P The l1 norm in the above equations is defined as jv j1 ¼ nj¼1 jv jj . The posterior PDF in this case has a complex form and closed-form solutions cannot be derived. Therefore, we use the RML sampling approach to find realizations from the posterior PDF. This is achieved by solving the following minimization

n o T i i i i i min Jðv i Þ ¼ min ðgðv i Þ  dobs Þ C1 d ðgðv Þ  dobs Þ þ kjv  v uc j1 : vi

vi

ð13Þ to generate the realization i from the posterior distribution. While for nonlinear problems the parameter k is difficult to exactly specify, a small number of trial and error attempts can be used to identify a reasonable value for k. We follow the iteratively reweighted least-square (IRLS) method proposed in [36] to minimize the above objective function. We write the IRLS objective function as i

T

i

i i i Jðv i Þ ¼ ðgðv i Þ  dobs Þ C1 d ðgðv Þ  dobs Þ þ k SPðv  v uc Þ:

where

ð14Þ

at

iteration (k + 1) of the minimization, we have P i;j i;j 2  v ¼ nj¼1 W i;j SPðv k  ððv kþ1  v uc Þ þ eÞ with the weights de1   2 2 ¼ v i;jk  v i;juc þ e . Note that at iteration (k + 1) fined as W i;j k i kþ1

i uc Þ

the weight of the sparsity term is written in terms of the solution at the previous iteration (k). The parameter e is generally a small positive constant included to insure numerical stability of the solution since the derivative of the l1-norm is not defined at the origin. In this study, however, we follow the approach in [36] and set these parameters to kgðv

2 obs kC d kÞ

k Þd

SPðv

ekþ1 ¼ eðv k Þ ¼ kgðv k Þ  dobs k2C1 and kkþ1 ¼ kðv k Þ ¼ d

. The formulations presented in Eqs. (13) and (14) tend

to promote sparse updates, i.e. minimize jv i  v iuc j1 , to the unconditional realizations. However, since the unconditional realizations are themselves sparse, the solution is also expected to be sparse. For minimizing the above objective function, we use a gradientbased optimization method, known as the iteratively reweighted least-squares approach, as described in [31,36]. Finally, to implement the RML algorithm, we need to define the basis functions f/i gi¼1;2;:::;n . While generic compression transforms such as discrete cosine transform (DCT) [3,4,25–27] or discrete wavelet transform (DWT) [28,39] or the classical PCA or KLT parameterization can be used as sparse basis functions, in this study we use a recently introduced geologically-trained sparse K-SVD dictionary [2,31,32]. The K-SVD algorithm and its important properties for subsurface flow inverse problems have been discussed in our recent publications [31,32]. The K-SVD dictionary elements f/i gi¼1;2;...;n are derived from a prior training dataset UnL ¼ fui gi¼1;2;...;L and have the property that they can approximate any model similar to those in the prior library with a small number of dictionary elements, s  n, in the linear expansion. Construction of a sparse geologic dictionary is performed by providing L prior model realizations as columns of a the matrix UnL ¼ fui gi¼1;2;...;L , and specifying a dictionary size K (with K elements) and the sparsity level S. In practice, S and K are specified by the user based on the complexity of the model and the desired compression level, see [31] for additional details on specifying K and S. The algorithm will compress the input model to generate an S-sparse dictionary UnK and the corresponding sparse representations of the prior models VKL ¼ fv i gi¼1;2;...;L such that only S (different) elements are needed to provide an approximate representation for each prior model realization. The following optimization problem can be solved to generate the dictionary elements:

^ ¼ arg min kU  UVk2 ; ^ Vg fU; F U;V

s:t: kv i k0 6 s;

i ¼ 1; 2; . . . ; L ð15Þ

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

We note that if instead of the minimum support constraint kv i k0 6 s; i ¼ 1; 2; . . . ; L the orthonormality constraint is imposed, the Singular Value Decomposition (SVD) of the prior library would result. The above optimization problem can be efficiently solved using a two-step iterative K-SVD algorithm [2]. Each iteration of the KSVD algorithm consists of two steps: the first step (known as sparse coding) is used to find sparse representations V for the initial dictionary U and prior models U by solving the sparse reconstruction sub-problems of the form

v^ i ¼ argmin kui  Uv i k2F ; v i

s:t: kv i k0 6 s

ð16Þ

The second step updates the dictionary by fixing V and U and solving the sub-problem

U ¼ arg min kU  U

UVk2F

2    K X  T ¼ arg minU  /i v i  U   i¼1

ð17Þ

F

where /i and v Ti denote the ith-column of U and the ith row of V, respectively. The K-SVD algorithm is initialized by using K members of the initial prior library as the dictionary. This initial dictionary undergoes a series of dictionary update and sparse coding steps, using Eqs. (16) and (17), such that any element in the prior library can have a S-sparse approximation in the final converged dictionary. Since the K-SVD algorithm decouples the sparse dictionary reconstruction into two steps it belongs to the class of greedy algorithms that do not have a formal convergence proof. However, the K-SVD algorithm provides a practical and efficient approximate solution [2]. There are several intrinsic differences between the K-SVD and SVD algorithms that result in different behavior and performance for the two methods (see Khaninezhad et al. [31,32] for details). First, the SVD basis is constructed by decomposing the images into

27

their spectral content, whereas the K-SVD dictionary elements are initialized with the prior library and gradually updated to generate a sparse dictionary. As a result, the dictionary elements in the KSVD algorithm tend to inherit the spatial connectivity of the prior library, which can be quite important in constructing complex non-Gaussian random fields, especially in an ill-posed inverse problem where prior information plays an important role. Second, the SVD construction tends to concentrate a large fraction of the energy (variance) of the images into a small number of sorted leading basis images. On the other hand, the dictionary elements in KSVD do not follow a particular order and, on average, have similar contributions to reconstruction of the prior models. The ordering of the SVD basis implies that the first S leading vectors can be selected a priori to provide S-term approximations to the prior models. In general, non-leading SVD basis images tend to represent the details (high frequency information) in an image that have relatively negligible contribution to representation of natural images. The a priori selection (truncation) of the SVD basis elements is convenient and is known to provide the lowest approximation RMSE across the entire library; however, a priori truncation is not optimal for approximating each individual element in the prior library, even in RMSE sense. On the other hand, the K-SVD algorithm solves a sparse coding problem (minimization problem) by searching over the entire dictionary to find S individual elements that have significant contribution to S-term approximation of each individual element. As such, unlike the truncated SVD where the selected S terms are fixed regardless of which image in being approximated, the selected S terms in K-SVD change from one image to another. In approximation theory, the former approach is referred to as linear approximation to differentiate it from the latter method that is known as non-linear approximation. These differences between the two methods can explain the better performance of the K-SVD algorithm in solving ill-posed inverse problems. It is also noteworthy that this improved performance comes with some cost; mainly the K-SVD dictionary does not possess the orthogonality

Coef. Quantiles

Coef. Quantiles

(a) Q-Q plot for Gaussian distribution

Standard Gaussian Quantiles

Coef. Quantiles

Coef. Quantiles

Standard Gaussian Quantiles

Standard Gaussian Quantiles

Standard Gaussian Quantiles

Coef. Quantiles

Coef. Quantiles

(b) Q-Q plot for Laplace distribution

Laplace Quantiles Coef. Quantiles

Quantiles

Laplace Quantiles

Laplace Quantiles

Laplace Quantiles

Fig. 1. The Q–Q plot of the sparse K-SVD coefficients againsts the standard Gaussian distribution (a) and the Laplace distribution (b) for a one-dimensional example with discrete facies maps. The plots show that the Laplace distribution is more appropriate than the Gaussian distribution for describing the sparse K-SVD coefficients.

28

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

Well Setup

Reference Log-Perm

(a)

(s1)

(s2)

(s3)

( s 4)

M e an

Va r ia n ce

(b)

(c)

Fig. 2. (a) Field setup (left) and reference log-permeability (right); (b) four unconditional log-permeability realizations (Columns 1–4), unconditional sample mean (Column 5) and variance (Column 6); (c) conditional log-permeability realizations (Columns 1–4) corresponding to (b), and the conditional sample mean (Column 5) and variance (Column 6). The sample mean and variance are calculated from 50 realizations.

of the SVD basis, is computationally more demanding to construct, can lead to non-unique reconstruction results, and requires a sparse reconstruction formulation to solve the inverse problem. In this paper, we use the K-SVD as the sparsifying expansion and refer interested readers to [31,32] for additional details about the K-SVD algorithm and its properties. The prefix K in the K-SVD algorithm refers to the size of the dictionary and is used in the name to reflect the K times rank-one SVD operation that is implemented in the algorithm. To test the validity of the Laplace prior on the K-SVD coefficients, we performed a Q-Q plot test with onedimensional discrete (non-Gaussian) facies models. For this test, we generated 40,000 one-dimensional realizations with the SNESIM algorithm and specified a dictionary size of K = 100 and a sparsity level of S = 20. The K-SVD sparse coefficients were generated using the Orthogonal Matching Pursuit (OMP) algorithm with the same sparsity level. The Q–Q plots for four of the coefficients against the Normal and Laplace distributions are shown in Fig. 1a and b. The linear trend in Fig. 1b suggests that the Laplace prior is a better assumption than the multi-Gaussian prior, which is consistent with our formulation. In addition, investigation of the covariance matrix of the K-SVD coefficients revealed that the coefficients have little correlations with each other. 3. Numerical experiments We use numerical experiments from two-phase flow in hydrocarbon reservoirs to evaluate the performance of the SpRML algorithm. Application of the method to three numerical experiments involving two two-dimensional synthetic models and the threedimensional PUNQ-S3 reservoir model is considered. In each case, we start with nens = 50 prior log-permeability realizations and condition these samples using the RML algorithm. In this paper, we do not apply any importance sampling approach to accept or reject the solutions obtained from the RML algorithm. We note, however, that larger sample sizes are necessary for approximating the statistics of the resulting conditional distribution. The first example resembles a multi-Gaussian-type distribution for the variability in the spatial permeability map, which is amenable to classical variogram-based geostatistical description methods. The second example represents a meandering channel shape connectivity that is not amenable to variogram-based representations. These two

examples are used to compare the performance of the conventional model calibration methods with the sparsitypromoting algorithm using the K-SVD dictionary. Finally, the third example demonstrates the application of the method to a more realistic-looking PUNQ-S3 3D example. The immiscible two phase fluid flow equations (e.g., oil and water) in porous media can be derived from the mass conservation principle and Darcy’s law and can be compactly represented as:

r

r





   ko @ So þ qo kðrPo  co rZÞ ¼ / @t Bo Bo

ð18Þ

   kw @ Sw þ qw kðrPw  cw rZÞ ¼ / @t Bw Bw

ð19Þ

where k and / are the rock permeability tensor and porosity, respectively; B, k and c are formation volume factor and phase mobility and density, respectively. The pressure and saturation are denoted as P and S while the source/sink terms are denoted as q. The subscripts o and w stand for oil and water. Assuming a diagonal tensor for permeability for a two dimensional domain in Cartesian coordinates, the above equations reduce to

      @ ko @Po @ ko @Po @ So þ ¼ þ qo kxx kyy / @x Bo @y Bo @t @x @y Bo

ð20Þ

      @ kw @Pw @ kw @Pw @ Sw kxx kyy / þ ¼ þ qw @x Bw @y Bw @t @x @y Bw

ð21Þ

These two equations can be combined with two constitutive equations, namely the capillary pressure equation and the physical saturation constraint, that is, Table 1 Geostatistical simulation parameters for the two-dimensional reservoir in Example 1. Property name

Consistent library

Kriging type Maximum conditioning data Variogram type Ranges max, med, min Angles azimuth, dip, rake Mean value (LogPerm) Standard deviation (LogPerm)

Ordinary kriging 12 Spherical 50, 10, 5 30, 0, 0 4 1

29

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

(a) K-SVD Dictionary Elements

(b) Unconditional Samples

(c) Updates Applied

(s1)

(s2)

(s3)

(s4)

Fig. 3. (a) 25 Sample K-SVD dictionary elements (out of K = 2500) with sparsity level S = 20; (b) sparse coefficients representing the four unconditional realizations in Fig. 1; (c) updates applied to unconditional samples to generate the conditional realizations.

(a)

(b)

(c)

Fig. 4. Data match and production forecasts for all 50 realizations at Producer # 2 (first row) and Producer # 3 (second row); Columns (a)–(c) show the well watercut, the oil production rate (STB/D) and the total oil produced (STB), respectively.

30

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

(b)

(a)

(c)

Fig. 5. (a) 25 Meandering fluvial channels used as prior unconditional models; (b) 25 K-SVD dictionary elements (with K = 1000 and S = 50); (c) 25 truncated SVD bases elements for S = 90;

Well Setup

Reference Log-Perm

(a)

(s1)

(s2)

(s3)

(s4)

Mean

Variance

(b)

(c)

(d)

Fig. 6. (a) Field setup (left) and reference log-permeability (right); (b) four unconditional log-permeability realizations (Columns 1–4), unconditional sample mean (Column 5) and variance (Column 6); (c) conditional log-permeability realizations from the SpRML method (Columns 1–4) corresponding to (b), and the conditional sample mean (Column 5) and variance (Column 6); (d) the conditional log-permeability realizations from the SVD-RML with S = 90 (Columns 1–4) corresponding to (b), along with the conditional sample mean (Column 5) and variance (Column 6). The sample mean and variance are calculated from 50 realizations.

Producer 1

Producer 2

Producer 3

(a)

(b)

Fig. 7. Data match and production forecasts for all 50 realizations at Producer # 1 (first column), Producer # 2 (second column) and Producer # 3 (third column); (a) and (b) show the well watercut for SpRML and SVD-RML methods respectively.

31

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

P0  P w ¼ Pc ðSw Þ Sw þ S0 ¼ 1

ð22Þ

to close the system of equations. The notation Pc represents the capillary pressure. The resulting equations can be solved numerically by discretizing the domain and specifying the initial and boundary conditions. In this paper, we assume a uniform (constant) initial pressure and oil saturation and consider no-flow boundary

Layer 1

Layer 2

Layer 3

condition at all boundary gridblocks. The boundary conditions are specified as source/sink terms in gridblock that contain injection and production wells. 3.1. Example 1: 2D multi-Gaussian case The two-dimensional model in this study is a 50  50 reservoir with heterogeneous permeability under a conventional 9-spot

Layer 4

P3

(a)

I2

P2 P4 P1

I3 I1

Layer 5 Injector Producer

I4

(b1)

(b2)

(b3)

(b4)

(c)

(d)

Fig. 8. (a) Reference log-permeability map for the PUNQ-S3 model; (b1)–(b4) four sample unconditional log-permeability realizations; unconditional sample mean (c) and variance (d) are calculated from 50 realizations. In all rows, Columns 1–5 show the five layers of the model.

32

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

waterflooding scenario. Fig. 2a shows the well configuration (left) and the reference log-permeability model (right). The total simulation time was 6 years, of which the simulated data from the first 4 years were used for model calibration while the data from the last two years were used to evaluate the predictive performance of the calibrated models. The reference log-permeability model and the unconditional samples were assumed to have the same

Layer 1

Layer 2

variogram (covariance) model. The geostatistical parameters used to generate the samples are shown in Table 1. The first four columns of Fig. 2b show four sample (out of 50) unconditional realizations of the log-permeability field. The fifth and sixth columns in Fig. 2b display the unconditional sample mean and variance, respectively. While the unconditional realizations have identical variograms, they generally exhibit different

Layer 3

Layer 4

Layer 5

(a)

(b1)

(b2)

(b3)

(b4)

(c)

(d)

Fig. 9. (a) Reference log-permeability map for the PUNQ-S3 model (repeated from Fig. 4); (b1)–(b4) four sample conditional log-permeability realizations corresponding to unconditional samples in Fig. 4; the conditional sample mean (c) and variance (d) are calculated from 50 realizations. In all rows, Columns 1–5 show the five layers of the model.

33

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37 Table 2 Geostatistical simulation parameters for PUNQ-S3 reservoir (Example 2). Property name

Layer 1

Layer 2

Layer 3

Layer 4

Layer 5

Grid dimensions (ft) Dx, Dy, Dz Kriging type Maximum conditioning data Variogram type Ranges max, med, min Angles azimuth, dip, rake Mean value (LogPerm) Standard Deviation (LogPerm)

180,180,1 Ordinary kriging 12 Spherical 1000, 285.7, 1 30, 0, 0 0.14 0.11

180,180,1 Ordinary kriging 12 Spherical 750, 750, 1 30, 0, 0 0.08 0.04

180,180,1 Ordinary kriging 12 Spherical 1500, 375, 1 45, 0, 0 0.14 0.11

180,180,1 Ordinary kriging 12 Spherical 750, 375, 1 150, 0, 0 0.10 0.06

180,180,1 Ordinary kriging 12 Spherical 1250, 416.6, 1 60, 0, 0 0.14 0.11

variability and connectivity in different regions in the reservoir. As a result, the unconditional sample mean does not show any particular trend in the log-permeability field. Fig. 2c displays the corresponding plots after conditioning with the SpRML. The conditional realizations show substantially reduced ensemble variability. The identified local connectivity patterns in the conditional samples (and sample mean) are consistent with the reference model. The main uncertainty in the variance map corresponds to the regions that are far from the wells where the observations are located. We note that the mean and variance are calculated from 50 realizations and are expected to have errors. This is especially true about the variance map. Fig. 3a shows 25 samples from the (2500) sparse K-SVD dictionary elements. In this case, the sparsity of the dictionary is S = 20. The linear expansion of the four unconditional realizations over the learned K-SVD dictionary (the expansion functions) results in the sparse representations shown in Fig. 3b. Fig. 3c depicts the updates applied to the set of sparse coefficients corresponding to the unconditional realizations to generate the conditional samples in the sparse domain (the spatial representations of the unconditional and conditional realizations are shown in Fig. 2b and c, respectively). These figures verify that the applied updates to the unconditional realizations are sparse, which is expected in SpRML implementation (Eq. (13)). Fig. 4 summarizes the data match (for the first four years) and future predictions (for the last two years) for the watercut, oil production rate and total oil production in Producer 2 (top row) and Producer 3 (bottom row). The figures clearly demonstrate the improvements that are obtained after calibration with the SpRML. It is noteworthy that even though in Producer 3 (second row of Fig. 4) water breakthrough occurs after 4 years (calibration period); the calibrated models can successfully forecast the water breakthrough time and the shape of the breakthrough function. 3.2. Example 2: 2D Non-Gaussian case (meandering channels) The second example is also a 2D reservoir with a more complex heterogeneity, representing a meandering fluvial channel formation. Unlike Example 1, the underlying geologic connectivity in this case is more complex and cannot be represented with the traditional variogram-based description methods. Furthermore, preserving such complex connectivity structures during model calibration is challenging. In this example, we aim to evaluate the performance of the K-SVD dictionary and the proposed sparse RML formulation for detecting such complex geologic connectivity by inverting the dynamic flow data. The reservoir in this case is discretized into a 100-by-100 grid system. The dictionary size is K = 1000 and the sparsity level is S = 50. Fig. 5a shows 25 sample permeability realizations. The corresponding 25 K-SVD dictionary elements and the leading 25 truncated SVD basis functions are shown in Fig. 5b and c, respectively. A total of S = 90 leading left singular vectors were used for the inversion with the SVD method. A clear distinction between the K-SVD and truncated SVD elements

in Fig. 5b and c is the appearance of the meandering channels in the K-SVD elements and the loss of the meandering features in the SVD basis. Considering that the training dataset for both K-SVD and SVD elements is identical, the two methods clearly show contrast in behavior. The K-SVD dictionary elements maintain a strong resemblance to the patterns in the prior library. During the sparsity promoting inversion, the most relevant features in the dictionary are picked and combined to reconstruct the solution. In the case of truncated SVD, the basis elements are generated by spectral decomposition of the features in the training data and lose the shape of the training patterns. During model calibration, the leading elements of the basis are combined to form the patterns in the solution. However, since the problem is ill-posed, the data required to retrieve the complex meandering channel is not available and the obtained solutions fail to represent the complex channel structure. Fig. 6a shows the well configuration for this example. A 13-spot well setup with 4 injectors and 9 producers is used to simulate a water-flooding scenario. The reservoir is initially filled with oil and water is injected to push the oil toward the producers and to maintain high reservoir pressure. The reference model is shown on the left side of Fig. 6a. Fig. 6b displays four sample initial permeability models (Columns 1–4) with the corresponding mean and variance maps (Columns 5–6). With a similar order to Fig. 6b, Figs. 6c and d display the resulting calibrated permeability maps using the K-SVD and SVD parameterization methods, respectively. The last two columns of Fig. 6c and d show the mean and variance maps for the calibrated realizations. It is evident from these results that the sparse K-SVD parameterization is better able to preserve the connectivity in the meandering channel whereas the SVD parameterization (with S = 90 basis elements) fails to reconstruct the meandering shape of the channel. This can be attributed to the second-order (covariance-based) nature of the SVD parameterization, and the flexibility of the sparsity-promoting algorithm for basis selection. We note that in the case of K-SVD, the calibration process uses all the K = 1000 elements with a sparsity promoting inversion while for calibration with the truncated SVD, the leading S = 90 singular vectors are included in a reduced-order parametric description. Finally, Fig. 7 shows the well watercut (WWCT) data match and prediction performance for the SpRML (Fig. 7a) and the SVD-based RML (Fig. 7b) methods for three different wells. In each case the first four years of data have been used for model calibration while the data from the last two years have been used for assessing the forecasts. A particularly important part in these plots is the water breakthrough time. In Producers 2 and 3, the watercut occurs during the last two year (prediction part) while in Producer 1 the water breakthrough happens during the first four year (model calibration part). Both SVD and K-SVD based methods can reproduce the early water breakthrough time in Producer 1. However, in Producers 2 and 3 (second and third columns of Fig. 7), where the water breakthrough time happens later during the prediction part, the predictions with the K-SVD based methods are evidently

34

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

(a)

(b) Unconditional Samples

(c) Updates Applied

(s1)

(s2)

(s3)

(s4)

Fig. 10. (a) 15 Sample K-SVD dictionary elements (out of K = 2660) with sparsity level S = 50; (b) sparse coefficients representing the four unconditional realizations of Fig. 4 in the K-SVD domain; (c) updates applied to the (sparse) unconditional realizations to generate the conditional samples.

superior to those from the SVD-based method. Overall, the results clearly indicate a better model calibration and prediction performance for the sparse RML approach. 3.3. Example 3: 3D PUNQ-S3 case We now consider the PUNQ-S3 3D benchmark model. The PUNQ-S3 case study is developed for uncertainty quantification in a real field setting. It is a small-size industrial reservoir model with 19  28  5 gridblocks, of which 1761 blocks are active. This

example contains geologic features that are more challenging than our previous example since it consists of channel-like permeability features is layers 1, 3, and 5 and low permeability in layers 2 and 4 that hinder vertical flow movement. The original model undergoes a water drive supported by a fairly strong aquifer to the north and west of the field. In this paper, we borrow the structure and reference permeability and porosity maps of the original PUNQ-S3, but study a waterflooding experiment with 4 injection and 4 production wells, as shown in Fig. 8a. The reservoir has five layers, each with a distinct heterogeneous permeability distribution.

35

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

(a)

(b)

(c)

Fig. 11. Data match (first four year) and production forecasts (last two years) for all 50 realizations at Producer # 1 (first row) and Producer # 3 (second row); Columns (a)–(c) show the well watercut, the oil production rate (STB/D) and the total oil produced (STB), respectively.

Table 3 Computation times for different examples in the paper. Example

Total computation time for K-SVD (min)

Computation time for 1 adjoint based forward simulation (min)

Total history matching time for one realization (min)

Example 1 multi-Gaussian Example 2 non-Gaussian Example 3: PUNQ-S3

568.5 828.5 189.9

0.91 1.89 0.99

18.3 29.1 30.3

The injection wells are completed in the bottom two layers while the producers are completed in the top three layers. The injectors are under rate control and the producers are under BHP (Bottom Hole Pressure) control. The injection BHP and production oil and water rates are collected and used as conditioning data in the first 4 years. The data collected from these wells during the fifth and sixth years are used to validate the prediction performance of the conditional realizations. We estimate the horizontal logpermeability field and derive the porosity and vertical log-permeability field from the relations log 10ðuXY Þ ¼ 9:02ðuÞ þ 0:77 and uZ ¼ 0:306  10uXY þ 3:124, respectively. These relations are borrowed from the original PUNQ-S3 model. Fig. 8b1–b4 shows four unconditional realizations that are generated from the prior model. Table 2 summarizes the geostatistical parameters used to generate the unconditional realizations. Fig. 8c and d display the sample mean and variance, respectively, generated from 50 unconditional realizations. The unconditional realizations and the corresponding sample mean do not contain the distinctive channel-like patterns that are observed in Layers 1, 3, and 5 of the reference model. Fig. 9 contains the conditional plots corresponding to Fig. 8. The plots in Fig. 9 show that conditioning with the SpRML algorithm results in a set of calibrated models that capture the main trends in the reference model. The main discrepancy in the calibrated models can be seen in the lower region of Layer 1. A closer examination reveals no observations from this region and unconditional samples all show high permeability for it. Hence, this discrepancy can be attributed to the bias in the unconditional realizations (prior model) and lack of observations to correct it. Overall, the calibrated log-permeability models capture the major trends in high and low permeability regions of the reference model fairly accurately. Fig. 10 contains the calibration results with the sparse (K-SVD) representation. Fig. 10a shows 15 (out of 2660) K-SVD dictionary elements (shown in 3D). The sparse representations of the four unconditional samples in Fig. 8b1–b4 using the K-SVD dictionary elements are shown in Fig. 10b. The sparsity level in this example is set to S = 50. This implies that only 50 non-zero coefficients (out of 2660) are needed to compactly approximate each unconditional sample in Fig. 10b. Fig. 10c contains the updates applied to the

unconditional realizations in Fig. 10b, to generate the conditional realizations. The sparse nature of these updates is consistent with the SpRML formulation (Eq. (13)). Finally, Fig. 11 depicts the production data match and forecasts with the unconditional and conditional log-permeability realizations. The plots of watercut and oil production rate, as well as total oil production for Producer 1 (first row) and Producer 3 (second row) are shown on the left, middle and right columns, respectively. It is clear from these figures that conditioning with the SpRML has led to substantial improvement in the data match and future production predictions. In terms of computation, the case studies included in this paper were run on a 64-bit machine with Intel core i7 CPU with clock time of 3.40 GHz. The total time for K-SVD dictionary generation for each case study is included in Table 3. The total time for running a forward flow simulation and the adjoint model for gradient computation, as well as the full inverse modeling time for one realization are also reported in Table 3. Construction of the K-SVD dictionary is computationally more demanding than generating the SVD basis. However, the dictionary construction is implemented off-line prior to solving the inverse problem. Liu and Jafarpour [38] discuss the complexity of the K-SVD algorithm and recommend reduced-order representation of the prior library to improve its computational complexity. The computational complexity of solving the inverse problem itself is comparable to that of standard gradient-based inversion methods.

4. Conclusion In this paper, we developed a RML implementation with sparse prior models that are not amenable to multi-Gaussian description. The SpRML provides a practical approach to formulate and solve sparse reconstruction inverse problems and to quantify the uncertainty in the solutions and future predictions. To this end, we first represent the reservoir property distributions in a sparse domain. Here, we use geologically-learned K-SVD dictionaries for this purpose. We then describe the sparse representations of the permeability models in the K-SVD domain with the Laplace prior PDFs. Combining the Laplace prior model and a Gaussian likelihood

36

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37

function to describe the observations, we derive the form of the posterior distribution that is not amenable to full analytical characterization. For nonlinear problems with multi-Gaussian priors, the RML has been proposed to approximately sample from the posterior distribution. Here, we follow the same approach and introduce the SpRML formulation as an approximate method to generate an ensemble of calibrated models for flow prediction and uncertainty quantification. While the RML sampling approach generally approximates the conditional realizations from the posterior distribution, it presents a practical model calibration and uncertainty quantification method. To evaluate its performance, we applied the SpRML method to three numerical experiments: (1) in a twodimensional model with a variogram-based multi-Gaussian permeability distribution; (2) in a two-dimensional model with a complex meandering fluvial channel connectivity generated using multiple point geostatistical simulation, i.e., not amenable to variogram-based representations, and (3) an adapted version of the three-dimensional PUNQ-S3 reservoir model. The results from these experiments suggest that the SpRML is a promising model calibration and uncertainty quantification method that can take advantage of sparse representations as a general attribute of spatially correlated reservoir properties. The SpRML algorithm is more advantageous than the deterministic sparse regularization techniques as it offers a systematic mechanism for including the uncertainty in the prior model and measurements, and for quantifying the resulting uncertainty in the estimated parameters and future predictions. Solving spatial inverse problems in geologically-inspired sparse domains transforms the model calibration into a low-order feature estimation problem, which is more consistent with the ill-posed nature of the problem. Sparse (low-rank) feature-based representation and estimation methods that take advantage of prior geologic connectivity information are better able to respect the expected geologic continuity during model calibration. They provide an effective framework for simplifying the description and estimation of complex (non-Gaussian) connectivity patterns by first applying sparse learning to encode such patterns and then recovering the solution by selecting and combining the relevant sparse features during the inversion. We introduced the SpRML as a practical uncertainty quantification approach for sparsity-promoting model calibration and uncertainty quantification under Laplace priors, and demonstrated its suitability using a series of numerical experiments for conditioning two-phase subsurface flow models on dynamic response data. In particular, our results show that, in comparison with the truncated SVD basis that is often used for reduced-order parameterization, learned (from prior geologic models) sparse dictionaries can be represented with Laplace distributions, which are more appropriate for sparse reconstruction and can offer superior performance in estimating complex non-Gaussian random fields.

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13] [14]

[15] [16]

[17] [18]

[19]

[20]

[21]

[22]

[23]

Acknowledgements This work is partially supported through a research grant awarded by the Department of Energy. The multiphase flow simulations were performed using Schlumberger’s Eclipse reservoir simulation software. The authors wish to thank the anonymous reviewers of this manuscript for their insightful comments and suggestions.

[24] [25]

[26]

[27]

References [28] [1] Aanonsen SI, Nævdal G, Oliver DS, Reynolds AC, Vallès B. The ensemble Kalman filter in reservoir engineering–a review. SPE J 2009;14(3):393–412. http:// dx.doi.org/10.2118/117274-PA. SPE-117274-PA. [2] Aharon M, Elad M, Bruckstein A. K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation. IEEE Trans Signal

[29]

Process 2006;54(11):4311–22 [see also Acoustics, Speech, and Signal Processing, IEEE Transactions on]. http://dx.doi.org/10.1109/TSP.2006.881199. Bhark E, Jafarpour B, Datta-Gupta A. A generalized grid-connectivity-based parameterization for subsurface flow data assimilation. Water Resour Res 2011;47:W06517. http://dx.doi.org/10.1029/2010WR009982. Britanak PC, Yip P, Rao K. Discrete cosine transform: general properties, fast algorithms, and integer approximation. Academic Press; 2006. http:// dx.doi.org/10.1016/B978-012373624-6/50003-5. Caers J, Zhang T. Multiple-point geostatistics: a quantitative vehicle for integrating geological analogs into multiple reservoir models. AAPG Mem 2004;80:383–94. Candès E, Wakin MB, Boyd SP. Enhancing sparsity by reweighted l1 minimization. J Fourier Anal Appl 2008;14(5):877–905. http://dx.doi.org/ 10.1007/s00041-008-9045-x. Carle SF, Fogg GE. Modeling spatial variability with one and multi-dimensional continuous Markov chains. Math Geol 1997;29(7):891–918. http://dx.doi.org/ 10.1023/A:1022303706942. Carrera J, Neuman SP. Estimation of aquifer parameters under transient and steady-state conditions, 1. Maximum likelihood method incorporating prior information. Water Resour Res 1986;22(2):199–210. http://dx.doi.org/ 10.1029/WR022i002p00199. Carrera J. State of the art of the inverse problem applied to the flow and solute transport problems. In: Groundwater flow and quality modeling, NATO ASI series, 1987. p. 549–585. http://dx.doi.org/10.1007/978-94-009-2889-3_31. Chavent G, Bissell R. Indicators for the refinement of parameterization. In: Tanaka M, Dulikravich GS, editors. Inverse problems in engineering mechanics. Elsevier; 1998. p. 309–14 (Proceedings of the third international symposium on inverse problems ISIP 98 held in Nagano, Japan). http:// dx.doi.org/10.1016/B978-008043319-6/50036-4. Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM J Sci Comput 1998;20:33–61. http://dx.doi.org/10.1137/ S1064827596304010. Chen Y, Zhang D. Data assimilation for transient flow in geologic formations via ensemble Kalman filter. Adv Water Resour 2006;29:1107–22. http:// dx.doi.org/10.1016/j.advwatres.2005.09.007. de Marsily G, Delhomme J, Delay F, Buoro A. 40 Years of inverse problems in hydrogeology. C R Acad Sci Ser II A Earth Planet Sci 1999;329(2):73–87. Deutsch CV, JournelL AG. GSLIB. Geostatistical software library and user’s guide. 2nd ed. Oxford, New York: Oxford University Press; 1998, ISBN 0 19 510015 8. Donoho D. Compressed sensing. IEEE Trans Inform Theory 2006;52(4):1289–306. http://dx.doi.org/10.1109/TIT.2006.871582. Elad M. Sparse and redundant representations: from theory to applications in signal and image processing. Springer; 2010. http://dx.doi.org/10.1109/ TIT.2010.2040892. Evensen G. Data assimilation: The ensemble kalman filter. New York: Springer-Verlag, Inc.; 2006, ISBN 3642037100. Franssen HH, Kinzelbach W. Real-time groundwater flow modeling with the ensemble Kalman filter: joint estimation of states and parameters and the filter inbreeding problem. Water Resour Res 2008;W09408. http://dx.doi.org/ 10.1029/2007WR006505. Gao G, Zafari M, Reynolds AC. Quantifying uncertainty for the PUNQ-S3 problem in a Bayesian setting with RML and EnKF, Paper SPE 93324 presented at the SPE reservoir simulation symposium, Houston, Texas, 31 January–2 February, 2005. http://dx.doi.org/10.2118/93324-PA Gavalas GR, Shah PC, Seinfeld JH. Reservoir history matching by Bayesian estimation. Soc Petrol Eng J 1976;16(6):337–50. http://dx.doi.org/10.2118/ 5740-PA. Gómez-Hernández JJ, Wen XH. To be or not to be multi-Gaussian? a reflection on stochastic hydrogeology. Adv Water Resour 1998;21(1):47–61. http:// dx.doi.org/10.1016/S0309-1708(96)00031-0. Guardiano FB, Srivastava RM. Multivariate geostatistics: beyond bivariate moments. In: Geostatistics tróia ’92. Quantitative geology and geostatistics, vol. 5, 1993. p. 133–44. http://dx.doi.org/10.1007/978-94-011-1739-5_12. Hendricks Franssen H-J, Alcolea A, Riva M, Bakr M, van der Wiel N, Stauffer F, Guadagnini A. A comparison of seven methods for the inverse modeling of groundwater flow and the characterisation of well catchments. Adv Water Resour 2009. http://dx.doi.org/10.1016/j.advwatres.2009.02.011. Jacquard P, Jain C. Permeability distribution from field pressure data. Soc Petrol Eng J 1965:281–94. http://dx.doi.org/10.2118/1307-PA. 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. http://dx.doi.org/10.2118/106453-PA. Jafarpour B, Goyal VK, Freeman WT, McLaughlin DB. Transform-domain sparsity regularization for inverse problems in geosciences. Geophysics 2009;74(5):R69–83. http://dx.doi.org/10.1190/1.3157250. 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;42(1):1–27. http://dx.doi.org/10.1007/s11004-009-9247-z. Jafarpour B. Wavelet reconstruction of geologic facies from nonlinear dynamic flow measurements. IEEE Trans Geosci Remote Sens 2010;49(5):1520–35. http://dx.doi.org/10.1109/TGRS.2010.2089464. 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. http://dx.doi.org/10.1029/2010WR009090.

M.M. Khaninezhad, B. Jafarpour / Advances in Water Resources 69 (2014) 23–37 [30] Jahns HO. A rapid method for obtaining a two-dimensional reservoir description from well pressure response data. Soc Petrol Eng J 1966;6(12):315–27. http://dx.doi.org/10.2118/1473-PA. [31] Khaninezhad MM, Jafarpour B, Li L. Sparse geologic dictionaries for subsurface flow model calibration: part I. Inversion formulation. Adv Water Resour 2012;39(0):106–21. http://dx.doi.org/10.1016/j.advwatres.2011.09.002. [32] Khaninezhad MM, Jafarpour B, Li L. Sparse geologic dictionaries for subsurface flow model calibration: part II. Robustness to uncertainty. Adv Water Resour 2012;39(0):122–36. http://dx.doi.org/10.1016/j.advwatres.2011.10.005. [33] Kitadinis P. Quasi-linear geostatistical theory for inversing. Water Resour Res 1995;31(10):2411–9. http://dx.doi.org/10.1029/95WR01945. [34] Laloy E, Rogiers B, Vrugt JA, Mallants D, Jacques D. Efficient posterior exploration of a high-dimensional groundwater model from two-stage Markov Chain Monte Carlo simulation and polynomial chaos expansion. Water Resour Res 2013;49. http://dx.doi.org/10.1002/wrcr.20226. [35] Lawson CL, Hanson RJ. Solving Least Squares Problems. Society for Industrial and Applied Mathematics; 1995. http://dx.doi.org/10.1137/1.9781611971217. [36] Li L, Jafarpour B. Effective solution of nonlinear subsurface flow inverse problems in sparse bases. Inverse Prob 2010;26(10):105016. http://dx.doi.org/ 10.1088/0266-5611/26/10/105016. [37] Li L, Jafarpour B. A sparse Bayesian framework for conditioning uncertain geologic models to nonlinear flow measurements. Adv Water Resour 2010;33(9):1024–42. http://dx.doi.org/10.1016/j.advwatres.2010.06.005. [38] Liu E, Jafarpour B. Learning sparse geologic dictionaries from low-rank representations of facies connectivity for flow model calibration. Water Resour Res 2013;49(10):7088–101. http://dx.doi.org/10.1002/wrcr.20545. [39] Mallat S. A wavelet tour of signal processing: the sparse way. 3rd ed. Academic Press; 2003. [40] Nowak W. Best unbiased ensemble linearization and the quasi-linear kalman ensemble generator. Water Resour Res 2009:W04431. http://dx.doi.org/ 10.1029/2008WR007328.

37

[41] Oliver D, Reynolds A, Liu N. Inverse theory for petroleum reservoir characterization and history matching. Cambridge University Press; 2008. http://dx.doi.org/10.1017/CBO9780511535642. [42] Oliver DS, Chen Y. Recent progress on reservoir history matching: a review. Comput Geosci 2011;15(1):185–221. http://dx.doi.org/10.1007/s10596-0109194-2. [43] Oliver DS, He N, Reynolds AC. Conditioning permeability fields to pressure data, Paper presented at the 5th European conference for the mathematics of oil recovery, 3–6 September, Leoben, Austria: 1996. http://dx.doi.org/10.1007/ BF02769620. [44] Shah PC, Gavalas GR, Seinfeld JH. Error analysis in history matching: the optimal level of parametrization. Soc Petrol Eng J 1978:219–28. http:// dx.doi.org/10.2118/6508-PA. [45] Strebelle S, Journel AG. Reservoir modeling using multiple-point statistics, Paper SPE 71324 presented at the 2001 SPE annual technical conference and exhibition, 30 September–3 October, New Orleans: 2001. http://dx.doi.org/ 10.2118/71324-MS. [46] Tarantola A. Inverse problem theory and methods for model parameter estimation. Philadelphia, PA: SIAM; 2005. http://dx.doi.org/10.1137/ 1.9780898717921. [47] Tikhonov AN, Arsenin VI. Solution of ill-posed problems. Washington, DC: Winston; 1977. http://dx.doi.org/10.1007/978-94-015-8480-7. [48] Yeh W. Review of parameter estimation procedures in groundwater hydrology: the inverse problem. Water Resour Res 1986;22:95–108. http:// dx.doi.org/10.1029/WR022i002p00095. [49] Wen X, Chen W. Real-time reservoir model updating using ensemble Kalman filter. SPE J 2006;11(4):431–42. http://dx.doi.org/10.2118/111571-PA. [50] 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. http://dx.doi.org/10.1029/2001WR001146.