Abstract 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.

Abstract 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.

Abstract A tightly integrated high performance parallel simulator for the management and optimization of multiple reservoirs, wells, and surface facilities is presented with subsea and giant field example applications cases. These cases have helped to develop and improve subsequent simulator releases. Specialized production management algorithms are implemented by flexible user-developed procedures written in a specialized object-based Fortran-like language. Tight coupling enables accurate prediction of all phase production rates and multiple gas and water injection rates. The new reservoir simulation tool was applied to the multireservoir Greater Plutonio (GtP) field in offshore Angola and to the giant Prudhoe Bay Unit (PBU) field in Alaska. These cases have provided the greatest challenges to the simulator in the BP system and have driven development. In both cases, the simulations were developed coupling the reservoir to the surface facilities, and user-developed procedures were included to manage and optimize production. Numerous algorithmic improvements were made to the base code to accommodate additional required functionality for the user management procedures and to achieve sufficient performance for the extensive PBU network model. For GtP, four fields were included in the full network model with multiple wells, risers, manifolds, and seabed pipeline connections. Procedures were used to model the gas injection, export, riser gas lift optimization, and facility and production constraints. Another procedure replicated the yearly production engineering optimization by shutting in wells by water-cut class for week-long tests each January. From these tests, the procedure selected the best combination of wells to accelerate oil production for the remainder of the year. Results have shown good agreement with actual production. For the PBU, a development campaign is underway to enable predictive simulation, including algorithms for predictive well management (PWM), seam tuning, and automatic drilling. The PWM algorithm for the optimal allocation of wells to low- and high-pressure systems in each of the six PBU flow stations will be implemented by user procedures. New functionality is also being developed to improve parallel operation of the network and to enable intelligent selective use of the network over time. The effort is expected to lead to substantially improved simulator performance in terms of flexibility and practicality. Many of the new features have benefited other assets. The updated simulator represents a step change in reservoir engineering management capability for the operator, resulting from increased accuracy, performance, and user flexibility. Tight coupling with scalable high performance enables accurate and stable prediction and operation. The inclusion of an extensive user procedure facility for customizable production management algorithms enables both mirroring of existing management practices and the testing of a variety of proposed alternatives. User procedures can be reused and shared among assets.

Abstract In a generalized reservoir simulator for compositional and black oil models, the primary governing equations are mass conservation equations and volume balance equations with pressure and component masses as the primary unknowns. The phase equilibrium is honored in the calculation of partial volume derivatives. The Newton-Raphson method is widely chosen to solve this non-linear system: the equations are first linearized then solved using direct or iterative linear solvers. Direct solvers are not practical for large systems; therefore, iterative solvers are used. Iterative solvers are not exact and are usually converged to a specified criterion. The stricter the criterion, the higher the computational effort required by the linear solver. On the other hand, a loose convergence criterion can lead to more Newton iterations or false physical solutions. In addition to the convergence tolerances for linear solvers, certain criteria must be applied to the non-linear Newton iterations. One of the main convergence criteria is that the volume error at any gridblock is less than a small fraction of its pore volume. The linear solver tolerance should be somehow related to this non-linear tolerance, or the linear solutions should be adjusted to help the convergence of the non-linear iterations. Because of the limited accuracy of iterative linear solvers, a two-pass procedure is proposed to improve the linear solutions by redistributing linearized mass balance errors. In the first step, all the gridblocks are divided into different groups according to whether the mass change is positive or negative and by the magnitude of contribution to local fluid volume change. Then a constrained optimization problem is defined and solved to adjust the linear solutions, satisfying global mass conservation. With these modified solutions, the linearized volume errors for all the gridblocks are calculated. In the second step, a criterion similar to the aforementioned volume convergence criterion for Newton iteration is imposed on the linearized volume error for each gridblock. With reasonable linear solutions, most of the gridblocks should satisfy this criterion. For gridblocks that do not satisfy this condition, another constrained optimization problem is defined to ensure local volume balance with local mass solutions further adjusted. With this two-pass procedure, global mass conservation and local volume balance are achieved for the implicit formulation of the governing equations, with little or no degradation of the Newton convergence, and usually without needing to tighten the linear solver convergence tolerances.

Abstract Equation-of-state (EOS) fluid characterization is used to model the behavior of hydrocarbon reservoirs when variation in the fluid composition has a significant influence on the recovery of hydrocarbons. Examples are miscible gas-injection processes, where mass transfer between the injected gas and in-situ hydrocarbons can result in the injected fluid developing into a fluid that is miscible with the in-situ hydrocarbons, or gas-condensate fluids, where the liquid yield varies with pressure and composition as the reservoir depletes. For many reservoirs, more than one recovery mechanism is employed. For example, insufficient supply of miscible injectant might result in part of a reservoir employing miscible injection, while part is under waterflood. In such cases, accurate prediction of recovery for the miscible injection area might require an EOS characterization with a large number of components, while the phase behavior in the waterflood area could be modeled with sufficient accuracy with many fewer components. Phase-behavior calculations become much more computationally expensive as the number of components increases, but current commercial reservoir simulators must use the same number of components everywhere in the reservoir model. Thus, the requirement to use a large number of components to model the recovery process in part of the reservoir results in a large computational penalty in regions of the reservoir where a smaller number of components would suffice. In this paper, we describe a method whereby the fluid can be locally lumped for phase-behavior calculations so that regions that are able to maintain sufficient compositional accuracy with fewer components can use less-expensive EOS calculations. Different lumpings can be used at different times in the life of the reservoir, as compositional effects become more or less important in different regions in the reservoir. Example problems that demonstrate the flexibility, efficiency, and consistency of the method are explored.

Abstract When multiple reservoirs are produced through a common facility network, the capability to integrate the modeling of surface and subsurface can be critical to field development and optimization. The shared facility network imposes constraints that the combined production cannot exceed, determines the pressure drop in the flow lines, and the composition and volume of the sales and reinjection streams. Pressure drop in flow lines is particularly important in deepwater field development, where flow lines are long, and production from multiple reservoirs can flow through the same riser. The most robust method for solving the combined surface-subsurface system is to fully couple and simultaneously solve the reservoir and facility equations. However, if the reservoirs fluids are represented with compositional equation of state models, and different pseudo-components are used in some or all of the reservoirs, then the reservoir fluid must be delumped into a common set of pseudo-components in the network. In this paper, we describe a method to consistently and efficiently model such a system with a fully coupled surface-subsurface simulator. At every point in the network, where fluid from only a single reservoir is present, the phase behavior calculations can use the fluid characterization for that reservoir and exactly reproduce the result that would have been obtained if the fluid had not been delumped into the network pseudo-components. However, at any point in the network, the characterization of the common network fluid can be used. In addition, if more compositional accuracy is required, the network fluid can be further delumped into more pseudo-components, or if accuracy is not critical, the fluid can be lumped into fewer pseudo-components for more computational efficiency. Examples demonstrating the flexibility, efficiency, and consistency of the method are provided.

The need for flexible and efficient multi-porosity, reservoir-simulation capabilities has never been greater. Much of the world’s oil reserves are contained in highly variable fractured reservoirs. Moreover, there is now unprecedented interest in the simulation of unconventional gas reservoirs, where up to four porosity types may be required. This paper discusses the design of an N-porosity, full-featured reservoir simulator capable of handling these diverse scenarios. Essential to the design is the programming data-structure paradigm of the SubGrid. SubGrids represent collections of simulation nodes within spatial regions of interest and encapsulate all data and variables required. Any number of associated SubGrids may coexist, representing N-porosities. Inter- and intra-SubGrid connections may be arbitrarily specified. Because the linear solver operates on the associated SubGrids as distinct entities, the design has tremendous flexibility. For example, regions with simple connectivity, such as in dual-porosity, single-permeability models, can have simple pre-solver elimination performed while other regions with uniform intra-porosity connectivity require a merged SubGrid solution. The scheme efficiently treats regions with different numbers of porosities, levels of interconnectivity, and levels of implicitness. Applications may be applied across multiple reservoirs simultaneously. Treatment of missing "fracture" zones has always been an important feature in the dual-porosity, single-permeability simulation. The traditional approach is to treat the matrix within missing fracture zones effectively as though it were an interconnected portion of the fracture porosity. This preserves connectivity within the fracture-free zones and the more computationally efficient single-permeability solver pre-elimination step. A second missing "fracture" scheme is introduced, which differs in the way that the connections between the fracture and matrix blocks are computed at the boundary of fracture-free zones. Flexible transfer-function application and porosity dependency during equilibrium initialization are discussed.

Abstract In a fully integrated reservoir and surface facility simulator, an overlapping multiplicative Schwarz method was chosen to solve the coupled system, using perforated grid blocks as the overlapping layer. But in many cases it was found that this overlapping approach did not improve global linear solver performance, and the matrix for the extended surface network, which included the perforated grid-blocks, was much larger and denser and lost its original tree structure, which taxed the linear solver for the network. Furthermore, the pressure solver for the reservoir domain was totally decoupled from the network domain. The loose coupling in the pressure solution led to the disappointing performance. In most cases, the accuracy of the pressure solution across the reservoir and surface network directly determines the performance of the global linear solver, so it is crucial to define an appropriate global pressure matrix to represent the flow exchange between these two domains. Taking advantage of a new formulation for a generalized network model of wells and facilities in which node-based variables, pressure and component compositions, are chosen as the primary variables, instead of mixed node and connection- based variables, algebraic methods are designed to reduce the full system matrix, involving pressure and component masses for the reservoir domain and pressure and component compositions for the network domain, to a pressure-only matrix. This global pressure matrix works as the first stage preconditioning matrix in a two-stage solution method. For the reservoir domain, a widely-used IMPES-like reduction method was implemented. This paper focuses on methods to construct the pressure matrix for the network domain and coupling matrices between the domains. Performance comparisons between different approaches are presented.

Summary There is increasing interest in modeling networks of wells, including subsurface components of complex wells and surface facilities. Such modeling requires setting constraints at various points in the network. Typical constraints are maximum phase flow rates and minimum flowing pressures. A major difficulty in network calculations is determining which of these constraints is active. This paper presents a method that uses slack variables in determining active constraints. The linearized equations of interest generally come in pairs, with each pair consisting of a base equation and a constraint equation. The base equation is the equation that normally applies. The constraint equation replaces it if the constraint is active. Normally, only one of these two equations can be satisfied. The slack variable provides a way to ensure that both are satisfied, regardless of which is active. If the constraint is inactive, the slack variable is added to the constraint equation and accounts for the slack, which by definition is the amount by which the inactive equation is not satisfied. On the other hand, if the constraint is active, the slack variable is instead added to the base equation, and the constraint equation as originally written is satisfied. To obtain this behavior, we define a parameter w and add w times the slack variable to the base equation and (1 − w) times the slack variable to the constraint equation. Thus, if w = 1, the slack variable is added to the base equation, and the constraint is active. On the other hand, if w = 0, the slack variable is added to the constraint equation, and the base equation is active. The slack is always in the inactive equation. There is a w associated with each slack variable. Determining the parameter w is an iterative process. The efficiency of the process is improved by manipulating the network matrix such that we can create a Schur complement that has the slack variables as its unknowns and contains the only references to the w's To determine the slack variables, we need only to work with this matrix, which typically is much smaller than the network matrix. The resulting method is implemented within a general purpose reservoir simulator. Testing of the method in more than 700 cases has shown it to be much more robust than an earlier heuristic procedure.

Summary Artificial lift by means of gas injection into production wells or risers is frequently used to increase hydrocarbon production, especially when reservoir pressure declines. We propose an efficient optimization scheme that finds the optimal distribution of the available gas lift gas to maximize an objective function subject to surface-pipeline-network rate and pressure constraints. This procedure is a nonlinearly constrained optimization problem solved by the generalized reduced-gradient (GRG) method. The values of objective function, constraint functions, and derivatives needed for optimization can be evaluated through two methods. The first method repeatedly solves the full-network equations using Newton iteration, which takes into account the flow interactions among wells; however, this method can be computationally expensive. The second and more efficient method is a new approach proposed in this paper. It constructs a set of proxy functions that approximates the objective function and constraints as functions of gas lift rates. The proxy functions are obtained by solving part of the network that consists of a gas lifted well or riser, assuming a stable pressure at the terminal node where the partial network is decoupled from the rest of the network, and are used to inexpensively evaluate the objective function, constraints, and necessary derivatives for the optimizer. A procedure to predict the proxy functions on the basis of previous values can be used to reduce the number of partial-network solves, and the partial-network solution has been parallelized for faster simulation. These two methods can be applied at different timesteps during the course of the simulation. The proposed methods are implemented within a general-purpose black-oil and compositional reservoir simulator and have been applied to real-field cases.