ECMOR XIV - 14th European Conference on the Mathematics of Oil Recovery

Numerical simulation of multiphase flow in porous geological formations is widely used for subsurface flow management, including oil and gas production. Large-scale simulation is often limited by the computational speed, where nonlinear convergence is one of the main bottlenecks. The nonlinearity of flow properties, the coupling of driving forces for fluid migration, and the heterogeneities of the formation are three main causes for convergence difficulties. Here, we analyze the nonlinearity of two-phase transport in porous media, and we propose an efficient nonlinear solver based on this understanding.

We focus on immiscible, incompressible, two-phase transport in the presence of viscous and buoyancy forces. We investigate the nonlinearity of the discrete transport equation obtained from finite-volume discretization with Single-Point Upstream weighting (SPU), which is the industry standard. In particular, we study the discretized numerical flux expressed as a function of the upstream and downstream saturations of the total velocity. We analyze the locations and complexity of the unit-flux, zero-flux, and inflection lines of the numerical-flux saturation space. The unit- and zero-flux lines, referred to as kinks, correspond to a switch in the flow directions of the different phases, and if SPU is used, then the numerical flux is not differentiable at those points. These kinks and inflection lines are major sources of nonlinear convergence difficulties, especially when their locations in the numerical flux depend on both saturations in the upstream and downstream cells. Our analysis of the discrete (numerical) flux offers a theoretical basis of the convergence challenges associated with multi-cell problems and serves as a foundation for developing efficient nonlinear solvers.

With this understanding, we propose a nonlinear solution scheme that is a significant refinement of the works of Jenny et al. (2009) and Wang and Tchelepi (2013) . We divide the flux function into saturation ‘trust regions’ delineated by the kinks and inflection lines. Determining the boundary of each trust region is straightforward, and it only needs to be computed once in a preprocessing step before performing a simulation. The Newton updates are performed such that two successive iterations do not cross any trust-region boundary. If a crossing is detected, the saturation value is chopped back to the boundary. Our saturation chopping algorithms captures the inflection lines of the numerical flux much more accurately than the treatment of Wang and Tchelepi (2013) . The nonlinear convergence behavior is analyzed using numerical examples, and significant improvements over existing trust-region nonlinear solvers are demonstrated.

Fully Implicit discretization of flow and transport equations gives rise to a system of coupled nonlinear equations that is typically solved using standard Newton method. For a given timestep size, even if the Newton-based iterative procedure converges, the cost associated with updating all the unknowns simultaneously can be quite expensive. Conventional sequential-implicit strategies can be used to reduce the cost, but they suffer from severe restrictions on the allowable timestep size.

In this paper we formulate, verify and analyze the computational efficiency of a new nonlinear solution strategy. The crux of the proposed algorithm is the use of a hybrid strategy to treat the co-current and counter-current flow regions differently. At each Newton iteration, we first update the pressure variables by solving the Schur complement reduced form of the equations. To avoid costly computation of the reduced Jacobian, we employ phase-based potential ordering, giving rise to a lower triangular Jacobian that can be used to compute the reduced Jacobian efficiently. After updating pressures, we decompose the domain into components in a way that can be sorted and traversed from upstream to downstream. A component in the co-current flow regions is made up of a single cell, whereas a component made up of multiple connected cells indicates counter-current flow in which pressure and saturation are tightly coupled. Marching down from upstream to downstream, for a single-cell component we update saturation nonlinearly by solving the scalar saturation equation(s). For a multi-cell component we discard the pressure of the corresponding cells obtained in the first step, and then perform a simultaneous linear update of the saturation and pressure variables by solving the local linear system. Once all of the components are visited and updated, the iteration is over.

We present a variety of challenging numerical examples in 2-D and 3-D in the presence of strong gravity and heterogeneity. Our results show that as compared with standard Newton method, the proposed hybrid solver has a lower computational cost per iteration without compromising the allowable timestep size. The computational efficiency gain depends on the number and size of the components which vary over the course of iterations. At best, pressure and saturation updates are fully decoupled and saturation variables are updated sequentially one at a time. We present a rigorous complexity analysis of the algorithm for linear solvers with arbitrary order of complexity.

The scalability of the Algebraic Multiscale Solver (AMS) ( Wang et al., 2014 ) for the heterogeneous pressure system that arises from incompressible flow in porous media is analyzed and experimentally demonstrated in parallel computing environments. The major steps of AMS are highly parallel, but the solver overall scalability is strongly tied to the choice of parameters and algorithms involved in each step. These choices additionally impact the convergence properties of the solver. The balance between computational scalability and convergence rate is carefully considered, to ensure high overall performance while maintaining robustness.

The basis-function kernel, which dominates the setup phase, and the local smoother, which dominates the solution phase, are studied in detail. In addition, the balance between convergence rate and scalability as a function of the coarsening ratio, Cr, is investigated.

Based on this analysis, the performance of a scalable AMS implementation is tested using highly heterogeneous problems derived from the SPE10 benchmark ( Christie et al., 2001 ) and geostatistically generated. The problems range in size from a few million to more than a 100 million cells. The solver robustness and scalability is demonstrated on modern multi-core systems and compared with the widely used parallel SAMG solver ( Stüben, 2012 ).

Locality is inherent to all transient flow and transport phenomena. Moreover, the superposition of the two disparate spatiotemporal scales that underlie flow and transport leads to a problem of dimensionality. Numerous Adaptive discretization methods have been devised to exploit an a priori understanding of locality. While such methods have provided various degrees of success, they are fundamentally restricted by the fidelity of the discretization under aggressive adaptivity. This work seeks a novel class of nonlinear solvers which are proposed to adapt the level of computation to precisely match that of the underlying spatiotemporal change, without affecting the underlying discretization model. We devise a class of nonlinear iteration that on the one hand, converges as rapidly as the best available safeguarding method (e. g. trust-region), while on the other, performing a number of operations that is at most equal to the number of cells that experience changes over an iteration.

We achieve this by developing an Asynchronous Multirate numerical integration of the Newton Flow differential system. At each nonlinear iteration, the domain of interest is partitioned adaptively in two or more levels of disjoint subdomains on the basis of the predicted rate of change of state variables. On one extreme, there are only two levels of partitions; cells that will experience a nonzero change, and cells will not. This two-level solver is equivalent to the adaptively localized solution of the linear Jacobian system. On the other extreme, the domain is decomposed into multiple disjoint subdomains. Under this scheme, the subdomains are solved sequentially using a dynamic partitioning strategy in the order from fastest to slowest. We present detailed computational results focused on general multiphase flow models. The performance improvement directly depends on the extent of locality present in the model. On the two-level end of the spectrum, the convergence rate of the proposed method is unadulterated while the performance is improved by an order of magnitude in computational time.

On the multilevel end of the spectrum, while additional performance gains are obtained for transport components, the convergence rate for pressure requires more costly synchronization strategies. This method looks very promising for the simulation of extremely complex models where well controls change dynamically.

An adaptive Algebraic Multiscale Solver for Compressible (C-AMS) flow in heterogeneous oil reservoirs is developed. Based on the recently developed AMS [ Wang et al., 2014 ] for incompressible linear flows, the C-AMS extends the algebraic formulation of the multiscale methods for compressible (nonlinear) flows. Several types of basis functions (incompressible and compressible with and without accumulation) are considered to construct the prolongation operator. As for the restriction operator, C-AMS allows for both MSFV and MSFE methods. Furthermore, to resolve high-frequency errors, Correction Functions and ILU(0) are considered. The best C-AMS procedure is determined among these various strategies, on the basis of the CPU time for three-dimensional heterogeneous problems. The C-AMS is adaptive in all aspects of prolongation, restriction, and conservative reconstruction operators for time-dependent compressible flow problems. In addition, it is also adaptive in terms of linear-system update. Though the C-AMS is a conservative multiscale solver (i.e., only a few iterations are employed infrequently in order to maintain high-quality results), a benchmark study is performed to investigate its efficiency against an industrial-grade Algebraic Multigrid (AMG) solver, SAMG [ Stuben, 2010 ]. This comparative study illustrates that the C-AMS is quite efficient for compressible flow simulations in large-scale heterogeneous 3D reservoirs.

Unconventional Reservoir simulations involve several challenges not only arising from geological heterogeneities, but also from strong nonlinear physical coupling terms. All exiting upscaling and multiscale methods rely on a classical sequential formulation to treat the coupling between the nonlinear flow-transport equations. Unfortunately, the sequential strategies become severely inefficient when the flow and transport equations are strongly coupled. Examples of these cases include compositional displacements, and processes with strong capillarity effects. To extend the applicability of the multiscale methods for these challenging cases, in this paper, we propose a Constrained Pressure Residual Multiscale (CPR-MS) method. In the CPR-MS method, the CPR strategy is used to formulate the pressure equation, the approximate conservative solution of which is obtained by employing a few iterations of the iterative multiscale procedure. Several local- (ILU(k), BILU(k), etc.) and global-stage (Multiscale Finite Volume, MSFV, and Multiscale Finite Element, MSFE) solvers with different localization conditions (Linear BC, Reduced Problem BC, etc.) are employed in order to find an optimum strategy for the highly nonlinear compositional displacements. Numerical results for a wide range of test cases are presented, discussed and future studies are outlined.

In many reservoir-modeling applications, geological uncertainty is treated by considering multiple realizations generated from a specified geological scenario. In reality, however, the geological scenario itself is uncertain, and the use of qualitative criteria for its specification may lead to inaccuracy in flow predictions. In this work, we introduce a systematic procedure for the determination of the most likely a posteriori geological scenario. As is common in geomodeling applications, the scenario is defined in terms of a training image (TI). We introduce continuous parameterizations for uncertain TI attributes such as channel thickness and orientation, and then determine optimal values for these and other key model quantities. Optimality is defined here in terms of the level of agreement between observed data and flow results for realizations generated from a given geological scenario. The optimum scenario is determined using Particle Swarm Optimization. In the second step of the methodology, a set of specific realizations, which provide closer agreement with observed data, is identified using a rejection-sampling method. The workflow is applied to a synthetic channelized system, and the procedure is repeated using several different ‘true’ reservoir realizations to gauge its performance with data from realizations that are more or less representative of the true scenario. Values for TI parameters found by our procedure are shown to be in reasonable agreement with those of the true scenario in all cases considered. In addition, following the identification of specific realizations using rejection sampling, predicted flow results are shown to be of similar quality to those from the true scenario.

Smart wells are advanced operation facilities used in modern fields. Typically, a smart well is equipped with downhole sensors that collect and transmit, for instance, pressure and temperature data in order to monitor well and reservoir conditions in the field. For economical reasons, however, the number of downhole sensors is limited. Therefore, they may not be able to provide complete information about the properties of the fluids, e.g., the flow rates, in places other than the locations of the sensors. In order to evaluate fluid properties in the well, one needs to estimate them based on the collected data from the sensors. Such an exercise is often called “soft sensing” or “soft metering” (see, for example, Lorentzen et al., 2010 ).

In this work the authors study the multiphase flow soft-sensing problem based on the framework used in Lorentzen et al. (2013). There are three functional modules in this framework, namely, a transient well flow model that describes the response of certain physical variables in a well, for instance, temperature and pressure, to the flow rates entering and leaving the well zones; a Markov jump process that is designed to capture the potential abrupt changes in the flow rates; and an estimation method that is adopted to estimate the flow rates in the Markov jump process, based on the measurements from downhole sensors.

In Lorentzen et al. (2013), the variances of the flow rates in the Markov jump process are chosen manually. To fill this gap, in the current work two approaches are proposed in order to optimize the variance estimation. Through a numerical example, we show that, when the estimation framework is used in conjunction with these two proposed variance-estimation approaches, it can achieve good performance in terms of matching both the measurements of the physical sensors and the true underlying flow rates.

Geologic facies distributions are commonly represented in geomodels by categorical variables that are intrinsically non-Gaussian and thus difficult to calibrate in ensemble Kalman filter-like algorithms.

For certain types of stochastic models such as the truncated plurigaussian, it is possible to directly update model variables in such a way that the resulting realizations appear to be samples from the posterior. For other types of models, this is not possible. One common approach has been to invert flow data using the ensemble Kalman filter (EnKF) to obtain “probability maps” which are then used to condition facies realizations. Data matches obtained in this method are generally poor, however, because the probability map neglects important joint probabilities of model parameters imposed by flow data. In this paper, we propose a data assimilation method with a post-processing step that resembles the post-smoothed maximum-likelihood (ML) reconstruction method described in Nuyts et al. (2005 ). Disregarding the categorical feature of the facies model, reservoir properties are first updated using an EnKF-like assimilation method to honor flow data. In the post-processing step a penalty term forcing model variables to take discrete values is jointly minimized with the distance to the posterior realizations to solve for facies models that match data. The distance to posterior realizations is quantified using the ensemble representation of the posterior covariance, which represents the joint probability of model parameters. The matrix inversion lemma is used in solving the minimization problem to avoid inversion of the covariance.

The ability of the ensemble to accurately represent information in data is demonstrated on two linear examples and a nonlinear reservoir flow example. Comparison is made with approaches that use only the probability map to represent the assimilated data. The results show better data matches obtained with the proposed method and reflect the importance of the information captured by the updated ensemble from the data with respect to the joint probabilities of model variables.

With the advent of Ensemble methods for Assisted History Matching of reservoir models, it has become increasingly common to use initial ensembles of reservoir models consisting of a fixed (certain) grid geometry and variable (uncertain) 3D facies distributions and 3D petro-physical property distributions conditional upon facies. In such process, 3D facies realizations are often derived by modifying facies boundary locations around certain known 3D features (often from seismic interpretation).

In the context of a fixed grid, distributions of petro-physical parameters are biased because all the models in the initial ensemble are such that facies boundaries are located at grid boundary locations; in reality, the facies boundary can be located anywhere, and the petro-physical properties of cells containing facies boundaries differ from strictly facies-dependent statistics. The bias results, typically, in excessive variance of petro-physical parameters in initial ensembles in cells which possibly belong to different facies. Such variance bias has detrimental consequences on iterative inversion processes, incrementally introducing a cumulative error on the match term.

A 3D Facies Boundary Filter has been developed to statistically correct this bias over ensembles of initial models. The filter relies on the construction of a refined grid and simple averaging of petro-physical properties and is based on a hypothesis of piecewise linearity of facies boundaries between cell centres. In any ensemble approach, a filtered ensemble would be used to compute unbiased updates and such updates applied upon un-filtered ensembles. The filter helps obtain better matched models (closer to initial ensemble and/or, in synthetic cases, the truth) and ultimately determines better performance forecast. The bias will be illustrated through a simple 1D example, a 3D synthetic dataset and a real 3D field model. The principle of the filter will be detailed. The impact of the filter on the petro-physical statistics will be illustrated as well as its impact within the context of EnKF Assisted History Match.