In contrast to the widespread use of Monte Carlo simulations (MCS), stochastic equations have shown remarkable potential to develop efficient, accurate and physical insightful models. The main challenge of solving stochastic equations is to develop tractable descriptions of the system response in terms of stochastic differential operators and random fields. In particular, the solution of stochastic PDEs in porous media flow raises unexplored challenges to the solution of the discretized system of equations due to the significant size of standard reservoir problems combined with the uncertainty induced by error measurements, discontinuities, nonlinearities in the reservoir parameters at different spatial and temporal scales.

This work introduces a Krylov-Karhunen-Loeve moment equation (KKLME) approach for the solution of stochastic PDEs arising in large-scale porous media flow applications. The approach combines recent developments in Karhunen-Loeve moment equation methods with a block deflated Krylov iterative solution of a sequence of deterministic linear algebraic equations sharing the same matrix operator but different right-hand sides (RHSs).

In this approach, the log of the random field (i.e., log of permeability) describing the transmissibility coefficients is decomposed as the sum of a deterministic average log field plus a mean-zero random fluctuation. The covariance function associated with the fluctuation component is further decomposed by means of the Karhunen-Loeve (KL) expansion scheme. This expansion consists of modes (i.e. stochastic orders) of increasing frequency but decreasing magnitude. It has been shown that the KL expansion is of mean square convergence and summation thus implying significant computational savings. Subsequently, a mixed finite element procedure (equivalent to a cell-centered finite difference scheme under a suitable quadrature rule) is employed to derive a system of linear random algebraic equations.

In order to compute higher-order approximations for the different pressure moments, a perturbation approach followed by an expansion of orthogonal random variables is performed to express the variability of pressures with respect to the random field. The algebraic manipulation of modes and moments results in a sequence of deterministic linear systems with multiple RHSs sharing the same matrix operator. This matrix operator corresponds to the discretization of the average random field that, in general, has better algebraic properties than the operator associated with the original random field. A set of independent RHSs becomes available when a lower moment is computed, and each moment involves the solution of several modes or combinations of them with previous moment solutions. Since the associated average system is generally large and sparse, a Krylov subspace iterative solver is suitable in this case. The availability of RHSs describes a particular computational pattern that is amenable for highly scalable implementations and, at the same time, imposes particular challenges for an efficient Krylov subspace implementation.

With regard to the solution of systems with multiple RHSs, significant advances have been made in recycling the information generated by Krylov iterative methods. Deflation methods have been shown to be very effective in these circumstances as they ``remove'' components of the solution associated with extreme eigenvalues that prevent or slow down convergence. To perform deflation the Krylov subspace is selectively augmented with approximate eigenvectors that are either computed during the iteration or somehow known a priori from some geometrical/physical mean.

On the other hand, in order to reduce the negative effects of sequential inner products on memory performance, it is advisable to rely on block implementations or, in other words, process a subset of right-hand-sides in a simultaneous fashion. A novel block deflation Krylov iterative method is proposed. The main feature of this implementation is to enable changes both in the size of the block and in the augmented subspace as the iteration proceeds. This ensures both numerical stability and flexibility in replacing deflated RHSs (i.e., associated with converging solutions) with new RHSs.

Numerical experiments show that the KKLME approach requires significant less computing time than MCS to converge to the different statistical moments for the pressure response. Results are shown for single phase flow and for pressure system arising in two-phase fully implicit formulations using a block deflation CG iterative method. In light of these preliminary results, extensions to relate random fields with spectrum information are required in order to exploit further efficiencies.