{'fraction_non_alphanumeric': 0.034481218455491534, 'fraction_numerical': 0.011523516012327482, 'mean_word_length': 4.454323995127893, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Conditioning by Projection for the Sampling from Prior Gaussian Distributions

In this work we are interested in the (ill-posed) inverse problem for absolute permeability characterization that arises in predictive modeling of porous media flows. We consider a Bayesian statistical framework with a preconditioned Markov Chain Monte Carlo (MCMC) algorithm for the solution of the inverse problem. Reduction of uncertainty can be accomplished by incorporating measurements at sparse locations (static data) in the prior distribution. We present a new method to condition Gaussian fields (the log of permeability fields) to available sparse measurements. A truncated Karhunen-Lo\`eve expansion (KLE) is used for dimension reduction. In the proposed method the imposition of static data is made through the projection of a sample (expressed as a vector of independent, identically distributed normal random variables) onto the nullspace of a data matrix, that is defined in terms of the KLE. The numerical implementation of the proposed method is straightforward. Through numerical experiments for a model of second-order elliptic equation, we show that the proposed method in multi-chain studies converges much faster than the MCMC method without conditioning. These studies indicate the importance of conditioning in accelerating the MCMC convergence.


Introduction
Subsurface formations have spatial variability in their hydraulic properties in multiple length scales. It has been established that such variability plays an important role in determining fluid flow patterns in heterogeneous formations [24,25,18,4]. The characterization of such formations is essential for predictive simulations of multiphase flows in problems such as oil recovery, CO 2 sequestration, or contaminant transport. Procedures for characterization, such as a Bayesian framework using the Markov chain Monte Carlo (MCMC) algorithm, entails the solution of an ill-posed inverse problem. From static and dynamic data, the unknown coefficients of a governing system of partial differential equations that model fluid flow have to be recovered. In the inversion process, the available static data (for instance, measured absolute permeability values at sparse locations) have to be incorporated in the sampling of the (unknown) permeability field so that uncertainty is reduced as much as possible. This way, a conditioning procedure has to be applied such that all samples considered by the inversion method honor available measurements.
MCMC methods have been applied in several areas of science and engineering, and have attracted the attention of many research groups [19,33]. There are many recent and important developments related to gradient-based MCMC procedures. We mention, among many others, the Hamiltonian Monte Carlo [8], and the active subspace methods [10]. There are also Hessian-based procedures [34]. However, for subsurface flow problems gradient and Hessian calculations are computationally very expensive. Thus, for the problems we have in mind we consider the preconditioned (or two-stage) MCMC introduced in [16,9]. Some developments of this method include the multiphysics version [23], the parallelization in multi-core devices [21] and its extension to deformable subsurface formations [5]. This procedure has also been successfully applied in geophysics [44]. The Bayesian framework for inversion that we will use in this work has been carefully discussed in [1].
Our focus in this work is the conditioning of samples to sparse measured data. The importance of conditioning in reservoir simulation has been recognized for many years (see, for instance, [26]). The sequential Gaussian simulation (SGS) [13] has been widely used for this purpose. More recently a truncated Karhunen-Loève expansion (KLE) [29,16,22] has been applied both for dimensional reduction of the (computationally infeasible) large stochastic dimensions of fine grid discretizations of subsurface flow problems as well as for conditioning. In [37] one can find a discussion of the advantages of the KLE over the SGS. There are two lines of work for the conditioning of samples using the KLE. In [30] the basic KLE construction is altered to incorporate static data. In another line of work the vector of independent, identically distributed (i.i.d.) N (0, 1) Gaussian variables that enter the KLE construction is modified to take data into account while keeping the KLE eigenvectors unaltered. This approach can be seen in [36] where a search through a set of matrices of size determined by the number of measurements is performed, and they take the one with the best condition number. In [37] the authors use the expressions for conditional means and covariances of Gaussian random variables [46] that require the inversion of a data covariance matrix.
The method that we introduce here does not alter the basic KLE construction and is computationally very inexpensive. By using the measured data in the KLE we define a data matrix that models Gaussian perturbations on top of a surface defined by kriging [12] of the available data. The proposed procedure consists of replacing the vectors of i.i.d., N (0, 1) Gaussian variables by the closest vector to them (in the least-squares sense) in the nullspace of the data matrix. This way, the Gaussian fields produced by the KLE honor exactly all the data available. We remark that all one needs to implement the proposed method is a basis for the nullspace of the data matrix. Such a basis is used to construct the projection onto the nullspace of the data matrix. Besides introducing this new method, we also present multi-chain MCMC studies for a model elliptic equation. We show convergence of this method for large dimensional problems. Moreover, these studies are used to illustrate the importance of conditioning for the convergence of MCMC methods. Although outside the scope of this work we believe the method presented here can also be applied to other sampling strategies such as the hierarchical method discussed in [6]. This topic is currently being considered by the authors.
This work is organized as follows. In section 2 we discuss the model problem, kriging and dimensional reduction. Section 3 introduces the new method for conditioning by projection. In section 4 we explain the Bayesian framework for inversion. We describe the preconditioned MCMC method with and without conditioning for the construction of permeability fields. Also, the method we use for convergence assessment of MCMC methods is presented. Section 5 is dedicated to a discussion of numerical results from our experiments. Conclusions and some remarks appear in section 6.

The Model Problem
In the modeling of incompressible subsurface flow problems, such as in two-phase flows in oil reservoirs [14,38], the system of governing partial differential equations consists of a second order elliptic equation coupled to a hyperbolic dominated equation that describes fluid transport. In this work, we illustrate the proposed method for conditioning Gaussian fields in terms of the elliptic equation that enters in this system.
Let us consider Ω ⊂ R 2 , a bounded domain with Lipschitz boundary, let v(x) and p(x) denote the Darcy velocity and the pressure of the fluid, respectively, where x ∈ Ω, and k(x) is the absolute permeability. Using Darcy's law, we have the following model elliptic problem: v where f ∈ L 2 (Ω) is a given source term, together with appropriate Dirichlet and Neumann type boundary conditions p = g p ∈ H 1 2 (∂Ω p ) and v ·n = g v ∈ H − 1 2 (∂Ω v ), the pressure boundary data and normal velocity data, respectively, with ∂Ω = ∂Ω p ∪ ∂Ω v and ∂Ω p ∩ ∂Ω v = ∅,n is the exterior unit normal vector. The numerical approximation of the above elliptic problem is derived from the weak formulation of the velocitypressure system. The weak formulation can be found in [17].

Kriging
Kriging is a type of interpolation method, which uses a local estimation technique that gives the best linear unbiased estimator of unknown values [12,27]. The main idea of kriging is to estimate the values of a function (or random field) in a domain by calculating a weighted average of a given point in some neighborhood. The kriging takes into account (i) the distances between the estimated points in the domain and the given data points, and (ii) the distances between the given data points themselves. Thus, the kriging gives more weights to the nearest data points. When a point x 0 coincides with a given data point x i , the kriging gives the exact value and the variance is zero. In the following subsection, we discuss the KLE, and in Section 3, we combine both kriging and KLE.

Dimensional Reduction
We use the KLE to reduce the large dimensional uncertainty space describing the permeability field k(x). In this subsection, we present a brief description of the KLE [29,20]. Let x ∈ Ω, and suppose log [k( can be written as the following: (2) can be written as where λ i and ϕ i satisfy equation (3). We arrange the eigenvalues in the descending order. The infinite series in equation (4) is called the Karhunen-Loève expansion. In general, we only need first n dominating terms such that the energy E is more than 95% [28], i.e.
Then, a truncated KLE is defined by From now on we will use Y (x) instead of Y k n (x) for samples generated by the truncated KLE.

Conditioning by Projection
Suppose the Gaussian field Y (x) defined in equation (6) is a Gaussian perturbation on top of a kriged fieldŶ (x) discussed in section 2.2. That is, This can also be written as where D is a diagonal matrix of the dominating n eigenvalues and θ θ θ = (θ 1 , θ 2 , . . . , θ n ).
are known. Thus, by substituting (9) in equation (8) we have This honors m data values at the corresponding sparse locationsx. We can rewrite (10) as follows The solution of this homogeneous system of equations (11) is the nullspace N (A), a subspace of R n . Given an arbitrary vector θ θ θ ∈ R n of i.i.d., N (0, 1) Gaussian variables that would be used in the KLE (7), in order to honor the given m data we select the closest vector to θ θ θ in the subspace N (A). This vector is obtained by projecting θ θ θ onto the nullspace of the data matrix A.
Assume that A has rank r, we can easily find an orthonormal basis β = {q 1 , q 2 , . . . , q n−r } for the nullspace of A. We assume that r < n, so that A has a non-trivial nullspace. The projection matrix P onto the nullspace of A is given by [43] where Q is a matrix that has q 1 , q 2 , . . . , q n−r as columns. The projection vector of θ θ θ onto the nullspace of the data matrix A is given bŷ   1 illustrates the projection of the i.i.d., N (0, 1) Gaussian variables θ θ θ vector onto the nullspace of the data matrix A. The error vector defined as e e e = θ θ θ −θ θ θ is minimized in the least squares sense. It is perpendicular to the projection vectorθ θ θ. Finally, we can rewrite equation (7) as whereθ i , i = 1, . . . , n are the components ofθ θ θ in equation (13). If x =x i , i = 1, 2, . . . , m, in equation (14), the second term on the right-hand side is zero. This leads to Y (x i ) =Ŷ (x i ), and thus, we honor the exact values of the permeability field at the known locationsx 1 , . . . ,x m .

Posterior Exploration
Our primary interest in this article is to characterize the permeability field of our domain of study. Here, we discuss the characterization of the permeability field conditioned on available pressure data. We use half (in a chessboard pattern) of the total pressure data in the discretized domain [45]. We denote the (log of the) permeability field and the reference pressure data by η η η and R p , respectively. A Bayesian statistical approach along with a preconditioned MCMC method with and without conditioning is used to solve the inverse problem. The posterior probability is calculated using the Bayes' theorem with respect to R p as follows: where P (η η η) represents the Gaussian prior distribution. The normalizing constant is disregarded in this investigation as we do an iterative search in our MCMC method. We construct the permeability field η η η(θ θ θ) using the KLE, where θ θ θ is generated by the MCMC method. The likelihood is considered to be a Gaussian function as in [16], and as a result of that we have our likelihood function as where R η η η denotes the simulated pressure data. We set the covariance matrix Σ as Σ = I I I/2σ 2 R , where I I I and σ 2 R represents the identity matrix and the precision parameter, respectively.
Below we describe the preconditioned MCMC method with and without conditioning.

MCMC algorithms with and without conditioning
In our study, we use the MCMC method with and without conditioning. The following random walk sampler (RWS) [11] is used for both MCMC cases: where θ p denotes the current proposal and θ represents the previously accepted proposal. The symbol β denotes the tuning parameter and represents N (0, 1)-random variable. At η η η(θ θ θ)generate θ θ θ p using equation (18).

7:
Solve the forward problem on the coarse-scale to get R c .

9:
if η η η p is accepted then 10: Use η η η p in the fine-scale simulation to get R f .

12:
if η η η p is accepted then η η η = η η η p . 13: end if 14: end if 15: end for We first discuss the algorithm of the MCMC method without conditioning [16,9]. The filtering step of this method is based on a coarse-scale model approximation of the governing equation (1). The coarse-scale discretization is similar to the fine-scale discretization and the permeability field η η η(θ θ θ) is projected on the coarse-scale. An upscaling procedure [15] is used to set an effective permeability field on a coarse-grid that provides a similar average response as that of the underlying fine-scale problem. The numerical simulator is run on the coarse-scale model and produces the coarse-grid pressure field R c . The coarse-scale and fine-scale acceptance probabilities are estimated as α c (η η η, η η η p ) = min 1, I(η η η|η η η p )P c (η η η p |R p ) I(η η η p |η η η)P c (η η η|R p ) , and where P c and P f represent the posterior probabilities computed at coarse-and finescale, respectively. In the MCMC method with conditioning, we construct the permeability field η η η(θ θ θ) using equation (14). The algorithm of the MCMC method without conditioning is described in [2] and the algorithm of the MCMC method with conditioning is presented in Algorithm 1.

Convergence Assessment
In this subsection, we discuss the convergence of the MCMC methods for a reliable characterization of the permeability field. Determining the stopping criterion in the MCMC simulations is a challenging task. For this reason, researchers use several MCMC convergence diagnostics to decide when it is safe to stop the simulations. Different diagnostic measures are discussed in [40,31,35,41]. A common approach is to start more than one MCMC chain using distinct initial conditions. In this paper, we study the convergence of the proposed MCMC methods using the Potential Scale Reduction Factor (PSRF) and its multivariate extension MPSRF [7]. Note that the PSRF does not take all the parameters into account and thus may not provide a very reliable indicator of the convergence of the parallel chains. However, the MPSRF takes all the parameters and their interactions into consideration. Therefore, the MPSRF is a better measure to determine the convergence of several chains that explore a high dimensional parameter space. Let us consider l posterior draws of θ (dimension n) for each of the k > 1 parallel chains. Let θ c j be the generated samples at iteration c in jth chain of k parallel MCMCs. Then, the posterior variance-covariance matrix is computed by The within-chain covariance matrix W and between-chain covariance matrix B are defined by and respectively.θ j. andθ .. represent the mean within the chain and the mean of the k combined chains, respectively, and T stands for transpose. The PSRFs are computed by The maximum of PSRF values, which is closer to 1, indicates that MCMC proposals sample from approximately the same posterior distribution. The MPSRF is defined as [32] where λ is the largest eigenvalue of the positive definite matrix W −1 B/l. If the proposals in k chains are sampled from almost the same distribution, the between chain covariance matrix becomes zero. Consequently, λ and the MPSRF approach to 0, and 1, respectively, for sufficiently large l. Thus, we achieve the required convergence of the chains.

Setup of the problem
In this subsection, we present the details of the simulation of our problem of interest. We consider a unit square-shaped physical domain and solve the elliptic equation (1) in the domain. The simulated elliptic solutions are used in the MCMC algorithms. We then present a comparative study between the MCMC methods with and without conditioning. In the MCMC method with conditioning, we include the prior measurement of the permeability field at nine specific sparse locations. This permeability field is constructed using the KLE in which we consider the covariance function as where L x = 0.4 and L y = 0.8 are the correlation lengths and σ 2 Y = Var[(Y k ) 2 ] = 1. Fig. 3 illustrates the decay of the eigenvalues for these values. Considering the fast decay of the eigenvalues, in this paper, we take the first twenty dominating eigenvalues in the KLE, which preserves almost 100% (using the equation (5)) of the total variance or energy of Y in the equation (6). We set the source term f = 0 and impose Dirichlet boundary conditions, p = 1 and p = 0, on the left and right boundaries, respectively. We set a Neumann-type boundary condition to zero (i.e., no-flow) everywhere else on the domain boundaries. In our study, we use a synthetic reference permeability field that is generated on a fine-grid of size 16 × 16. We then run the forward problem using the numerical simulator and generate the corresponding reference pressure field. See Fig. 2 for both fields. We take the tuning parameter β = 0.85 in (18) for both MCMC methods. As a screening step in the MCMC algorithm, we introduce a coarse-scale filter of size 8 × 8. The coarse-scale simulation runs almost two times faster than the corresponding finescale simulation while maintaining the general trend of the fine-scale simulation. We change one stochastic parameter at every iteration to construct a new permeability field for both MCMC studies.

Convergence Analysis of the MCMC Method
We analyze here the convergence of the MCMC algorithm for both studies in terms of the maximum of PSRFs and MPSRF. We run four MCMCs starting from distinct initial seeds for the MCMC method without conditioning. We also run another four chains, but using the same corresponding initial seeds in the former runs, for the MCMC method with conditioning. In Table 1 we show the acceptance rates for these two MCMC studies. We take σ 2 F = 10 −4 and σ 2 C = 5 × 10 −3 for coarse-and fine-scale precisions, respectively, in equation (16).
MCMC without conditioning MCMC with conditioning σ 2 We observe from the Table 1 that the MCMC method with conditioning has a better acceptance rate than the MCMC method without conditioning. In the calculation of the PSRFs and MPSRF, we take into account the θ θ θ values that are accepted on the fine-scale.
If there is no change in the theta values from one state to another in the Markov chain, the repetition of the accepted theta values on the fine-scale is considered in the calculation. Note that a MCMC method converges to the target distribution if both PSRFs and MPSRF values are below 1.2 [42]. Fig. 4 depicts the maximum of the PSRFs and MP-SRF for both cases with the same number of proposals, 26000, after burn-in. In each chain, 6500 proposals are considered. We observe that the MCMC with conditioning achieves a full convergence. Towards the end of the curves, the maximum PSRF and MPSRF values are 1.14 and 1.16, respectively. Both PSRF and MPSRF curves for the MCMC method without conditioning are considerably away from the numerical value of 1.2. Thus, we conclude that the MCMC method with conditioning explores the posterior distribution much faster than the MCMC method without conditioning. Moreover, Fig. 4 demonstrates that the maximum of PSRFs is bounded above by the MPSRF as shown in [7]. In the next paragraph, we discuss the importance of this convergence in order to recover the underlying permeability field.  Now, we discuss the permeability fields that are accepted in the MCMC methods. In Fig. 5 we show twelve permeability fields (three per chain) at 40, 5000, and 10000 iterations, respectively, in four chains for the MCMC method without conditioning. In Fig. 6 we present the corresponding permeability fields in four chains obtained from the MCMC method with conditioning. The process of constructing the permeability fields for this MCMC method with conditioning has two steps: In the first step, we create a permeability field using the KLE in which we provide the projected proposal as an input. Then we add the kriging data obtained from the measurements at nine locations (See section 3). We observe that the MCMC method with conditioning recovers the reference permeability field in Fig. 2 much better than the MCMC method without conditioning. This result is consistent with the MCMC convergence discussed in the above paragraph. Thus, we conclude that the MCMC method with conditioning shows a better performance for the characterization of the permeability field in our problem.

Conclusions
In this work, we considered the use of MCMC methods for the characterizations of subsurface formations. Usually, sparse measurements of quantities of interest (such as the permeability field) are available in several field locations and they have to be incorporated in the characterization procedure to reduce uncertainty as much as possible.
We introduced a computationally efficient method for conditioning (the log of) permeability fields. The method is based on a projection onto the nullspace of a data matrix defined in terms of a KL expansion. We also presented multi-chain studies that illustrate the importance of conditioning in accelerating the MCMC convergence.
The authors and their collaborators are currently applying the Multiscale Perturbation Method [3,39] to further speed up MCMC studies.