Numerical Simulations of Polymer Flooding Process in Porous Media on Distributedmemory Parallel Computers
Abstract
Polymer flooding is a mature enhanced oil recovery technique that has been successfully applied in many field projects. By injecting polymer into a reservoir, the viscosity of water is increased, and the efficiency of water flooding is improved. As a result, more oil can be recovered. This paper presents numerical simulations of a polymer flooding process using parallel computers, where the numerical modeling of polymer retention, inaccessible pore volumes, a permeability reduction and polymer absorption are considered. Darcy’s law is employed to model the behavoir of a fluid in porous media, and the upstream finite difference (volume) method is applied to discretize the mass conservation equations. Numerical methods, including discretization schemes, linear solver methods, nonlinearization methods and parallel techniques are introduced. Numerical experiments show that, on one hand, computed results match those from the commercial simulator, SchlumbergerEclipse, which is widely applied by the petroleum industry, and, on the other hand, our simulator has excellent scalability, which is demonstrated by field applications with up to 27 million grid blocks using up to 2048 CPU cores.
keywords:
High performance computing, Polymer flooding, Water flooding, Scalability1 Introduction
The emergence of parallel computers compels parallel computation techniques into an array of application areas, including groundwater flow, contamination transport modeling, geothermal engineering, multiphase flow, carbon dioxide sequestration and nuclear waste storage Coumou2008 (). Beginning in the mid1970s, supercomputers were introduced to accelerate reservoir modeling problems through vectorization, and computations could be completed at an advanced speed. However, models and programs had to be reorganized and reworked to take advantage of leveraged computational power through vectorization. Besides, the program performance deteriorated once a CPU number went beyond a specific number (usually 4, 8 or 16) Dogru2000 (). Except on a shared memory system, parallel computations can also be carried out on distributed memory clusters. Over the past few decades, significant progress has been made in developing high performance modeling tools for distributed memory systems. However, they have not been widely applied in reservoir simulations. The reasons are complicated but can be interpreted as follows: First, the need for high performance computing for a petroleum reservoir simulation model is dominated by the physics of the underlying process. The severely nonlinear nature of a physical process challenged the development of efficient parallel schemes for reservoir simulations that result from the multiphase, multicomponent flow through heterogeneous porous media with complex phase behavior equilibrium calculations. Second, the rapid advance in computing hardware strengthened the traditional oneprocessor simulators, such as vectorization techniques in sharedmemory machines Coats1987 (). It hindered the development of high performance parallel schemes, since PCs are much more affordable and available than a cluster. The last but not the least, the earlier developed parallel algorithms depended on a machine structure and were difficult to implement on different computers or architectures Killiough1995 (); Wu2002 ().
On the other hand, the demand for modeling capability has increased rapidly in recent years with an increase in computational efforts. More complex geological, physical and chemical features are modeled through reservoir simulations to assess new exploration and production technologies, such as enhanced recovery processes. In addition, the traditional serial simulators have reached their simulation capability limits. The high performance simulation technology has been progressively viewed as an important, alternative modeling approach to solve largescale simulation problems with multimillion and even multibillion block models Edwards2012a ().
Simulation of multiphase flow involves the discretization of conservation and constraint equations on a grid that is constructed from a set of blocks in the domain of interest. Finite difference methods (FDM) are widely adapted by commercial software due to their simplicity and efficiency, while finite volume methods (FVM) are also favoured by many engineering simulations due to their excellent conservation property and ability to handle complex geometry. In either case, their discretization procedure is designed for consistency, such that the errors incurred will vanish rapidly as gridblock sizes approach zero. Nevertheless, the total number of grid blocks that can be handled depends on the capacity of the available computing hardware and is necessarily limited. Clearly, this constraint also applies to the need for resolution of smallscale phenomena. For example, the number of grid blocks required to solve water front propagation phenomena in a secondary recovery model may not be affordable in the context of a complete reservoir simulation and alternative strategies, such as an adaptive grid system, may be required.
Initially, hydrocarbons are displaced from a reservoir by the natural reservoir energy, as the formation pressure is considerably higher than the bottomhole pressure inside a wellbore. However, as the formation pressure declines because of production, the recovery stage reaches its economic limit at a too low production rate, or too high proportions of gas or water in the production stream. This stage is named the primary production that produces only a small percentage of the initial hydrocarbons in place, typically around 10% for oil reservoirs. The second stage of hydrocarbon production maintains reservoir pressure and displaces hydrocarbons toward a wellbore by injecting an external fluid, such as water or gas through injection wells. However, when considerable amounts of injected fluid are produced from production wells, the production is no longer economical. About additional 15% to 40% of the original oil can be recovered after the secondary recovery method. Due to unfavorable mobility ratios and reservoir heterogeneities, early breakthrough during the secondary recovery process prevents the injected water from sweeping the oil efficiently. However, it is often not feasible to change the properties of the displaced fluid or the permeability to the displaced fluid. Most mobility control processes of current interest involve addition of chemicals to the injected fluid. These chemicals increase the apparent viscosity of the injected fluid and/or reduce the effective permeability to the injected fluid. The chemicals used are primarily polymers when the injected fluid is water Green1998 (), which are the most popular choice due to their lower costs compared to other types of additives. Other chemicals, such as surfactant or alkali, can also be injected alternately or simultaneously to enhance oil recovery Lake1992 (), but they are beyond the scope of this research.
By applying polymer, the viscosity of the water phase is increased, and, as a result, the mobility of the water phase is reduced, which results in a more favorable fractional flow curve and then leads to a more efficient sweeping pattern and reduced viscous fingering. The mobility reduction of the injected water is due to two main effects. First, the viscosity of a polymer solution is higher than that of pure water (the viscosity of a polymer solution increases with raising polymer concentration). Second, the permeability to water is reduced after the passage of a polymer solution through rock materials (the permeability to oil is, however, largely unaffected). Both effects reduce the water mobility while the oil mobility is unaltered.
Polymer flooding holds a bright future because it can improve the area swept efficiency not only in the macro scale but also in the micro scale. The first polymer flooding application was reported in 1964 Pye1964 (); Sandiford1964 (). The development of polymer flooding boomed in the US during the 1970s and 1980s with several polymer flooding projects Lake1992 (). However, it declined in the late 1980s because of low oil prices. During the middle 1990s, polymer flooding was resumed in China to large extent. Especially, oil production from polymer flooding contributed to 22.3% of the total oil production in the Daqing oilfield by 2007 Wang2008 (); SaboorianJooybari2016 (). This significance has attracted the petroleum industry’s interest in using reservoir simulators as tools for reservoir evaluation and management to minimize operation costs and increase the process efficiency Chen2007b (); Chen2006a (). Reservoir simulators with special features are needed to represent coupled chemical and physical phenomena present in polymer processes.
Bondor Bondor1972 () presented the development of a threephase, fourcomponent, compressible, finite difference polymer simulator. The model represented a polymer solution as a fourth component that was included in the aqueous phase and was fully miscible with the water phase. Adsorption of polymer was represented as well as the permeability reduction of the water phase. An implicit pressureexplicit saturation (IMPES) procedure was used to solve the coupled system. Lutchmansingh Lutchmansingh1987 () extended this simulator by solving pressure and saturation distributions simultaneously and polymer concentration explicitly. Based on Lutchmansingh’s work, AbouKassem AbouKassem1996a () eliminated nonrelevant equations and unknowns by properly ordering the set of all equations and unknowns, thus providing significant savings in CPU time. Chang Chang1990a () implemented a thirdorder finite difference method to capture a physical dispersion effect which is normally smeared by artificial numerical dispersion. An IMPEC (implicit pressureexplicit concentration) scheme was adapted to solve an isothermal, threedimensional, miscibleflooding compositional model. The simulator is well known as UTCOMP. On the other hand, its run time can be tremendous because of a timestep restriction on the IMPEC form, if a large number of gridblocks are necessary for either simulation of a large reservoir or refinement of a model for more accurate simulation. To overcome this computational limitation, a fully implicit formulation has been adapted in the development of our simulator.
To capture finescale phenomena and optimize a polymer flooding process, largescale reservoir simulations with finescale grids are required. A parallel polymer flooding reservoir simulator has been developed to address these issues. In this paper, the mathematical model of polymer flooding is introduced, including the conservation laws for water, oil and polymer, mechanisms of polymer flooding, and well modelling. Numerical methods are presented. The upstream finite difference (volume) method is applied to discretize the model equations. The standard Newton and inexact Newton methods are applied for their highly nonlinear systems, and an algorithm for the inexact Newton method is introduced. Linear systems from polymer flooding are illconditioned, especially when a reservoir has heterogeneous porosity and permeability. In this case, the linear systems are difficult to solve. In our simulator, a multistage preconditioner is employed to speed up the system solution. Parallel implementations are also introduced. Different polymer flooding cases are used to illustrate the accuracy and scalability of our simulator. The results show that this polymer flooding simulator has good scalability and largescale reservoir models can be simulated.
2 Mathematical Model
The twophase oil and water model is applied, and a temperature change is not considered here. The oil component is assumed to stay in the oil phase, the water stays in the water phase, and polymer only distributes in the water phase. The following sections will present a short summary of all related mathematical models for rock, fluids and well handling.
2.1 Rock Model
When considering porous media at the macroscale, the flow is governed by volume averaged equations. Each computational block contains both solid and pore space which is filled with fluids, such as gas, oil and water. The percentage of pore space, which is called porosity, is defined as
where is the volume of the pore space and is the volume of a block. Porosity is a function of pressure (and temperature), and it can be modelled by the following equation:
(1) 
where is the compressibility factor of the reservoir, is pressure, and is the reference porosity at the reference pressure .
2.2 Fluid Model
The notion of saturation is introduced to define the ratio of the volume of phase to the pore space in a block:
The saturations of the oil phase () and the water phase () satisfy the following relationship:
(2) 
Darcy’s law is applied to handle the relationship among flow rates of a phase, reservoir properties, fluid properties and pressure in a reservoir, which is described as
(3) 
where is a crosssectional area in a flow direction, is a pressure difference, is the viscosity of a fluid, and is the length of a porous medium in the flow direction. is the effective permeability for the given phase, which is the product of absolute permeability and relative permeability . is defined as a tensor with respect to all the , and directions; mostly, it is a diagonal tensor: . Darcy’s law can also be rewritten, with Darcy’s velocity ,
(4) 
With gravity, the mass of each phase satisfies the following conservation law Chen2007 ():
(5) 
where is the phase density, is the source term that models the mass changes caused by injection or production wells, is the gravity, is the depth of a block, and stands for the relative permeability for the phase. In addition, when polymer exists in the water phase, the mass conservation law for the water phase becomes
(6) 
where is the permeability reduction factor caused by polymer and is the viscosity of a waterpolymer solution. The definitions of and will be introduced later.
The water phase pressure, , and the oil phase pressure, , are related by
(7) 
The pressure difference is called the capillary pressure, which usually depends on the saturations of the phases in porous media and is measured by lab experiments. If saturation and any phase pressure are known, the other phase pressure can be calculated by the above formula.
2.3 Polymer Model
The flow of polymer is assumed to act as a component dissolved in the water phase, which is modeled by the following equation:
(8) 
where is the concentration of the polymer in the water phase and is the polymer adsorbed by the reservoir.
When polymer molecules flow through porous media, part of them are restricted in pores, where only water or brine is allowed to pass by with a reduced mobility. As the polymer solution interacts with the reservoir rock, polymer is adsorbed or desorbed from the rock surface; this mechanism is known as polymer retention. There are two mechanisms during the polymer retention process, which are separated as adsorption of the polymer on rock surfaces and entrapment of polymer molecules in small pore space. Both these mechanisms increase the resistance of flow. These effects are modeled by reducing the permeability of the rock to water.
The long chains of polymer molecules can flow into a large pore opening and get trapped when the other end has a smaller opening. Entrapment can also take place when the flow is restricted or stopped. When this happens, the polymer molecules lose their elongated shape and coil up. Desorption of the polymer from the reservoir rock can also take place if sufficient polymer has already been adsorbed above a residual sorption level. It is difficult to quantify what percentage of injected polymer is adsorbed and what percentage is trapped in small pore spaces since only the produced polymer concentration can be measured. Both these mechanisms result in a loss of polymer to the reservoir.
The adsorption process causes a reduction in the permeability of the rock to the passage of the aqueous phase and is directly correlated to the adsorbed polymer concentration. The reduction factor, , is a function of polymer adsorption and the residual resistance factor (RRF), which is expressed as
(9) 
where is the cumulative adsorption of polymer per unit volume of the reservoir rock and represents the maximum value of , which denotes the maximum adsorptive capacity of polymer per unit volume of the reservoir rock. Both RRF and are functions of the reservoir rock permeability.
Assuming equilibrium sorption with the reservoir rock, the sorption phenomenon can be described as a function of polymer concentration only:
(10) 
This relationship is specified in the form of a table.
Not only is the rock permeability to water reduced after the passage of a polymer solution through porous media, but also the viscosity of the polymer solution is higher than that of pure water. Small concentrations of polymer, on the order of a few hundred to a few thousand ppm (by weight), increase the viscosity of an aqueous solution significantly Green1998 ().
The ToddLongstaff technique is used to calculate the effective viscosity that incorporates the effect of physical dispersion at the leading edge of a slug and also the fingering effect at the rear edge of the slug Lie2012 (). The viscosity of a fully mixed polymer solution, denoted by , rises as the polymer concentration () increases. The viscosity of the solution at the maximum polymer concentration is also specified that is denoted by . Then the effective polymer viscosity is taken to be
(11) 
where is the ToddLongstaff mixing parameter. The mixing parameter is useful in modeling the degree of segregation between water and the injected polymer solution. If , then the polymer solution and water are fully mixied. If , the polymer solution is completely segregated from the water.
A partially mixed water viscosity is calculated in an analogous manner using the fully mixed polymer viscosity and the pure water viscosity
(12) 
The effective water viscosity is calcualted by the partially mixed water viscosity and the effective polymer viscosity as a harmonic average:
(13) 
where is the effective saturation of the injected polymer solution within the total aqueous phase in a block.
The mixing of polymer and water modifies the solution viscosity as well. Since polymer has higher viscosity compared to pure water, no matter which mixing rule is selected, the mixture viscosity increases as a function of polymer concentration in the solution. Two commonly used mixing rules are used, which include a linear mixing rule:
(14) 
and a nonlinear mixing rule:
(15) 
where is a parameter dependent on polymer concentration given by
(16) 
A higher water viscosity and a reduction in permeability will result in an increase in the resistance to flow, and divert the polymer solution toward areas unswept by water. This mechanism is wellknown as mobility control. Directly, the wateroil mobility ratio is reduced to close to unity or less. Then the volumetric sweep efficiency is improved and higher oil recovery is achieved compared to conventional water flooding.
As mentioned above, polymer molecules can flow into large pore openings. However, there are also small openings which are not contacted by polymer molecules. To describe this phenomenon, an inaccessible pore volume (IPV) is used to measure all the pore space that may not be accessible to polymer molecules. The presence of IPV causes the polymer solution to travel at a greater velocity than inactive tracers embedded in water. This chromatographic effect is modeled by assuming that the IPV is constant for each rock type and either does not exceed the corresponding irreducible water saturation or is independent of the water saturation. The concept of IPV allows a polymer solution to advance and displace oil at a faster rate than predicted on the basis of total porosity.
2.4 Wellbore Models
A numerical simulation of fluid flows in petroleum reservoirs must account for the presence of wells. They supply a set of realistic boundary conditions for computations of pressure distributions Dumkwu2012 (). The fundamental task in modeling wells is to model flows into/from a wellbore accurately and to develop accurate well equations that allow the computation of the bottom hole pressure with a given production or injection rate, or the computation of a rate with known pressure Chen2006a ().
Peaceman Peaceman1978 (); Peaceman1983a () associated a steadystate pressure for an actual well with the computed pressure at a grid block through the concept of an equivalent radius . If a well was completed in more than one grid block, a well index (WI) was introduced to account for well pressure losses within the grid blocks due to the radial inflow into the well. The well index depends on the geometry of a grid block, location and orientation of the well segment in that grid block, anisotropic reservoir property and a skin factor STARS ():
(17) 
where is the well fraction that is evaluated by the angle open to flow and varies due to the well position in a grid block. It equals 1 for a well going approximately through the center of a grid block. represents a grid block thickness along the well direction, and is the grid block thickness factor. The current completion length in the current grid block is the product of and . is the geometric average permeability and estimates the formation’s absolute permeability perpendicular to the well direction. denotes the skin factor, which may also differ from one perforated block to another within the wellbore. This is especially true if different perforation densities and intervals exist within each individual simulation layer. is the wellbore radius. The equivalent radius and formation absolute permeability are computed according to the wellbore direction and the discretization procedures. For instance, if a well is parallel to the direction in a Cartesian grid, then
(18)  
Similar to the well fraction , the factor in (18) depends on the geometry of a grid. It equals 0.249 for a well going approximately through the center of a grid block. Detailed information can be found in STARS (); Palagi1994a ().
The flow rate, , for the phase in a perforated grid block is the product of the well index, the fluid mobility and the drawdown pressure Chen2007 (); Fung2005 ():
(19) 
The wellbore pressure at each grid completion () is different from one layer to another, depending on the existing pressure drop in a wellbore. It is calculated by the hydrostatic pressure difference drawn from the average density of the fluid mixture in the wellbore:
(20) 
where is the fluid unit weight which depends on the fluid mixture density in the wellbore and is the reference bottom hole pressure at reference depth Coats1980 (); Fung2005 ().
To optimize oil production and to reduce operation costs, various well operations may be employed at any time, such as fixed bottom hole pressure, a fixed oil production rate, a fixed water production rate, a fixed water injection rate, or a fixed liquid production rate. When the fixed bottom hole pressure well operation is applied to a well, the constraint for the well is described as
(21) 
where is a constant. The fixed water rate condition is the following equation:
(22) 
where is a constant. For the fixed oil rate production operation, the constraint is
(23) 
where is a fixed constant. For the fixed liquid production rate operation, the production of oil and water is fixed, and its constraint equation is
(24) 
where is a constant.
3 Numerical Methods and Parallelization
The reservoir model for polymer flooding is highly nonlinear, which is hard to solve analytically. In this paper, numerical solutions are obtained by the fully implicit method thanks to its unconditional stability. Figure 1 shows the details of the solution flowchart, which includes the following steps:

Loading model. In this step, a model file is loaded, which contains information for reservoirs, such as permeability, porosity and geometry, for oil, water and polymer, such as density, viscosity, and concentration, for well operations, and for numerical parameters. The model may have hundreds of parameters.

Grid generation and distribution. The model file defines a grid, such as dimensions in the , , and directions, sizes of each gridblock, and coordinates. Since the simulation is parallel, the grid must be distributed to each MPI, and a communication structure must be set up.

Initialization. This step sets the initial pressure, saturations of oil and water, concentration of polymer, bottom hole pressure, porosity and permeability of the reservoir, and other properties, such as relative permeability, density, viscosity, and well index.

Time discretization and time step selection. A time step is dynamically selected, which satisfies some conditions:

If Newton methods fail, the time step will be cut.

In one time step, well operations keep unchanged.

The time step has a maximal value, which is set by the model file.

The algorithm always attempts to increase a time step to reduce simulation time.


Newton iteration. In each time step, an nonlinear system is solved by the Newton method, which converts an nonlinear system to a linear system, . Inside this step, the properties for the reservoir and fluids must be computed, such as porosity, density and viscosity.

Linear iteration. This step solves the linear system . A proper linear solver, a preconditioner and solution parameters must be chosen, which can be input by the model file.
The MPI (Message Passing Interface) is applied to handle communication among computation nodes. When calculating properties, such as transmissibility, neighboring information is always required. To develop a scalable parallel application, communications should be minimized. In reservoir simulations, the communication pattern is determined by a grid distribution.
Grid partitioning algorithms aim to minimize the idle time and communication flow between different processors by dividing the blocks equally among partitions and minimizing the number of partition spanning edges. The quality of the partitioning plays a crucial role in the performance of applications. Reservoir simulation applies gridbased numerical methods, such as the finite volume and finite difference methods, to approximate a governing system. Grid blocks have a higher possibility to communicate with each other when they are closer to each other. Based on this assumption, we have introduced a modified Hilbert spacefilling curve partitioning method Liu2017a () in our inhouse simulator, which has shown promising loading balance and excellent scalability.
After the grid partitioning subroutine occurs, the subregion topology graph and block level connection lists are set up. The block level connection lists record the geometry information both in the global and local indices. Here, the global index is normally defined by the natural ordering. It labels the grid blocks in the , and directions, respectively. Once the subregion configuration is fixed, the dataset is read and scattered into each processor. A parallel IO algorithm through MPIIO is implemented to read the data files.
A set of data structures have been designed to store distributed data, such as DOF (Degrees of Freedom), VEC, MAT, SOLVER, and SOLVER_PC prsi (). A DOF is defined on a grid, which can be defined to store blockbased data, such as porosity, density and viscosity, and on a well, which can be defined for bottom hole pressure, a well rate, and a well index. The VEC and MAT are distributed vectors and matrices. The reader refers to our previous paper prsi () for more details of parallelization.
3.1 Time Discretization
In reservoir simulation, to study well performance and optimize production, a simulation period is usually long, such as over 10 years. Here the backward Euler method is applied to discretize a time derivative. Let be a function, be the solution of at a time step , and be an nonlinear system. For the following differential equation:
(25) 
the backward Euler method computes the numerical solution by
(26) 
where is a time step. This method is implicit, and the above system is still nonlinear, which must be solved at each time step.
3.2 Spatial Discretization
The oil phase pressure, water saturation, polymer concentration and well bottom hole pressure are chosen as the primary unknowns, and other unknowns are functions of these primary unknowns.
When fluids move in a reservoir, there may be fluid exchange in two neighboring blocks. The term, transmissibility, is defined to describe the amount of fluid exchange. Here, let be any space direction and be the area of a face in the direction; then the transmissibility term for phase is defined as
(27) 
where is the cell length along the direction, is the permeability, is the relative permeability, is the viscosity, and is 1 for the oil phase and for water and polymer.
The transmissibility is defined on each face of a block. For any two neighboring blocks, since they share a face, the value of the transmissibility term is the same for these two blocks. Different weighting schemes must be applied to average different properties at an interface for them to make a physical sense. The geometric properties, such as a cell length and a cross area , and rock permeability are harmonically averaged. Since the fluid properties, such as and , do not change much, they are averaged by a pore volumeweighted arithmetic average scheme. An upstream weighting technique is applied to calculate a relative permeability, which means that a relative permeability at the interface of two neighboring blocks is calculated from the block that has a higher potential. For example, the relative permeability at the interface is defined as
(28) 
Other higherorder upstream weighting techniques can also be used for relative permeabilities Chen2007 ().
3.3 Newton Methods
The standard Newton method is the usual way to solve nonlinear systems. In our implementation, an inexact Newton method inexactNewton () is also applied, which relaxes the termination tolerance to reduce computation time. Its algorithm is described in Algorithm 1. The only difference from the standard Newton method is that the latter uses a fixed termination tolerance for a linear solver while the inexact Newton method uses a dynamic tolerance, which is computed automatically by equation (30).
(29) 
is the dynamic termination tolerance for the linear system . The choice of this parameter is designed to speed the convergence of the Newton method, to avoid oversolution of a linear system, and to reduce computation time. Three popular algorithms are provided as follows:
(30) 
where is the residual vector of the th Newton iteration and is the righthand side vector. The third formula is applied in our simulator, where is 0.5 and is . is also tailored to meet this contition: .
3.4 Linear Solver
A Jacobian matrix is nonsymmetric and highly illconditioned. The Krylov subspace solvers are applied to solve the linear system . In real application, a preconditioner is always applied to solve an equivalent linear system . A family of scalable CPR methods CPRH () have been developed to handle linear systems from reservoir simulations, and the CPRFP method in CPRH () is used as the preconditioner.
The system has four unkonws, oil phase pressure (), water saturation (), polymer concentration () and well bottom hole pressure ), which are vectors. There are two common strategies to arrange the unknown , which are pointwise and blockwise, respectively. For the pointwise strategy, is written as
(31) 
while, for the blockwise strategy, is written as
(32) 
where is the number of grid blocks and is the number of wells. In real simulations, the preconditioner using the blockwise strategy has better convergence than that using the pointwise strategy. However, if the pointwise strategy is applied, the Jacobian matrix J has a clear block structure as
(33) 
where is the matrix corresponding to the oil phase pressure, is the matrix corresponding to the water saturation, is the matrix corresponding to the polymer concentration, is the matrix corresponding to the well bottom hole pressure, and other matrices are coupled terms.
If we define a restriction operator from to , then its formal formula can be written as
(34) 
A prolongation operator from to can be defined as
(35) 
The preconditioning linear system must be solved in each iteration. The CPRFP method CPRH () can be described by Algorithm 2, where the first step is to solve an approximate solution using restricted additive Schwarz (RAS) method and the third step is to solve the subproblem by algebraic multigrid method (AMG). It is well known that RAS method and AMG method are scalable for parallel computing, so the CPRFP method is also scalable.
4 Numerical Experiments
4.1 Model Validation
In this section, we present several results to validate the application of our simulator (BOS) in different polymer flooding projects. First, a homogeneous, onedimensional polymer flooding case with 15 blocks is presented. The polymer concentration varies in different production stages. Second, we consider a twodimensional polymer flood with multiple constraints on the injection wells. Finally, we consider a threedimensional case with a stratified reservoir. It has a number of layers with significant or complete lateral continuity. The layer permeabilities vary greatly and adjacent layers communicate vertically. Simulation results show that water advances rapidly in highly permeable layers and slowly in tight layers as the driving fluid to displace oil.
ECLIPSE100, Schlumberger, is a fullyimplicit, threephase, threedimensional, general purpose black oil simulator that can model different chemical EOR (enhanced oil recovery) processes, including polymer and surfactant flooding. It can be run in fully implicit and adaptive implicit modes, and is widely considered as a reference and benchmark standard, which has been successful in reproducing laboratory measurements and field applications. It is widely used in the oil industry to evaluate different production processes. Our simulator will be compared against it to show the accuracy.
4.1.1 Polymer flood in a onedimensional geometry
Here we consider a special case of polymer flood in onedimension geometry. The reservoir contains 15 grid blocks along the direction where the properties are uniformly distributed at each grid block with porosity and intrinsic permeability . A polymer adsorption curve and water viscosity multiplier are displayed in Figure 2, and the fluid properties are summarized in Table 1.
Model parameter  Value 
Spatial grid size ()  100 
Initial resident water saturation  0.4 
Initial polymer concentration ()  0.1 
Water  
Mass density ()  1025.18 
Formation volume factor ()  1.0 
Compressibility ()  3.03e6 
Viscosity ()  0.5 
Oil  
Mass density ()  832.96 
Formation volume factor ()  1.0 
Compressibility ()  1.0e5 
Viscosity ()  0.5 
Rock  
IPV  0.15 
RRF  2.67 
Maximum polymer adsorption (kg/kg)  0.0035 
We inject water with a fixed constant volume rate at by the injection well located at the first grid block. The production well is perforated through the 15th grid block, where a liquid volume well constraint was operated at . The polymer injection concentration varies along the production process. At the beginning, only pure water is injected for , and then polymer is added with concentration at to improve the volumetric sweep efficiency. Finally, we stop injecting polymer after injection, and simulation is terminated after total .
In Figure 3, the oil and water production rates are compared with those from SchulumbergerEclipse that are represented by the green line, while our numerical solutions are plotted in red color and line markers. After about of water injection, the water front has broken through, and the oil production rate reduces until the whole reservoir is flooded. A good agreement is found in Figure 3 with some minor differenced due to numerical diffusion and different ways to approximate the water phase viscosity.
4.1.2 Polymer flood in a twodimensional geometry
Next, we simulate a polymer flood case on a grid and the formation is initially fully saturated with oil with porosity and intrinsic permeability . The fluid properties are summarized in Table 2.
Model parameter  Value 
Spatial grid size ()  75 
Initial polymer concentration ()  0. 
Water  
Mass density ()  64 
Formation volume factor ()  1.0 
Compressibility ()  3.03e6 
Viscosity ()  0.5 
Oil  
Mass density ()  52 
Formation volume factor ()  1.0 
Compressibility ()  1.0e5 
Viscosity ()  2 
There is no polymer in the reservoir initially. An injection well locates at the left corner of the grid and a production well locates at the other end of the diagonal. Water and a polymer solution are injected at a maximum injection rate of with polymer concentration at . At the same time, the injection well is constrained with its bottom hole pressure, with no more than . There is no constraint on the production rate, but the bottom hole pressure of the production well is fixed at . After production, there is no polymer injected into the formation. The simulation is terminated after total .
The same case was carried out by SchulumbergerEclipse, and the oil and water production rates are compared in Figure 4 by the red and green colors that represent our simulator and Eclipse, respectively. The same production pattens are shown with different simulators, and the difference between them is negligible.
4.1.3 Polymer flood in a threedimensional geometry
In this case, we extend the twodimensional case into a threedimensional problem. The same as in the twodimensional problem, two vertical wells locate at the two ends of the diagonal of a plane and act as an injection well and a production well. Both wells have a multiple layer perforation that is different from the twodimensional case. For multiple perforated layer wells, a reference bottom hole pressure is selected at a datum depth. The fluid density is calculated at that depth and the initial pressure in the reservoir is determined by marching down the reservoir in small steps by recalculating the density at each step. It is treated explicitly through the Newton iteration. This algorithm is more complicated when more than one phase is present. Starting at the datum depth, the hydrostatic pressure for the datum phase can be calculated by marching up and down the reservoir. Pressures in the other phases can then be determined at the contact depths, and then the hydrostatic pressures can be determined throughout the reservoir by marching up and down again. Once the phase pressures are known, the phase saturations can be determined at each depth so that the hydrostatic pressure variation is balanced by the capillary pressure between the phases.
The agreement of our simulator with commercial software is shown in Figure 5 for the oil and water production rates.
4.2 Performance Test
In this section, the performance of our simulator is reported in terms of speedup. The speedup, , is defined as the ratio of the elapsed time when executing a simulation by cores to the execution time when processors are used:
(36) 
Independent timings are performed in different running cases by our simulator. In the ideal scenario, the amount of time reduces linearly with the number of CPUs rising. From a performance point of view, our simulator is demonstrated to provide significant speedup with respect to multicore CPUs through our tests.
4.2.1 Low resolution reservoir model
The first performance test is carried out on a polymer flooding project with grid blocks in the , and directions, respectively, on a IBM System x3750 M4 Intel Xeon processorbased, entrylevel 4socket server. The server consists of 32 3.3GHz SMP CPU cores with 512 GB local RAM and 10GbE networking options in a 2U form factor.
The reservoir properties, such as porosity and permeability, vary along different layers, as Table 3 shows. The initial reservoir pressure is at a reference depth of . The simulation was based on water flood by a vertical and horizontal well patten for , followed by polymer flood for almost . Four vertical injection wells are distributed at the corners of the reservoir, while a horizontal producer is drilled at the center of the reservoir.
Layer  Porosity  Permeability (mD)  

X  Y  Z  
1  0.17393  326.4  980.6  163.4 
2  0.1694  445.3  1335.8  222.6 
3  0.25714  148.9  446.6  74.4 
4  0.17344  118.8  356.4  59.4 
5  0.1187  71.2  213.6  35.6 
The results reported in Figure 6 demonstrate the performance of our simulator. Simulations were performed from 32 (the reference) to 256 cores. It is shown in Figure 6 that the simulator approximately has the linear speedup performance.
In the case of water flood, the water phase advances through the heterogeneous formation much faster compared to the polymer flood. This is due to the fact that pure water being less viscous is more mobile. On the other hand, the polymer flood is able to achieve a much more efficient sweep of the reservoir. Although this makes polymer flooding a slower recovery process compared to waterflooding, the net oil recovery can be substantially more in polymer flood. The polymer flood results clearly show that a larger area has been swept and also the higher water saturation levels have penetrated deeper into the domain compared to the water flood.
4.2.2 High resolution reservoir model
The second performance test is carried out by Graham cluster supplied by Compute Canada. The Graham consists of 35,520 cores and 320 GPU devices, spreaded across 1,107 nodes. All nodes except bigmem3000 have Intel E52683 V4 CPUs, running at 2.1 GHz. This test was run on the nodes with 128GB memory where a lowlatency highbandwidth Infiniband fabric connects all nodes and scratch storage. Nodes configurable for cloud provisioning also have a 10Gb/s Ethernet network, with 40Gb/s uplink to scratch storage. The same reservoir description as in the last case is used except with the refined grid size of . Simulations were performed from 128 (the reference) to 2048 cores.
This example tests the scalability of the platform as well as the solution of linear systems (including the restarted GMRES iteration solver and the CPR preconditioner). Figure 7 shows that when processors are doubled, the elapsed time costed by the platform and linear solver is cut by half, which means that our simulator has excellent scalability. This example also shows that the simulator can deal with more accurate geological models at multimillion and even multibillion grid blocks under accepted run time if there are enough computation resources available.
5 Conclusions
The results of this work demonstrate that coupled fluid flow and polymer adsorption phenomena can be effectively simulated and distributed among multicore CPUs. High performance computing techniques are applied through our inhouse simulator, which has the capability to build and solve largescale reservoir systems, especially reservoir models with high resolutions and complexity. The rapid advancements in parallel computer hardware and software technologies are adapted to improve the efficiency of polymer flooding processes. A satisfactory parallel efficiency of our simulator is demonstrated for a complex field application with up to 27 million grid blocks by 2048 cores. Almost linear speedup is displayed.
6 Acknowledgments
The support of Department of Chemical and Petroleum Engineering, University of Calgary and Reservoir Simulation Group is gratefully acknowledged. The research is partly supported by NSERC/Energi Simulation Chair, Alberta Innovate Chair, IBM Thomas J.Watson Research Center, and the Frank and Sarah Meyer FCMG Collaboration Centre for Visualization and Simulation. The research is also enabled in part by support provided by WestGrid (www.westgrid.ca), SciNet (www.scinethpc.ca) and Compute Canada Calcul Canada (www.computecanada.ca).
References
References
 (1) D. Coumou, S. Matthäi, S. Geiger, T. Driesner, A parallel FEFV scheme to solve fluid flow in complex geologic media, Computers and Geosciences 34 (12) (2008) 1697–1707. doi:10.1016/j.cageo.2007.11.010.
 (2) A. H. Dogru, Megacell Reservoir Simulation, Journal of Petroleum Technology 52 (5). doi:10.2118/57907MS.
 (3) K. H. Coats, Reservoir Simulation, in: Petroleum Engineering Handbook, Society of Petroleum Engineers, 1987, Ch. 48, pp. 20–48.
 (4) J. E. Killough, The application of parallel computing to the flow of fluids in porous media, Computers {&} Chemical Engineering 19 (67) (1995) 775–786. doi:10.1016/00981354(94)000808.
 (5) Y.S. Wu, K. Zhang, C. Ding, K. Pruess, E. Elmroth, G. S. Bodvarsson, An efficient parallelcomputing method for modeling nonisothermal multiphase flow and multicomponent transport in porous and fractured media, Advances in Water Resources 25 (3) (2002) 243–261. doi:10.1016/S03091708(02)000064.
 (6) D. A. Edwards, D. Gunasekera, J. Morris, G. Shaw, K. Shaw, P. A. Fjerstad, J. Kikani, J. Franco, V. Hoang, L. Quettier, Reservoir Simulation : Keeping Pace with Oilfield Complexity, Oilfield Review 23 (4) (2012) 4–15.
 (7) D. W. Green, G. Willhite, Enhanced Oil Recovery, Society of Petroleum Engineers, 1998.
 (8) L. Lake, R. Schmidt, P. Venuto, A niche for enhanced oil recovery in the 1990s, Oilfield Review 4 (1) (1992) 55–61.
 (9) D. J. Pye, Improved secondary recovery by control of water mobility, Journal of Petroleum Technology 16 (08) (1964) 911–916.
 (10) B. Sandiford, Laboratory and field studies of water floods using polymer solutions to increase oil recoveries, Journal of Petroleum Technology 16 (08) (1964) 917–922.
 (11) D. Wang, R. S. Seright, Z. Shao, J. Wang, Key aspects of project design for polymer flooding at the daqing oil field, SPE Reservoir Evaluation & Engineering (2008) 1117–1124.
 (12) H. SaboorianJooybari, M. Dejam, Z. Chen, Heavy oil polymer flooding from laboratory core floods to pilot tests and field applications: Halfcentury studies, Journal of Petroleum Science and Engineering 142 (2016) 85–100.
 (13) Z. Chen, Y. Ma, G. Chen, A sequential numerical chemical compositional simulator, Transport in Porous Media 68 (3) (2007) 389–411.
 (14) Z. Chen, G. Huan, Y. Ma, Computational methods for multiphase flows in porous media, Society for Industrial and Applied Mathematics, 2006.
 (15) P. Bondor, G. Hirasaki, M. Tham, Mathematical Simulation of Polymer Flooding in Complex Reservoirs, Society of Petroleum Engineers Journal 12 (05) (1972) 369–382. doi:10.2118/3524PA.
 (16) P. M. Lutchmansingh, Development and application of a highly implicit, multidimensional polymer injection simulator, Ph.D. thesis, The Pennsylvania State University (1987).
 (17) J. H. AbouKassem, M. E. Osman, A. M. Zaid, Architecture of a multipurpose simulator, Journal of Petroleum Science and Engineering 16 (4) (1996) 221–235.
 (18) Y.B. Chang, G. A. Pope, K. Sepehrnoori, A higherorder finitedifference compositional simulator, Journal of Petroleum Science and Engineering 5 (1) (1990) 35–50.
 (19) Z. Chen, Reservoir Simulation Mathematical Techniques in oil recovery, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2007.
 (20) K. A. Lie, H. Nilsen, A. Rasmussen, X. Raynaud, An unconditionally stable splitting method using reordering for simulating polymer injection, ECMOR 2012  13th European Conference on the Mathematics of Oil Recovery (2012) 10–13.
 (21) F. A. Dumkwu, A. W. Islam, E. S. Carlson, Review of well models and assessment of their impacts on numerical reservoir simulation performance, Journal of Petroleum Science and Engineering 8283 (2012) 174–186. doi:10.1016/j.petrol.2011.12.005.
 (22) D. W. Peaceman, Interpretation of WellBlock Pressures in Numerical Reservoir Simulation(includes associated paper 6988 ), Society of Petroleum Engineers Journal 18 (3) (1978) 183–194. doi:10.2118/6893PA.
 (23) D. W. Peaceman, Interpretation of WellBlock Pressures in Numerical Reservoir Simulation With Nonsquare Grid Blocks and Anisotropic Permeability, Society of Petroleum Engineers Journal 23 (3). doi:10.2118/10528PA.
 (24) CMG, STARS User’s Guide (2013).
 (25) C. L. Palagi, K. Aziz, Modeling Vertical and Horizontal Wells With Voronoi Grid, SPE Reservoir Engineering 9 (1) (1994) 15–21. arXiv:24072, doi:10.2118/24072PA.
 (26) L. S. K. Fung, H. Su, C. T. Tan, K. Hemanthkumar, J. A. Pita, A FullyImplicit FullyCoupled Well Model for Parallel MegaCell Reservoir Simulation, in: SPE Technical Symposium of Saudi Arabia Section, 1416 May, Dharan, Saudi Arabia, 2005.
 (27) K. H. Coats, InSitu Combustion Model, SPE Society of Petroleum Engineers of AIME 20 (December) (1980) 533–554. doi:10.2118/8394PA.
 (28) H. Liu, K. Wang, B. Yang, M. Yang, R. He, L. Shen, H. Zhong, Z. Chen, Dynamic Load Balancing Using Hilbert SpaceFilling Curves for Parallel Reservoir Simulations, in: SPE Reservoir Simulation Conference, 2022 February, Montgomery, Texas, USA, 2017. arXiv:1708.01365.
 (29) H. Liu, K. Wang, Z. Chen, K. Jordan, J. Luo, H. Deng, A parallel framewrok for reservoir simulators on distributedmemory supercomputers, in: SPE176045MS,SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, Nusa Dua,Â Indonesia, 20 â 22 October, 2015.
 (30) T. Chen, N. Gewecke, Z. Li, A. Rubiano, R. Shuttleworth, Y. Yang, X. Zhong, Fast computational methods for reservoir flow models (2009).
 (31) H. Liu, K. Wang, Z. Chen, A family of constrained pressure residual preconditioners for parallel reservoir simulations, Numerical Linear Algebra with Application 23 (1) (2016) 120–146.