bluuebunny/crossref_metadata_2025 huggingface dataset:
A Discontinuous Galerkin Finite Element Framework for the Direct Numerical Simulation of Flow on High-Resolution Pore-Scale Images
Authors: Frank Florian , Liu Chen , Alpak Faruk O. , Araya-Polo Mauricio , Riviere Beatrice
Published: February 2017

Advances in pore-scale imaging, increasing availability of computational resources, and developments in numerical algorithms have started rendering direct pore-scale numerical simulations of multiphase flow on pore structures feasible. In this paper, we describe a two-phase flow simulator that solves mass and momentum balance equations valid at the pore scale, i.e. at scales where the Darcy velocity homogenization starts to break down. The simulator is one of the key components of a molecule-to-reservoir truly multiscale modeling workflow.

A Helmholtz free-energy driven, thermodynamically based diffuse-interface method is used for the effective simulation of a large number of advecting interfaces, while honoring the interfacial tension. The advective Cahn–Hilliard (mass balance) and Navier–Stokes (momentum balance) equations are coupled to each other within the phase-field framework. Wettability on rock-fluid interfaces is accounted for via an energy-penalty based wetting (contact-angle) boundary condition. Individual balance equations are discretized by use of a flexible discontinuous Galerkin (DG) method. The discretization of the mass balance equation is semi-implicit in time; momentum balance equation is discretized with a fully-implicit scheme, while both equations are coupled via an iterative operator splitting approach.

We discuss the mathematical model, DG discretization, and briefly introduce nonlinear and linear solution strategies. Numerical validation tests show optimal convergence rates for the DG discretization indicating the correctness of the numerical scheme. Physical validation tests demonstrate the consistency of the mass distribution and velocity fields simulated within our framework. Finally, two-phase flow simulations on two real pore-scale images demonstrate the utility of the pore-scale simulator. The direct pore-scale numerical simulation method overcomes the limitations of pore network models by rigorously taking into account the flow physics and by directly acting on pore-scale images of rocks without requiring a network abstraction step or remeshing. The proposed method is accurate, numerically robust, and exhibits the potential for tackling realistic problems.
Simulation of Matrix-Fracture Interaction in Low-Permeability Fractured Unconventional Reservoirs
Authors: Ding D. Y. , Farah N. , Bourbiaux B. , Wu Y-S. , Mestiri I.
Published: February 2017

Unconventional reservoirs, such as shale-gas or tight-oil reservoirs, are generally highly fractured, including hydraulic fractures, stimulated and non-stimulated natural fractures of various sizes, embedded in low-permeability formations. One of the main production mechanisms in unconventional reservoirs is the flow exchange between matrix and fracture media. However, due to extremely-low matrix permeability, the matrix/fracture exchange is very slow and the transient flow may last several to tens of years, i.e. almost the entire production life. The commonly-used dual-porosity modelling approach involves a computation of pseudo-steady-state matrix-fracture transfers with homogenized fluid and flow properties within the matrix medium. This kind of model clearly fails to handle the long-lasting matrix/fracture interaction in very-low permeability reservoirs, especially for multi-pahse flow with phase change problems. Moreover, a dual-porosity model is not adapted for the simulation of matrix/fracture exchange when fractures are described by a DFN (Discrete Fracture Network). This paper presents an EDFM (Embedded Discrete Fracture Model) based on the MINC (Multiple INteracting Continua) proximity function to overcome this insufficiency of the conventional dual-porosity model.
Multirate Coupling for Flow and Geomechanics Applied to Hydraulic Fracturing Using an Adaptive Phase-Field Technique
Authors: Almani Tameem , Lee Sanghyun , Wheeler Mary F. , Wick Thomas
Published: February 2017

We present and analyze a multirate fixed stress split iterative coupling scheme for coupling flow and geomechanics in a poroelastic medium involving fracture propagation modeled with a phase field approach. The novelty of this work lies in the efficient integration of the fixed-stress split coupling scheme with phase-field fracture propagation models. The multirate coupling algorithm utilizes different time-scales of the flow and mechanics problems, by allowing for multiple finer time steps for flow within one coarse mechanics time step. When applied to production scenarios, the multirate scheme results in massive reductions in the number of mechanics linear iterations, without jeopardizing the accuracy of the obtained results. A number of numerical simulations substantiate our algorithmic developments. These tests include prototype computations, multiple propagating fractures, and fractures initialized by a microseismic probability map.
Field-Scale Assisted History Matching Using A Systematic, Massively Parallel Ensemble Kalman Smoother Procedure
Authors: Lin Binghuai , Crumpton Paul I , Dogru Ali H.
Published: February 2017

This work presents a systematic and rigorous approach of reservoir decomposition combined with the ensemble Kalman smoother to overcome the complexity and computational burden associated with history-matching field-scale reservoirs in the Middle East. The results demonstrate the performance and feasibility of this integrated method.

The ability of the ensemble smoother to quantify the uncertainty of complex nonlinear models and handle a high dimensional parameter space has made it an extremely useful algorithm. However, for large fields with their associated amount of observational data the expense of this method can be prohibitive. To overcome this issue, the reservoir is decomposed into overlapping domains based on the prior knowledge of the geological and fluid flow features of the reservoir, and an improved iterative ensemble Kalman smoother is applied to each domain. For each iteration the ensemble updating is regularized to control the parameter increments and avoid filter divergence. Localization is employed to minimize spurious long-distance correlations and reduce the boundary effects between the domains. Each sub-domain exploits a parallel cluster, both to run each simulation within the ensemble independently, and then to run each reservoir simulation in parallel. The Kalman update itself cannot be implemented in a naive way because of the memory consideration of the large matrices involved. Instead, the HDF5 file format is utilized for efficient out-of-core memory reading/writing of the large matrices and facilitating the matrix operations.

The method is applied to a 60-million-cell Middle East field model. The systematic approach here decomposes the original reservoir into six smaller pieces. The ensemble smoother, with an ensemble size of 50, is then applied to the individual reservoir domains to generate history-matched models; this exploits a 3000-node cluster with 24 cores on each node. Each of the simulation runs takes 2-3 hours on 500 cores; the overall history matching process can be performed in about one week. The results of the models are then combined together to reconstruct the original reservoir. Compared to history-matching results from months of effort of the engineers, the ensemble smoother takes much less manpower to achieve better match to the data.

This parallel iterative ensemble smoother has demonstrated the feasibility and performance of history-matching large-scale complex real field reservoirs. By combining domain decomposition, the systematic history-matching workflow can also be applied to large Middle-East reservoirs, which generally cannot be achieved using other methods.
Improving the Computational Efficiency of Approximate Gradients Using a Multiscale Reservoir Simulation Framework
Authors: Moraes R. , Fonseca R. M. , Helici M. , Heemink A. W. , Jansen J. D.
Published: February 2017

We present an efficient workflow that combines multiscale (MS) forward simulation and stochastic gradient computation - MS-StoSAG - for the optimization of well controls applied to waterflooding under geological uncertainty. The Iterative Multiscale Finite Volume (i-MSFV), a mass conservative reservoir simulation strategy, is employed as the forward simulation strategy. The i-MSFV method has the ability to accurately capture fine scale heterogeneities, and thus the fine scale physics of the problem. This allows for accurate and efficient simulations with fine-scale error estimates in a more computationally efficient coarse-scale simulation grid. The Stochastic Simplex Approximate Gradient (StoSAG) method, a stochastic gradient technique, inspired from Ensemble Optimization (EnOpt), is employed to compute the gradient of the objective function. Stochastic methods have recently been shown to be very efficient in terms of gradient quality and computational cost especially for robust optimization problems. They rely on the reservoir simulator as a black-box and perform a large number of forward simulations to approximate a gradient. In our numerical experiments, we investigate the impact of MS parameters such as coarsening ratio and heterogeneity contrast on the quality of the approximate MS-StoSAG gradient. Our experiments illustrate that i-MSFV simulations provide accurate forward simulation responses for the gradient computation, with the advantage of speeding up the workflow due to faster simulations. Speedups up to a factor of five on the forward simulation were achieved for the examples considered in the paper. The combination of speed and accuracy of MS forward simulation with the flexibility of the StoSAG technique (several different parameters/responses can be easily considered) allows for a flexible and efficient optimization strategy suitable for large-scale problems.
Mechanistic Simulation Workflow in Shale Gas Reservoirs
Authors: Yan Bicheng , Mi Lidong , Wang Yuhe , Tang Hewei , An Cheng , Killough John E.
Published: February 2017

Shale gas reservoir is comprised of highly heterogeneous porosity systems including hydraulic/secondary fractures, inorganic and organic matrix. Multiple non-Darcy flow mechanisms in the shale matrix further bring challenges for modeling. In this paper, we developed a framework combining a multi-physics compositional simulator with Multi-Porosity Modeling preprocessor for gas storage and transport in shale. A Triple-Porosity Model is used to characterize the three porosity systems in shale gas reservoirs. In the fracture porosity the heterogeneous impact of secondary fractures distribution on matrix-to-fracture fluid transfer is revealed by shape factor distribution. They are upscaled with superior accuracy from a detailed Discrete Fracture Network Model (DFN) sector model, where orthogonal hydraulic fractures are explicitly discretized. With the occurrence of nano-pores in shale matrix, the interaction between pore-wall and gas molecules is considered via Knudsen diffusion and gas slippage. Gas adsorption on the pore-wall of organic matrix is modeled by extended Langmuir isotherm. The inter-porosity and intra-porosity connectivities in the Triple-Porosity Model are flexibly controlled by arbitrary connections. Our results show that gas production in the Triple-Porosity Model with shape factor upscaled from DFN exhibits different production performance from models with uniform shape factor distribution. The deviations are caused by the dominance of different regions at different production periods. Moreover, different combinations of flow and storage mechanisms are investigated. We show that Langmuir desorption maintains reservoir pressure, but gas slippage and Knudsen diffusion accelerate the pressure drop. Both mechanisms contribute to improve the gas production and the consideration of them simultaneously improve gas production most.
Evaluation of Reservoir and Production System Integration in Production Strategy Selection
Authors: von Hohendorff Filho João Carlos , José Schiozer Denis
Published: February 2017

In offshore petroleum field studies, the integration between reservoir and production system simulators may improve production forecasts. Therefore, it is important to check if this integration is necessary and a methodology to do it adequately. In this work, we analyzed the influence of this integration in the development of a petroleum field focusing on the effects on the production strategy parameters for a benchmark model based on an offshore field in Brazil.

The methodology is based on a 12-step procedure for closed loop reservoir development and management. The production facilities are integrated with the reservoir in Step 11 and we propose to re-optimize the production strategy using information about pressure gradients inside the satellite well pipes for production forecast, while evaluating net present value as the objective-function. We adapted an assisted optimization workflow to include the optimization of new variables such as, pipe diameters of production system, platform position and artificial lift application, and compared this with the production strategy obtained from the same benchmark in a nonintegrated approach.

Comparing nonintegrated and integrated production strategies, we observed changes that indicate the need to integrate reservoir and production systems in the decision making process. (1) Platform position was more affected by economic indicators than production loss due to pressure drop in flowlines. (2) Economic return and field oil recovery were significantly affected by the diameters of the production pipes. Artificial lift application with high gas injection rates enhanced well performance. (3) A re-optimization of the number of wells had little effect on the cumulative oil production but larger influence on net present value, indicating potentially unnecessary wells in a nonintegrated project. (4) The integration demanded further optimization of well placement, significantly impacting both net present value and the recovery factor.

The integrating production facilities with the reservoir, generated a significant increase in the net present value, maintained the same oil recovery factor and required less investment.
Multiscale Gradient Computation for Multiphase Flow in Porous Media
Authors: Moraes R. , Rodrigues J. R. , Hajibeygi H. , Jansen J. D.
Published: February 2017

A multiscale gradient computation method for multiphase flow in heterogeneous porous media is developed. The method constructs multiscale primal and dual coarse grids, imposed on the given fine-scale computational grid. Local multiscale basis functions are computed on (dual-) coarse blocks, constructing an accurate map (prolongation operator) between coarse- and fine-scale systems. While the expensive operations involved in computing the gradients are performed at the coarse scale, sensitivities with respect to uncertain parameters (e.g., grid block permeabilities) are expressed in the fine scale via the partial derivatives of the prolongation operator. Hence, the method allows for updating of the geological model, rather than the dynamic model only, avoiding upscaling and the inevitable loss of information. The formulation and implementation are based on automatic differentiation (AD), allowing for convenient extensions to complex physics. An IMPES coupling strategy for flow and transport is followed, in the forward simulation. The flow equation is computed using a multiscale finite volume (MSFV) formulation and the transport equation is computed at the fine scale, after reconstruction of mass conservative velocity field. To assess the performance of the method, a synthetic multiphase flow test case is considered. The multiscale gradients are compared against those obtained from a fine-scale reference strategy. Apart from its computational efficiency, the benefits of the method include flexibility to accommodate variables expressed at different scales, specially in multiscale data assimilation and reservoir management studies.
A Framework for Mechanistic Modeling of Alkali-Surfactant-Polymer Process in an Equation-of-State Compositional Simulator
Authors: Nghiem L. , Skoreyko F. , Gorucu S. E. , Dang C. , Shrivastava V.
Published: February 2017

Alkaline/Surfactant/Polymer (ASP) is an important chemical EOR process that involves the generation of in-situ soap through the reaction of the acid component in the oil with the alkali, in conjunction with intra-aqueous reactions, mineral dissolution and precipitation, micro-emulsion behavior and salinity gradient. The mechanistic simulation of ASP is very complex and has been carried out with specialized chemical flood simulators. Current trends in EOR processes show an interest in hybrid methods where chemical flooding is combined with other EOR methods such as low salinity waterflood, foam and gas/CO2 injection. Thus, it is beneficial to develop such capabilities in an Equation-of-State (EOS) compositional simulator for screening and combining different EOR processes.

This paper presents a framework for mechanistic modeling of the ASP process within an EOS compositional simulator. A new approach for modeling the Winsor Type I, II and III micro-emulsion phase behavior is introduced based on laboratory solubility data. In the Type III system, the emulsion is distributed judiciously between the oil and water phase without the introduction of a third liquid phase. The optimal salinity variation with the soap/(soap + synthetic surfactant) mole fraction is modeled. This feature allows the design of salinity gradient, an essential requirement for a successful ASP flood. The above physics are coupled with comprehensive geochemistry calculations (intra-aqueous reactions and mineral precipitation/dissolution reactions) and three-phase oil/gas/water flash calculations with an equation of state and Henry's law. The whole system of associated equations is solved simultaneously with the flow and energy equations using Newton's method, making the simulator one of the most robust and comprehensive simulators for EOR methods.

The simulator is validated with ASP core flood experiments. The optimization of the ASP process for a typical field is illustrated and discussed with regard to alkali, synthetic surfactant and polymer injection with decreasing salinity.
On the Convergence of System-AMG in Reservoir Simulation
Authors: Gries Sebastian
Published: February 2017

System-AMG provides a flexible framework for linear systems in simulation applications that involve various different types of physical unknowns. Reservoir simulation applications, with their driving elliptic pressure unknown, are principally well-suited to exploit System-AMG as a robust and efficient solver method.

However, in order to be efficient and robust, the coarse grid correction process of AMG on the one hand needs to be possible, i.e., the matrix needs to fulfill certain requirements. It has been demonstrated earlier how to ensure this by the dynamic rowsumming (DRS) method.

On the other hand, the coarse grid correction must be physically meaningful in order to speed up the overall convergence. It has been common practice in CPR-type applications to use an approximate pressure-saturation decoupling to fulfill this requirement. This, however, can have drastic impacts on AMG’s applicability and, thus, is not performed by the DRS-method.

In this work, we are going to see that the pressure-saturation decoupling indeed is not necessary for ensuring an efficient interplay between the coarse grid correction process and the fine-level problem. We will find that a comparable influence of the pressure on the different involved PDEs is much more crucial.

As an extreme case w.r.t. the outlined requirement, we will discuss linear systems from compositional simulations under the volume balance formulation. In these systems, the pressure typically is associated with a volume-balance, rather than a diffusion process. The corresponding coarse grid correction does, hence, not provide any benefit regarding the overall convergence: the other PDEs involve pressure-based diffusive parts that have a drastically different structure than the volume-balance has. We will see how System-AMG can still be used in such cases.
A Parallel Solution for Large Surface Networks in a Fully Integrated Reservoir Simulator
Authors: Fleming Graham , Wang Qinghua
Published: February 2017

Integrated modeling of surface and subsurface is necessary to ensure physical accuracy in reservoir simulation with complex wells and surface facilities. In this fully implicit coupled approach, two physical domains, reservoir and surface network, are assembled into a single global system and solved simultaneously. However, flows in these two domains behave very differently. The surface network, comprised of wells and surface facilities, usually requires many more Newton iterations than the reservoir domain; therefore, the global Newton iteration begins with a standalone network solution as a preconditioner, followed by a global linear system solution. For large surface networks, the computational effort necessary to solve the surface network equations can be significant. Parallelization of the reservoir domain is a well-studied topic, but little consideration has been given to parallel solutions for integrated models with large facilities. In this paper, an approach to parallelization of the surface network is presented.

The network solution is involved in two stages, the linear global solution and the nonlinear standalone preconditioning. In the global solution stage, a multiplicative Schwarz method is chosen to solve the global linear system, and a parallel sparse direct solver is applied to solve the surface network matrix. Before the global linear solution, the coefficient matrix for the surface network is reduced to a Schur complement with connection-based variables eliminated, including total flow rate on network connections and component flow rates on perforation connections. This elimination is also performed in parallel. The parallelization of the nonlinear solution of the surface network in the preconditioning stage is more challenging. Similar to the reservoir domain, a domain decomposition method with a bottom-up approach is used to partition the surface network to substructures. Fortunately, the surface network often has a sparse tree structure. Most of the network computations, constraint setting, targeting, pressure, volume, temperature (PVT), and hydraulic calculations can be parallelized readily with a message passing interface (MPI); however, load imbalance can be a problem because the network structure changes over time, and the cost of the nonlinear solutions involved in each substructure can vary dramatically.

This paper demonstrates the performance of the network parallelization for a compositional model with a large surface network. The scalability of the standalone network solution, the linear global solution, and the full system are presented and compared to a serial approach for the network.

To achieve a reasonable scalability in reservoir models with large surface facilities, new algorithms and data structures are designed to parallelize the solutions for the coupled system of reservoir and surface network.
Flow Diagnostics on Fully Unstructured Grids
Authors: Zhang Zhao , Geiger Sebastian , Rood Margaret , Jacquemyn Carl , Jackson Matthew , Hampson Gary , Carvalho Felipe Moura , Silva Clarissa Coda , Silva Julio Machado , Sousa Mario Costa
Published: February 2017

Flow-diagnostics are a common way to rank and cluster ensembles of reservoir models based on their approximate dynamic behaviour prior to commencing full-physics reservoir simulation. Traditionally, flow diagnostics are carried out on corner-point grids inherent to geocellular models.

The novel "Rapid Reservoir Modelling" (RRM) concept enables fast and intuitive prototyping and updating of reservoir models. In RRM, complex reservoir heterogeneities are modelled as discrete volumes bounded by surfaces that can be modified using simple sketching operations in real time. The resulting reservoir models are discretized using fully unstructured 3D meshes where the grid conforms to the reservoir geometry.

This paper presents a new and computationally efficient numerical scheme that enables flow diagnostic calculations on fully unstructured grids. Time-of-flight and steady-state tracer distributions are computed directly on the grid. The results of these computations allows us to estimate swept reservoir volumes, injector-producer pairs, well-allocation factors, flow capacity, storage capacity and dynamic Lorenz coefficients which all help approximate the dynamic reservoir behaviour.

We use the Control Volume Finite Element Method (CVFEM) to solve the elliptic pressure equation. A scalable matrix solver (SAMG) is used to invert the linear system. A new edge-based CVFEM is developed to solve hyperbolic transport equations for time-of-flight and tracer distributions. An optimal reordering technique is employed to deal with each control volume locally such that the hyperbolic equations can be computed in an efficient node-by-node manner. This reordering algorithm scales linearly with the number of unknowns.

The total CPU time, including grid generation and flow diagnostics, is typically below 3 seconds for grids with 50k unknowns. Such fast calculations provide, for the first time, real-time feedback on changes in the dynamic reservoir behaviour while the reservoir model is updated.
Numerical Modeling of Unstable Water Floods and Tertiary Polymer Floods into Highly Viscous Oils
Authors: Loubens R. de , Vaillant G. , Regaieg M. , Yang J. , Moncorgé A. , Fabbri C. , Darche G.
Published: February 2017

The saturation distribution after unstable water flooding into highly viscous oil may have a decisive impact on the efficiency of tertiary polymer flooding, in particular due to hysteresis effects associated with oil banking. In this work, we model water flood and tertiary polymer flood experiments performed on Bentheimer sandstone slabs with heavy oils of about 2000 cP and 7000 cP, and compare the numerical results with experimental production, pressure and X-ray data.

The unstable water floods are initially simulated in 2D with our parallel in-house research reservoir simulator using a high-resolution discretization. In agreement with existing literature, we find that Darcy-type simulations based on steady-state relative permeabilities - inferred here from a 3D quasistatic pore network model (PNM) - cannot predict the measured water flood data. Even qualitatively, the viscous fingering patterns are not reproduced. An adaptive dynamic pore network model is then applied on a 2D pore network constructed from the statistics of the 3D network. If the fingering patterns simulated with this 2D PNM are qualitatively in good agreement with the experimental data, a quantitative match cannot be obtained due to the limitations of 2D modeling. Although 3D dynamic PNM at the slab scale would currently lead to prohibitively high computational cost, it has the potential to address the deficiencies of continuum models at highly unfavorable viscosity ratio.

For the tertiary polymer floods characterized by a much more favorable mobility ratio, Darcy-type modeling is applied and history matching is conducted from the end of the water floods. We find that unless hysteresis due to oil banking is accounted for in the relative permeability model, it is not possible to reconcile the experimental data sets. This hysteresis phenomenon, associated with oil invasion into previously established water channels, explains the rapid propagation of the oil bank. For the considered experiments, a simultaneous history match of good quality is obtained with the production and pressure data, and the simulated 2D saturation maps are in reasonable agreement with X-ray data.

This paper addresses the challenges in modeling highly unstable water flooding, using both a conventional Darcy-type simulator and adaptive dynamic PNM, by confronting the simulated results with experimental data including saturation maps. It also highlights the important role of relative permeability hysteresis in the tertiary recovery of viscous oils by polymer injection.
An Improved Slack Variable Method to Determine Active Constraints in Large Surface Networks
Authors: Wang Qinghua , Fleming Graham
Published: February 2017

In a fully integrated reservoir model consisting of one or more reservoirs and a surface network of complex wells and facilities, determining the active constraints in the network, while ensuring that the network equations are nonsingular, is a challenge. This paper introduces an improved slack variable method to address this difficulty.

The high solution cost of the matrix equations associated with the slack variable method makes it inefficient for simulation of reservoirs with a large surface network. After elimination of the nonslack variables, the resultant Schur complement is often a dense matrix and loses the typical sparse tree structure of the original network Jacobian, which makes it costly to solve. This paper revisits the slack variable method and designs two approaches to address the solution cost. The first approach, the pseudo slacks, conserves the sparse pattern of the surface network by introducing several pseudo variables that are added to the slack matrix; it works best for cases with targeting or constraints for a large group of wells. The second approach, minimal slacks, uses a simple heuristic technique to remove slack variables unless they are required to ensure the equations are nonsingular. This method reduces the slack matrix dimension but can increase the number of iterations required to solve the network.

Numerical experiments show significant performance improvements compared to the previous slack variable method in reservoir models with large surface networks.
Joint Optimization of Economic Project Life and Well Controls
Authors: Shirangi Mehrdad G. , Volkov Oleg , Durlofsky Louis J.
Published: February 2017

A new methodology for the joint optimization of economic project life and time-varying well controls is introduced. The procedure enables the maximization of net present value (NPV) subject to satisfaction of a specified modified internal rate of return. Use of this framework allows an operator to avoid situations where NPV continues to increase in time, but the late-time cash flows are negligible (in terms of an appropriate financial metric) relative to the capital invested in the project. The optimization is formulated as a nested procedure in which economic project life is optimized in the outer loop, and the associated well settings (time-varying bottomhole pressures in the cases considered) are optimized in the inner loop. The inner-loop optimization is accomplished by use of an adjointgradient-based approach, while the outer-loop optimization entails an interpolation technique. We demonstrate the successful application of this framework for production optimization for two-and three-dimensional reservoir models under waterflood. The tradeoff between maximized NPV and rate of return is assessed, as is the impact of discount rate on optimal operations. We believe this to be the first production optimization formulation that explicitly incorporates both NPV and rate of return. As such, this approach may represent an alternative to existing treatments that entail the bi-objective optimization of long-and short-term NPV.
A Double Control Volume Finite Element Method with Dynamic Unstructured Mesh Optimization
Authors: Salinas P.. , Pavlidis D.. , Xie Z.. , Pain C. C. , Jackson M. D.
Published: February 2017

A novel high-order, Double Control Volume Finite Element Method in which pressure is discretized using finite volumes is presented to model multiphase flow porous media simulations. This approach shows a better behavior to calculate the pressure field than the well-known Control Volume Finite Element Method. The latter, as many other formulations, have problems to calculate the pressure field when the elements conforming the mesh, triangles (2D) or tetrahedra (3D), have very obtuse angles. This limits the advantages that can be obtained by using dynamic adaptive meshing methods, as the meshes are constrained to have good quality elements rather than just optimizing the elements to adapt to the peculiarities of the flow. This paper presents a novel formulation that shows to be more resilient with regard to elements with very obtuse angles, enabling dynamic mesh adaptivity to be more effective.
Interactions Between Viscous Fingering and Channeling for Unstable Water/Polymer Floods in Heavy Oil Reservoirs
Authors: Luo Haishan , Delshad Mojdeh , Pope Gary A. , Mohanty Kishore K.
Published: February 2017

Unstable floods and resulting viscous fingers remain a big challenge for reservoir simulation as the gridblock size is usually many orders larger than the viscous finger wavelength. This problem becomes especially pronounced with increasing applications of polymer and other chemical floods in the development of heavy oil reservoirs. Traditional reservoir simulators do not consider sub-grid viscous fingering effects and tend to overestimate the waterflood oil recovery. Using extremely fine grid models with centimeters size is unrealistic for field-scale simulations.

While some researchers disregard viscous fingering by claiming that channeling dominates at the large scale for heterogeneous reservoirs, they miss the existence of viscous fingering at the small scale, which affects the displacement efficiency. To overcome this limitation, an effective-fingering model was developed to upscale fingering effects. The model divides each gridblock into three dynamic regions: two-phase flow, single phase oil flow, and bypassed-oil regions. Model parameters represent the maximum fraction of viscous fingering and the growth rates of different regions, which are used to modify flow functions. Model parameters from history match of a set of laboratory experiments show clear power-law correlations with a dimensionless viscous finger number, a function of viscosity ratio, velocity, permeability, interfacial tension, and core cross-sectional area.

The correlation was achieved in the lab scale by considering homogeneous cores, and we extended it further to the field scale by performing high-order spatial accuracy numerical simulations at the intermediate scale using fine gridblock sizes roughly the same as that of the core. Geostatistical realizations of the permeability field were generated with various variances and correlation lengths. In a statistical way, we were able to quantify the viscous finger number valid for a gridblock at the field scale affected by various heterogeneities using the effective-fingering model. We also observed that channelized permeability distributions increase the viscous finger number drastically, showing the important role of channeling in such cases. This new model was applied to a field case with high heterogeneity undergoing water/polymer floods. We observed that the oil recovery was improved by the polymer slug because of the enhancement in both local displacement efficiency and sweep efficiency.

In summary, we developed an upscaling model that provides a fresh-new insight on how to simulate unstable water/polymer floods at the field scale, which effectively accounts for the interplay of viscous fingering and channeling.
M3E_LINSOL: A Scalable Linear Solver Package for Reservoir Simulations on Massively Parallel Computers
Authors: Ferronato M. , Paludetto Magri V. A. , Janna C. , Franceschini A.
Published: February 2017

A numerical package called M3E_LINSOL for the solution of large linear systems of equations arising from reservoir simulations is presented. This suite includes Krylov-based solvers combined with a set of Factorized Sparse Approximate Inverse (FSAI) preconditioners specifically designed for massively parallel architectures. The computational efficiency and scalability of these preconditioners is evaluated according to challenging real-world reservoir problems.
A New Multiscale Computational Model for Low Salinity Alkaline Waterflooding in Clay-Bearing Sandstones
Authors: Lima Sidarta A. , Murad Marcio A. , Domingues Ramon
Published: February 2017

We develop a new multiscale model to compute effective properties such as relative permeabilities in low salinity EOR processes in sandstones containing reactive surfaces of kaolinite clays. The approach proposed herein considers the local nanoscale description ruled by the electro-chemistry of an aqueous electrolyte solution containing Na+, Ca2+, H+ and OH- ions lying between a crude-oil droplet and clay substrate. By considering concomitant surface complexation geochemical reactions between the ionic species in the invading water and the electrically charged kaolinite surface and oilwater interface we build-up the local electric double layer problem for the electric potential based on a generalized Poisson-Boltzmann equation with nonlinear boundary conditions ruled by the geochemical reactions. By considering the local mechanical equilibrium of the electrolyte solution and solving numerically the Poisson Boltzmann problem by finite elements we compute the disjoining pressure for a wide range of brine compositions. Moreover by invoking the Frumkin/Derjaguin wetting equation, the balance between the local disjoining and capillary pressures provides highly accurate computations of the contact angle along with its dependence on pH and salinity. Such methodology allows to compute wettability alteration of the rock surface from preferential oil to water wetness. The upscaling of this result allows to compute precisely the additional dependency of relative permeability on brine composition. Computational simulations illustrate the potential of the model proposed herein in bridging the pore-scale description, ruled by the local electro-chemistry of the thin liquid film, and the constitutive law for the relative permeabilities in low salinity alkaline waterflooding.
Model Order Reduction and Control Polynomial Approximation for Well-Control Production Optimization
Authors: Sorek N. , Zalavadia H. , Gildin E.
Published: February 2017

The objective of this paper is to reduce the computational effort in reservoir flooding optimization problems by a combination of different optimization parametrization methods and model order reduction techniques. We compare three different parametrization methods that reduce the cardinality of the original infinite set of control-decision variables to a finite set. The three methods include a traditional piece-wise constant (PWC) approximation, a polynomial approximation by Chebyshev orthogonal polynomials and a piece-wise polynomial approximation by cubic Spline interpolation. The Proper Orthogonal Decomposition with Discrete Empirical Interpolation Method (POD-DEIM) accomplishes the reduced order modeling (ROM)..

We compare a gradient-free global stochastic search approach and a gradient-based local search approach. We used Particle Swarm Optimization (PSO) as a gradient-free algorithm and Interior-Point Optimization (IPOPT) with L-BFGS method as a gradient-based algorithm. First, we compare the performances of the three parametrization methods solved by each optimizer, using fine scale simulations for an increasing level of parametrization refinement. Then, in the second part of this paper, we combine the parametrization methods with the reduced modeling workflow. For a given level of parametrization refinement, we compare the performance of each parametrization method coupled with POD-DEIM, and solved by each optimizer. In this part, we introduce an online training procedure, where the first optimization iteration is used to construct the snapshot matrix.

The results demonstrate how refining the control approximation with more decision variables per well lead to better NPV values, but with a higher computational cost. The best NPV was achieved using the highest refining level with Chebyshev polynomial approximation. Both polynomial and piece-wise polynomial approximations served as better training sets for POD-DEIM leading to a more accurate and fast reduced model. With the strategy proposed, POD-DEIM showed the best optimization accuracy for Chebyshev polynomial with the gradient-free optimizer, thus permitting the use of the model reduction methodology for global-stochastic search methods. However, the gradient-based approach seems to consistently outperform the gradient-free approach in terms of NPV and number of iterations for the cases shown.
A New Four-Phase Flash Algorithm for Surfactant Flood Simulation
Authors: Han Choongyong , Chang Yih-Bor , Wolfsteiner Christian , Shi Xundan , Schrader Mark
Published: February 2017

In the recent surfactant flood simulations, there have been efforts to consider the well-known effect of oil composition on the surfactant phase behavior. In order to truly take this effect into account, hydrocarbon phase behavior needs to be coupled with the surfactant phase behavior to evaluate more accurate oil composition, especially when gaseous phase could exist in a reservoir during a simulation. Also, there is a strong need to have this coupling in a fully implicit, fully coupled (shortly, FIM) simulator in order to get all the benefits from the FIM numerical solution approach in surfactant flood simulations.

In this study, therefore, we propose a new gaseous/oleic/microemulsion/aqueous phase flash algorithm to couple the hydrocarbon phase behavior with the surfactant phase behavior in our FIM simulator using a new set of non-linear solution variables. Instead of using the phase molar fraction which is one of the conventional natural non-linear solution variables in a FIM simulator, we use various concentrations/fractions along with pressure and phase saturation as input for the four phase flash algorithm. The concentrations/fractions are the non-linear solution variables newly introduced to our FIM simulator.

With the new algorithm, we can easily consider the cation exchange on surfactant and the indirect effect of volume of surfactant adsorbed on rock in the four phase flash. We need little conversion work and don't need to locally solve equations for concentrations/fractions even if the variables are required. Also we can incorporate the original two-phase hydrocarbon flash and three-phase surfactant flash models into our four phase flash without modification. These are all the advantages over a four phase flash algorithm using the phase molar fraction as input from the non-linear solve. Moreover, as molar fraction of each component is a calculated intermediate input for the new flash algorithm, any simulator using non-linear solution variables related to component moles can easily adopt the algorithm. This actually allows us to implement the same algorithm with slight extension in one of our well and separator models which use the molar non-linear solution variables.

We have successfully tested and validated the new four phase flash algorithm using our FIM simulator with both conventional natural and mixed non-linear solution variables. Also, we show the algorithm works well with the effect of cation exchange on surfactant phase behavior. The algorithm has great extensibility to adopt any hydrocarbon and surfactant phase behavior models in our FIM simulator with the mixed non-linear solution variables if the models use their own conventional input and output variable types.
Multiscale Finite Volume Method for Discrete Fracture Modeling with Unstructured Grids
Authors: Bosma S. B. , Hajibeygi H. , Tene M. , Tchelepi H. A.
Published: February 2017

A multiscale method for Discrete Fracture Modeling (DFM) using unstructured grids is developed. The fine-scale discrete system is obtained by imposing tetrahedron (triangular for 2D domains) shaped grid cells, while lower-dimensional fractures are imposed at the grid interfaces. The DFM approach is then used to describe the transmissibility coefficients for all the interfaces, including those with the lower-dimensional fractures. On this fine-scale discrete system, a new algebraic multiscale formulation is developed, which first imposes two sets of coarseand dual-coarse grids. The former grid is essential for conservative multiscale formulations, while the second one is used for the calculation of local multiscale basis functions. The coarse-scale partitioning of the fracture and matrix domains is flexible and totally independent. Moreover, the multiscale basis functions are constructed for both the matrix and fracture domains. By construction, the basis functions are a partition of unity. For 2D and 3D test cases, the performance of the multiscale method is systematically assessed. It is shown that the method (with no multiscale iterations) provides accurate results, even for complex fractured systems. The presented multiscale method is a promising framework for real-field application of DFM models.
Robust Geochemical Simulation of Alkali/Surfactant/Polymer Flooding with an Equation of State
Authors: Khorsandi Saeid , Qiao Changhe , Johns Russell T.
Published: February 2017

Alkali-surfactant-polymer (ASP) floods can recover significant amounts of remaining oil after a water flood. Currently, commercial and academic chemical flooding simulators use Hand's model and linear interpolation rules based on a salinity scan to model micro-emulsion phase behavior. Those simple models have significant limitations when the reservoir conditions are different than experimental conditions. Recently an equation-of-state (EoS) based on HLD-NAC was developed that gives improved phase behavior predictions away from the tuned experimental data. The new EoS for SP modeling was incorporated into our in-house compositional simulator, PennSim. In this research, we expand PennSim with an updated flash calculation module that includes alkali and geochemistry. The new simulator can more accurately model ASP flooding since phase behavior changes consistently with aqueous chemistry, oil saturation, alkali concentration, salinity, pressure, and other key variables.

ASP phase behavior depends on many factors including pressure, temperature, oil EACN, and geochemistry. PennSim solves the mass conservation equations of oil acidic components, water, and aqueous solute along with phase equilibrium relations and chemical reactions using an IMPEC approach. The reactions between oleic acids, aqueous and solid species, surfactant adsorption and mineral dissolution reactions are considered. A scalable capillary desaturation curve (CDC) with permeability and porosity is used for heterogeneous field cases. An algebraic multigrid linear solver is used to speed up simulation.

The functionalities of the developed code are demonstrated with multiple injection scenarios, such as the combination of an alkali and polymer slug with and without added synthetic surfactant. Simulations are made to model 1D core floods and 2D heterogeneous reservoirs with multiple wells. The simulation results show that the main difference between the new CDC and traditional CDC curve is the size of the developed oil bank, where the traditional CDC predicts a larger oil bank in the high permeability channel.