Dynamic PVT model for CO2-EOR black-oil simulations

A well-planned CO2EOR operation can help meet an ever-increasing need for energy and at the same time reduce the total CO2 footprint from the energy production. Good simulation studies are crucial for investment decisions where increased oil recovery is optimized and balanced with permanent CO2 storage. It is common to use a compositional simulator for CO2 injection to accurately calculate the PVT properties of the mixture of oil and CO2. Compositional simulations have significantly increased simulation time compared to blackoil simulations. Large simulation studies where many simulations are used either to represent uncertainty and/or optimize the results can thus be unpractical. On the other hand existing blackoil formulations often poorly represent the PVT properties of the oil-CO2 mixtures. We therefore present an extended blackoil formulation with dynamic blackoil properties that depends on the fraction of CO2 in the cell. These properties represents the density and viscosity of the hydrocarbon - CO2 mixture more accurately and thus give results closer to the compositional simulator. The dynamic blackoil functions are calculated from numerical slim-tube experiments based on one-dimensional equation-of-state (EOS) simulations. The same simulations also gives estimates of the minimum-miscibility pressure (MMP). We present examples based on data from the Fifth Comparative Solution Project: Evaluation of Miscible Flood Simulators, but uses CO2 as the injection gas. These examples shows that the new blackoil model gives results that are closer to compositional simulations compared to existing blackoil formulations. The model is implemented in the Flow simulator. The Flow simulator is developed as part of the open porous media (OPM) project and is an openly developed and free reservoir simulator that is capable of simulating industry relevant reservoir models with similar single and parallel performance as commercial simulators.


Introduction
The need to reduce CO 2 emissions opens up new possibilities for utilization of CO 2 . Injection of CO 2 for enhanced oil recovery is an EOR technique with great potential. An early summary study of onshore CO 2 -EOR project in the US shows that between 5-22% of extra oil is produced as a result of CO 2 injection [1]. An important prerequisite for a high recovery rate is that the reservoir is sufficiently swept by the injected CO 2 . A challenge is therefore to achieve a good sweep. Which in practice proves to be difficult since CO 2 both as a separate phase, or mixed with oil, flows more easily than oil without Tor Harald Sandve tor.harald.sandve@norceresearch.no 1 NORCE Norwegian Research Centre AS, Oslo, Norway dissolved CO 2 . Injected CO 2 will therefore take the fastest route to the nearest producers. Density differences may also cause gravitational override so that lower lying oil remains unswept. For onshore USA, good sweep can be ensured by having small distance between the wells, offshore this is demanding due to significant costs associated with drilling new wells. Traditional CO 2 -EOR uses CO 2 from natural sources and thus focus on producing as much oil as possible without injecting too much CO 2 . In light of today's climate challenge, it is interesting to look at opportunities to combine increased oil recovery with CO 2 storage. Projects like Boundary Dam [2] and Petra Nova [3] are examples of this. A prerequisite for the success of such projects is good planning so that both a sufficient amount of extra oil is produced to justify the investments and a considerable amount of CO 2 is stored. Simulations provide important input into such planning, both for optimizing design and for reducing uncertainty regarding investment decisionmaking. Simulations can also be important in order to document CO 2 storage for tax-refunding. Representation of uncertainty in model parameters by simulating multiple realizations is standard for several oil companies. This places great demands on the robustness and speed of the simulation tools. Automatic optimization under uncertainty for robust decisions further increases the requirement for performance [4].
It is common to use compositional simulations for CO 2 -EOR so that the PVT properties of the mixture of CO 2 and oil are accurately represented in the model. With, for example, 10-15 components in the model, this makes the simulations considerably more time-consuming than standard black-oil models. Pseudoization of components can be used to reduce the number of components and thus also the running time of the simulations [5]. A simplified model based on an extension of the black-oil model was proposed in [6]. This model has later been expanded and improved in for instance [7,8]. A mixture of oil and CO 2 can give non-monotonous oil density as a function of CO 2 fraction (See Fig. 9). Density and viscosity in existing extended black-oil models are typically based on a fourth order mixing rule. Such simplified rules can not fully represent the complicated PVT behavior of CO 2 and oil mixture. Simulations reported in [9] however show that the simplified black oil based models produce similar results as the compositional simulation with some adjustments. The aim of this work is to further improve the black-oil based models so that they provide sufficiently accurate results to be useful in CO 2 -EOR studies. Especially for preliminary studies where optimization of different strategies under uncertainty is necessary, the performance improvement obtained by using black-oil based models can be critical for feasibility for the study.
The model we present in this paper builds on the model proposed in [6] in that there is a four component extension of the black-oil model where injected gas is represented by the fourth component. To improve the representation of the PVT properties of the mixture of oil and CO 2 , the traditional black-oil tables are extended to also depend on the CO 2 content. The calculations of the values in the black oil tables need the composition of the oil and gas components in the oil and gas phase. The composition will vary during the injection process as more and more of the light components are washed out by the CO 2 . The usual assumptions about fixed gas and oil compositions that are used in a black-oil model can therefore not be assumed. We therefore choose to use a selection of compositions from a one-dimensional simulation as the starting point for the black-oil tables. In the one-dimensional simulation, CO 2 is injected at one end and produced from the other end corresponding to a slim tube experiment. The compositions used in the black-oil tables will then represent different stages of the gas injection. These simulations depend only on the fluid properties and are thus independent of the reservoir and can be done as a preprocessing step prior to the simulations.
To evaluate the method presented in this paper, we use SPE5 [10] benchmark study. The SPE5 benchmark compares compositional and extended black-oil simulations for a small reservoir under varying miscibility. In the original SPE5 study, the injection gas consists mainly of methane. Although the method developed here is independent of injection gas, we choose to use CO 2 instead of hydrocarbon gas as the example.
The model is implemented in the OPM Flow simulator. Flow is an open-source reservoir simulator that can be downloaded for free (http://www.opm-project.org). It uses industry-standard input/output formats and has comparable single and parallel performance as commercial simulators [11]. The implementation in Flow gives immediate applicability on industrial relevant field scale models.
This article is organized as follows. First, the model and calculation of the CO 2 -dependent black-oil functions are presented. Then follows a series of numerical examples where the black oil based models is compared to compositional simulations. The article ends with a brief summary and evaluation of the model.

Method
The description of the improved extended black-oil model is divided into two parts. The first part describes the equations, the second part the computation of the black-oil functions.

The extended black-oil model
In this section the extended black-oil model formulation we use is described. The three first equation below are the standard black-oil equations, see e.g. [12], and the fourth equation is an extra mass conservation equation for the injected CO 2 component. Subscripts w, o, g are used for water, oil and gas phase respectively. We use capitalised subscripts W, O, G, CO2 for the water, oil, hydrocarbon gas, and CO 2 components.
Here φ is the porosity, Q β is the source term for component β, and S α denotes the saturation of phase α. The formation volume factors B α , the dissolution factors R s and R v , as well as the CO 2 volume fractions x and y are now not only functions of pressure but also functions of the CO 2 fraction z. Details on the computation of the black-oil functions for various pressures and CO 2 fraction z are described below using the SPE5 fluids as an example. The Darcy velocity v α for phase α is given as The permeability K, relative permeability k α , viscosity μ α , phase pressures p α and the gravity g follows the standard definition. The phase density ρ α depends on the phase volume factors and the dissolution factors and are thus function of both pressure and the CO 2 fraction. To complete the system of equations we use tabulated capillary pressure relations and let the saturation for the phases S α sum to one.
The four unknowns we solve for are the reference pressure (p) which typically is the oil pressure, the water saturation S w , the CO 2 fraction z and either S g , R s or R v depending on whether we have both gas and oil phase, only oil phase or only gas phase present.
For a description of the compositional model we use for comparison we refer to the technical manual of the Eclipse simulator [13]. The standard solvent model is described in [7,8,11] For relative permeability, standard three phase formulations are used. High content of CO 2 in the oil phase may affect the relative permeability functions. In the standard solvent model an additional relative permeability function can be supplied to model miscible conditions. The above formulation can be extended similarly, but this is not done for the simulations presented in this manuscript as the paper focus on the PVT modeling.

Computation of the dynamic black-oil functions for the SPE5 fluids
In this section we show how the dynamic CO 2 dependent black-oil functions are computed. We use the SPE5 fluid [10] as an example, but the methodology are applicable for other fluid compositions. The composition of the SPE5 oil is shown in Table 1, and the parameters of the Peng-Robinson equation of state are shown in Tables 2 and 3. We however start by computing the saturation pressure for added CO 2 and the minimum miscibility pressure. Saturation pressure are determined by solving the equilibrium equations Here, χ = [χ 1 , . . . , χ n ] T is the composition vector consisting of the mole fractions χ i of each species i. Further, K i = χ v i /χ l i is the mole fraction ratio and φ i is the fugacity coefficient of species i. Superscripts v and l refer to the vapor and the liquid phase, respectively. Figure 1 shows the saturation pressure of SPE5 oil as function of added CO 2 .
The minimum miscibility pressure (MMP) is an important determination factor for evaluation of CO 2 -EOR. The minimum miscibility pressure of a flooding is defined as Table 3 Interaction parameters at reservoir temperature T = 344.26 K for the SPE5 mixture the lowest pressure for which the injection phase in a onedimensional displacement forms a continuous phase from inlet to outlet. The recovery rate by gas injection strongly depends on the pressure, and for gas injection above the minimum miscibility pressure one can achieve significantly better recovery rate than with other floodings. For this reason, the minimum miscibility pressure of a gas injection is one of the main parameters of the flooding process.
When gas is injected into an oil reservoir, the flooding will be characterized by exchange of species between gas and liquid. Here, there are two possible processes. In one process, the most volatile species in the liquid evaporates into the injection gas and is carried away with it. Since the gas has substantially lower viscosity than the liquid, this yields a flooding that can lead to a very high recovery rate. In an alternative or simultaneous process, the heaviest species in the gas will condense into the liquid. Usually, both processes will take place at different places and times in a gas flooding.
Determination of minimum miscibility pressure for a gas injection has traditionally been done in laboratory experiments. In the experiments, oil is displaced with gas in a slim vertical tube filled with sand. The tube typically has a diameter of less than 1 cm and is several meters long. By measuring the recovery rate as function of the pressure, one can determine the pressure which gives almost full recovery (usually defined as a recovery factor of 0.97 after injection of 1.2 pore volume of gas). However, experiments show that minimum miscibility pressure is independent of the properties of the porous medium. The minimum miscibility pressure is thus a purely thermodynamic quantity, which can be determined through simulation with thermodynamic models. An important assumption then is that the selected thermodynamic description is representative for the current mixture of oil and injection gas. Also, note that by omitting porous media modelling, recovery rates and cell compositions obtained from simulations significantly below MMP will be less reliable as multiphase flow effects then become more important.
We compute the minimum miscibility pressure using a model that simulates one-dimensional gas injection without using flow equations for the flow [26,27]. Thereby, a purely thermodynamic determination of minimum miscibility pressure is obtained. For a more stringent approach, we refer to the tie-line method given by [28]. Our implementation contains the following steps: 1. Select a tentative pressure for recovery factor calculations. 2. Calculate equilibrium on a series of n cells after injection of a small gas volume. 3. Excess gas (and liquid) is transferred to the next cell in the series so that the volume of the cell is kept constant.
The transfer consists of gas if possible, however if not enough gas is present, the gas is supplemented with liquid so that the cell volume remains constant. See Fig. 2  4. Repeat the injection until 1.2 pore volume is injected into the cell series. 5. Calculate recovery factor as the volume ratio (at 1 atm, 288.15K) between the cumulative liquid production (after 1.2 pore volume CO 2 injected) and the total initial amount of liquid present in the cell-series. Cumulative liquid production is accumulated by logging the outflux from the last cell in a cell-series after each time step. By logging outflux also after the n/4 respectively the n/2 first cells of the n-cells-series, partial results from the computation for n cells provide results for n/4 and n/2 as if these series were computed independently. (If time Fig. 2 Illustration of a single cell in the one-dimensional simulation. The pressure and temperature is kept constant. The change of the equilibrium comes from change in the total composition due to added gas T corresponds to one pore volume injected for n cells, T /2 corresponds to one pore volume for n/2 cells etc.) 6. Increase the accuracy of the computed recovery factor by extrapolation to n = ∞ assuming the following asymptotic formula for the recovery factor: RF = a 0 + a 1 n −1/2 + a 2 n −1 + a 3 n −3/2 + · · · . In our implementation, rational extrapolation [29] is used, as is usual when solving ordinary differential equations, see also [30]. Rational extrapolation has the advantage that it is easily seen how many terms in the asymptotic formula there is a justification for using. 7. Repeat the above for other pressures. One should not select so high pressures that the recovery factor is above about 0.95. At such pressures, the recovery factor begins to flatten out, so that a good extrapolation is difficult. 8. Determine the pressure where the recovery factor is 0.97. This is done by logarithmic regression with the formula RF = β exp(αp), where p is the pressure. The miscibility pressure is then given as MMP = [ln(0.97/β)]/α. Figure 3 shows the determination of the minimum miscibility pressure (MMP) for the case with CO 2 injection in SPE5 oil. The minimum miscibility pressure is determined to MMP = 15.03 MPa.
To illustrate a process with injection of pure CO 2 in SPE5 oil above the miscibility pressure, we show the composition in a series of fifty cells, simulated at a pressure of 16 MPa. The composition is shown after injection of 0.12, 0.48, 0.84 and 1.2 pore volume, see Figs. 4-5.
The variation in density during a one-dimensional flooding is shown in Fig. 6. Here, the results of the flooding are shown in Fig. 4 (right).
Computing black-oil functions depending on the injected CO 2 , can be done in different ways. For injection of pure CO 2 , one can start from the oil composition and create new compositions in which CO 2 is added in appropriate quantities. For each such composition, one may then create black-oil functions. However, such functions are not always well suited for the simulation task since the functions do not capture the exchange of species which occurs between the liquid and gas phases. Results based on this approach is nevertheless included in the results section to illustrate under what conditions it can/cannot be applied. To capture this exchange it is therefore better to perform a onedimensional gas flooding as the MMP simulation described above (step 1-4) and then generate the black-oil functions for appropriate compositions in such a simulation.
We have chosen to pick out the compositions containing given CO 2 mole fractions at the time level after injection of 0.48 pore volume at a pressure of 16MPa, slightly above  MMP. This is shown in Fig. 4 (right). By interpolation we have selected compositions having CO 2 mole fractions 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95 and 0.99. Strictly speaking, the first of these compositions does not have zero CO 2 mole fraction, but has the CO 2 mole fraction in the unflooded part, i.e. about 0.000045. At 0.48 pore volume injected the cells closest to the injectors is almost completely filled with CO 2 while the cells furthest away are not jet reached. This makes it a suitable choice to pick compositions. Other choices around the same time will give similar results.
The black-oil functions of the oil phase in the reservoir are defined through Here v o,res is the molar volume of the oil phase in the reservoir, v l,o,sfc is the molar volume of the part of the oil Here v g,res is the molar volume of the gas phase in the reservoir, v v,g,sfc is the molar volume of the part of the gas phase which is gas (vapor) at surface conditions, v l,g,sfc is the molar volume of the part of the gas phase which is liquid (i.e. condensate) at surface conditions and λ g is the vapor fraction (mole fraction) of the gas phase in the reservoir at surface conditions. Figure 7 shows the functions B o and R s for the compositions mentioned above. Similarly, Fig. 8 shows the functions B g and R v for the same compositions. The surface conditions are T sfc = 288.15 K (15 • C) and p sfc = 0.101325 MPa (1 atm). The reservoir temperature is T res = 344.26 K. Figure 9 shows the density of the oil and the gas phase in the reservoir.
From the data set depicted in Fig. 7 undersaturated formation volume factor for the one phase liquid region are computed according to (Fig. 7, left). Here, p 40 = 40MPa is the largest pressure value tabulated. The saturation pressure p sat is computed from the primary variables (z, R s ) by inverse interpolation in the R s (z, p) table (Fig. 7, right). To obtain extra data points to provide stable extrapolation for R s beyond this table, we have also considered a leaner hydrocarbon mixture, where 20% mole fraction of light hydrocarbons has been combined with the base (as sampled from Fig. 4 (right)) hydrocarbon composition before the inclusion of CO 2 . Additional computational robustness is    The undersaturated formation volume factor for the onephase vapour region are considered to be of the form Here, (z, p, R v ) are the primary variables and B g (z, p) and R v (z, p) are sampled from data depicted in Fig. 8.
For viscosity of mixtures we use the so-called LBC model, [31] which is standard in most compositional simulators. For mixtures containing a large part of carbon dioxide, the modified LBC model is said to be significantly better than the LBC model [32]. To simplify the comparison with the compositional model we however stick to the standard model for viscosity. The viscosity is shown in Fig. 10.

Results
This section presents results from simulations where CO 2 dependent black-oil functions are used. The results are compared against compositional simulations and against a standard solvent model. Eclipse 300 [13] is used for the compositional simulations. The NOMIX and PSSTA option is used in Eclipse 300 to make the simulations more comparable. The NOMIX option avoids interpolation of the near critical relative permeability and the PSSTA option makes the detection of single phase regions more accurate in Eclipse 300. For the standard solvent model, OPM Flow [11] is used.
The test cases use fluid and grid data from SPE5 except that the injection gas consists of CO 2 . The reservoir is divided into 7 × 7 × 3 cells. An injector is located in [1,1,1] and a producer in [7,7,3]. For details of grid, and input data such as porosity, relative permeability, and permeability, we refer to [10]. The open data-sets including the PVT tables used in this study can be downloaded from opmpublications [33].
For the pressure dependent miscibility in the standard solvent model a linear ramp going from immiscible at 1500 PSI to fully miscible at 2300 PSI is used. The Todd-Longstaff parameter is kept at 0.6 for the standard solvent model. For the other models standard relative permeability functions without miscibility effects are used.
The SPE5 study consists of three different test cases. The original test cases are designed to cover different scenarios. The initial saturation pressure is around 2300 PSI. For the injection gas used in SPE5, the MMP is according to the SPE5 paper at around 3000-3200 PSI. For CO 2 the MMP is 15.3MPa, ie around 2200 PSI, see Fig. 3. To cover the same scenarios as in the original SPE5 case, we have therefore changed the setup slightly so that the pressure regime reflect the changed saturation pressure and MMP due to injection of CO 2 .
Case 1 starts with only production for two years so that the reservoir pressure drops below the saturation pressure. After two years, alternating WAG injection is started with a cycle of one year. The pressure now gradually rises above the saturation pressure and MMP. For Case 2, the goal is to exceed the saturation pressure during the entire simulation. Here injection and production start from the beginning. In Case 3, the pressure is first lowered by 5 years with only production and then stabilized near the MMP as the WAG injection starts. For Case 2 and 3, the WAG cycle is 3 months. All the simulations are for 20 years. The conditions used for the wells are listed for reference in Table 4.
In order to evaluate the different models, we select a cell [3,3,1] where we compare saturation and pressure. In addition, we compare production data for the different models. A comparison of cell values for saturation and pressure in Fig. 11 for Case 1 indicates that both models with CO 2 -dependent black-oil functions agree well with the compositional model. Since the pressure is below the saturation pressure and MMP, there will only be some exchange of the heavier components from the oil phase to the gas phase and the fixed composition version will give similar results as the dynamic composition version. The original solvent model gives greater deviations for the saturation. The visual deviations for the oil and gas saturation for the original solvent model is highly exaggerated due to the solvent model representing the injection gas with a separate saturation. The solvent saturation is shown in a separate plot in the same figure. Dissolved CO 2 will therefore not be part of the oil saturation directly for the standard solvent model Fig. 11 Comparison of cell values for the first case. For the standard solvent model the injected CO 2 is represented as a separate solvent phase and not part of the oil and gas phase but distributed between the phases using the miscibility function. Production data from Fig. 12 shows consistently larger deviations for the standard solvent model. This is especially clear for the gas production rate. The other models seem to agree more for the production data. Figures 13 and 14 show the results for Case 2. Looking at the oil saturation, we can see how the oil is gradually carried away as the gas saturation with a dominating CO 2 component increases. This is expected since we are well above MMP. If we compare the oil saturation we see that the standard solvent model as well as the model with CO 2 dependent black-oil functions with fixed base composition keeps some oil in the cell. After 5 years the cell enters a single phase region according to the composisional and the dynamic blackoil model. The solvent and the fixed model never forms a single phase oil region for this cell. The effect of this is most clearly observed in the water production where the solvent and the fixed model produce Fig. 12 Comparison of production rates for the first case. For the standard solvent model the injected CO 2 is represented as a separate solvent phase and not part of the oil and gas phase Looking at the oil saturation, we see how the oil saturation first gradually decreases during the depletion period of 5 years. When the water injection starts after 5 years it will push the oil to create the oil front we see in the Fig. 15. The oil is eventually washed out with the help of the CO 2 in the compositional model. The dynamic blackoil model is able to capture this effect, though two years later than the compositional model, the two other models are not.
The models seem to agree on the oil production rates, but as for Case 2 we see earlier water break-through for the compositional and the dynamic blackoil model. The gas production rate for the solvent model are significantly different from the other models. Note that the large oscillation in the gas production rate is caused by the WAG setup and is not a numerical artifact.

Conclusions
This article presents a new extended black-oil model that better represents the PVT properties of the mixture of hydrocarbons and CO 2 . In addition to the three common black-oil equations, an additional equation is added to represent the injection gas. Pre-processed black-oil tables that depend on the fraction of injection gas are used to represent the density and viscosity of the mixture. The hydrocarbon composition included in the calculation of the black-oil tables is sampled from a simulated slim-tube experiment. The exchange of hydrocarbon components that occurs when CO 2 is mixed into the oil is thus taken into account in the calculations of the density and viscosity of the phases. Comparison to the compositional model shows that the new extended black-oil model with CO 2dependent black-oil functions gives better consistency than a standard solvent model. Adding more data points in the Fig. 15 Comparison of cell values for the third case. For the standard solvent model the injected CO 2 is represented as a separate solvent phase and not part of the oil and gas phase table could further improve the accuracy of the proposed model. Still, a simplified black-oil type model will naturally not give the exact same answer as a full compositional model. However, by improving the accuracy of the density and viscosity calculations for the black-oil model, as well as by representing the important processes in the system such as component exchange, the goal is to create a model that provides accurate enough results while minimizing the simulation time.
Since the simulation time naturally scales with the number of equations and the non-linearity of the model the new extended black-oil model is expected to give longer simulation time compared to a standard black-oil model. Testing on a refined version of SPE5 indicates approximately a factor of 2 on the simulation time compared to a black-oil model where reservoir gas is injected instead of CO 2 . Of this increase in simulation time 70% is caused by the cost of the added equation, while the rest is due to Fig. 16 Comparison of production rates for the third case. For the standard solvent model the injected CO 2 is represented as a separate solvent phase and not part of the oil and gas phase increase in number of newton iterations. Testing on other cases indicates a estimate of a factor 2-3 on the simulation time when going from black-oil to the extended black-oil model. This is however significantly less than the expected increase in simulation time for a fully compositional model where the simulation time typically scales with the number of components in the model. The reduced simulation time compared to a fully compositional approach thus makes uncertainty and optimization studies needed for robust decisions more feasible.
The model developed herein can also be used for studies with injection gases other than CO 2 . Implementation in the open-source simulator OPM-Flow means that the model will be available for use and further development from the community.
Schlumberger for providing academic license for the Eclipse 300 compositional simulator.
Funding Open Access funding provided by NORCE Norwegian Research Centre AS.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.