An improved procedure for generating pseudorelative permeabilities for water flooding in stratified reservoirs

This paper presents an improved procedure for generating pseudorelative permeability curves for stratified water flooding using either constant pressures or constant flux at the reservoir/grid block boundaries. The concept of pseudorelative permeability reflects the generation of a relative permeability curve that can be used to represent the entire reservoir thickness, rather than a specific layer during reservoir simulation, thus saving computational time. In this paper, fractional flow theory is applied to the generation of pseudorelative permeability curves for (1) constant flow rate and (2) constant pressure boundary conditions. Previously, pseudorelative permeability curves were generated for constant flow rate only, since the analytical solutions for constant pressure boundaries were non-existant. In this paper, this restriction has been removed based on novel analytical solutions for constant pressure boundaries. The method within this paper also differs from previous methodologies and studies, which are all based on an approximation using piston-like displacement for water flooding. Instead, this new model uses fractional flow theory to its fullest extent to generate pseudorelative permeability curves that are physically more realistic. The solution is extended to generate pseudorelative permeability curves for waterflooding of a reservoir under the assumption of constant pressure boundaries which is an equally realistic assumption in comparison to constant flow rate. The generated pseudorelative permeability curves are used in a 2D areal reservoir model in a standard reservoir simulator to predict the behavior of the fully layered 3D reservoir model. It is found that there is considerably better agreement between the results obtained with this new method and the fully layered reservoir model compared to previous methods.


Introduction
Reservoir simulation is a powerful tool to study and predict oil recovery for various drive mechanisms. The results of such studies are used to make capital management decisions for field developments and future reservoir operating strategies.
Three-dimensional multi-phase fluid flow simulation through a complex reservoir with high grid block resolution can require substantial human effort and computer resources. Reducing the number of grid blocks, i.e., making the model coarser, is one option to reduce computational time; however, this approach may ignore important finer scale reservoir details, resulting in discrepancies between simulation results and actual production history.
Generally speaking, an order of magnitude of 100,000 grid blocks can be used in reservoir simulation (some reservoir models may use more grid blocks as several million grid blocks), typically 100 × 100 × 10 in the two horizontal and one vertical directions, respectively. This grid block volume potentially exhibits a large degree of heterogeneity and if a core flood measured relative permeability curve is used to simulate such a large grid block, the important effects of scale dependency on the relative permeability curves and also on reservoir heterogeneity are neglected, thus significant errors may occur in simulation results (Cao and Aziz 1999).
Relative permeability curves are the main parameters that govern simultaneous multi-phase flow in reservoir simulation. They are normally obtained from steady-or unsteadystate two-phase core floods. To consider the heterogeneity within a large grid block, pseudo(up-scaled)relative permeability curves can replace intrinsic relative permeability curves (Cao and Aziz 1999). Pseudorelative permeability curves attempt to capture the effects of scale and heterogeneity within the coarse grid blocks. The goal of up-scaling is to reduce the CPU time needed for calculation either "by reducing the number of grid blocks or reducing the number of dimensions of the problem, such as reducing a 3D field case model to a 2D areal model" (quote from Cao and Aziz 1999). In other words, through up-scaling, we hope to retain precise information in less computational time and effort (Dykstra and Parsons 1950;Coats et al. 1967Coats et al. , 1971Hearn 1971).
The vertical characterization of a reservoir (Tompang and Kelkar 1988) is one of the most important considerations in reservoir simulation. The most accurate way to consider vertical effects is to use a 3D reservoir simulator; however, due to the computational constraints discussed above for 3D simulation, the use of 2D simulators is an attractive alternative. As pointed out by Hearn (1971), 2D reservoir simulation implies uniform reservoir properties and fluid saturation throughout the reservoir thickness, an assumption that surely is not physically correct. Thus, the input data require tuning to approximate the vertical effects (Hearn 1971). Dynamic pseudofunctions such as (Jacks et al. 1973, Kyte and Berry, 1975and Stone 1991, flux weighted and pore volume weighted could be used to upscale relative permeability; however, these methods need fine grid simulation results. History matching using and tuning relative permeability curve using genetic algorithm is another way to approximate the vertical effects. However, this approach needs many simulation runs to fine-tune the relative permeability (Fayazi et al. 2016).
Previous authors (Coats et al. 1967(Coats et al. , 1971) have shown the validity and advantages of two-dimensional areal simulation through the use of appropriate pseudorelative permeability and capillary pressure functions when vertical equilibrium exists. This assumption is applicable only for highly communicating layers with low displacement velocity.
An alternative method is needed to calculate up-scaled pseudorelative permeability curves in reservoirs without vertical equilibrium, in which case the displacement process is dominated by viscous forces rather than by gravity forces. Such a method was presented by Hearn (1971), which forms the basis for this paper. In reality, the use of an areal model to represent the stratified system lies between the two extreme cases represented by vertical equilibrium (Coats et al. 1967(Coats et al. , 1971) and non-communicating layers (Hearn 1971).

The method by Hearn/Dykstra-Parsons
We now briefly discuss the method proposed by Dykstra and Parsons (1950) and Hearn (1971) to determine pseudorelative permeability curves to approximate viscous flow in a stratified reservoir under the assumption of a constant flow rate in each layer. This assumption is an approximation at best, since this assumption can never hold true even if the total flow rate over all layers is constant, unless the mobility ratio is unity. This is because the total mobility in a layer, in general, varies as the fluids propagate and saturations change from injection to production within the layer.
The stratified model used in this paper representing a simulation grid block is shown in Fig. 1.
The fractional flow function for each layer is given by where k rwi and k roi can, for example, be calculated using a Corey model, i.e. (1) We assume horizontal flow for the sake of simplicity in Eq. 1. However, this assumption can be relaxed by including a gravity term. We use the Corey model to model relative permeability, as typical for water wet rocks, with a wi = 0.3 , a oi = 0.8 and n wi = n oi = 2 for all layers i for simplicity.
In the Hearn/Dykstra-Parsons method, the layers are first sorted according to decreasing water breakthrough time to simplify the computational procedure. In the Hearn method, the ordering n = 1, … , N is determined using an iterative procedure and in Dykstra-Parsons method the layers are ordered in decreasing values of the ratio where * wi is endpoint water mobility, M i is endpoint mobility ratio and ΔS i is mobile water saturation and are calculated as equations Pseudorelative permeabilities for water and oil, respectively, are defined through linear interpolation between the points S wn ,k rwn N n=0 ; S wn ,k ron N n=0 , where subscript 0, n refers to before breakthrough of the first layer and after breakthrough of layer i, respectively. Here, S w refers to pore volume averaged water saturation at the outlet face of the grid block, and k r refers to the phase mobility averaged pseudorelative permeability. Note that since a piston-like (4) displacement is assumed, saturation values used are only layer connate water and residual oil saturations. These points are calculated as follows: The above procedure is based on the assumption of piston-like displacement in each layer, which ignores the smooth variations in saturations after breakthrough in a layer. In other words, the water saturation behind and ahead of the front in layer i is 1 − S ori and S wci , respectively. Equation (8) gives the pseudowater saturations for all layers as the average water saturation at connate water (wc), water saturation from connate water to residual oil (wn) and the water saturation at residual oil (wor) conditions. These saturations are then used to calculate the pseudorelative permeabilities k rw ,k ro for water and oil in Eqs. (9) and (10). The objectives of this paper are twofold. First, we modify the Hearn/Dykstra-Parsons method using fractional flow theory instead of the assumption of piston-like displacement. We do this by first presenting a generalized method for developing pseudorelative permeability curves using the assumption of constant flow rate in each layer, as in Hearn (1971). Second, the main result of this paper is presented which is to extend the generation of the pseudorelative permeabilities to constant pressure grid block boundaries. In this, no approximation is made since a closed formula solution for the flow rate as a function of time was presented in Johansen and James (2016), see also "Appendix A". The idea of generating pseudorelative permeability curves to simulate waterflooding in a reservoir with constant grid block pressure boundaries has not been previously studied, and as will (8) Oil & water production Water injection Fig. 1 Stratified reservoir model be demonstrated, it is a more accurate method compared to the constant flow rate-generated pseudofunctions.
The pseudorelative permeability curves generated under the constant pressure boundary condition is an equally realistic well-operating condition as constant flow rate. For example, unless the reservoir pressure is close to bubble point initially, the wellbore pressure usually will be kept constant above the bubble point. Then, injectors are also conveniently operated at constant wellbore pressure, and the pseudofunctions of this paper can also easily be used in a well-to-well calculation using, e.g., a 2D stream tube method.

Calculation of pseudorelative permeability curves using the new method
The stratified model is depicted in Fig. 1. The reservoir is divided into layers based on absolute permeabilities (from core data), with each layer characterized by a thickness (h), absolute permeability (K), porosity ( ) , and relative permeability curves, in particular, residual oil saturation (S or ) and connate water saturation (S wc ). Assumptions of this reservoir model are: • Oil and water are incompressible fluids. • Capillary forces are negligible compared to viscous forces. • Reservoir properties are homogeneous within each given layer. • Layers are non-communicating.
It should be noted that the work flow is the same for the constant pressure boundaries and the constant flow rate cases with the exception of the calculation of the exact flow velocity in each layer as a function of time for the constant pressure boundary case. Consequently, the breakthrough time for each layer is exact for the constant pressure boundary case (Shadadeh 2015).

Calculation of breakthrough time
For the case of constant flow rate in each layer, the fractional flow theory (Buckley and Leverett 1942) is used to calculate an approximate breakthrough time and average water saturation at the outlet after water breakthrough of each layer (see "Appendix A"). As already mentioned, this is at best an approximation since the velocity in each layer will vary as the flow progresses since the constant flow rate on the inlet face of the grid block will distribute differently at different times among the layers according to the changing mobilities in the layers.
Unlike the case using constant total flow rate boundary condition in the grid block, in the constant pressure boundary case the total Darcy velocity u T = q T ∕ A varies over time and consequently, the calculation of layer breakthrough times and locations of any given saturation value is different from the constant rate case. Since the outlet saturations and breakthrough times are different in the two cases, the pseudorelative permeability curves are also different. As we also will show, the constant pressure boundary generated pseudofunctions yield more accurate results than the Hearn/ Dykstra-Parsons pseudofunctions (Shadadeh 2015).

Calculation of pseudowater saturations
The exact breakthrough times are calculated according to Johansen and James (2016), "Appendix A". The layers are sorted in descending breakthrough times and the pore volume averaged water saturation in the reservoir is calculated before the first layer breaks through (n = 0) and immediately after a new layer breaks through at the outlet (n = 1, … , N) according to where S * wi is water saturation at the outlet face of each layer at the moment of water breakthrough, which is calculated by the Welge tangent procedure (1952) in Figure B.1. Using only the breakthrough times when calculating pseudosaturations may yield too few points on the pseudocurves, in particular, after all the layers have broken through (corresponding to high water saturations). However, any sequence of times t 0 = 0 < t 1 < ⋯ < t M can be selected and the corresponding layer saturations at the outlet be calculated according to Johansen and James (2016), "Appendix A". The pseudowater saturations are then calculated using Eq. (11), where in this equation S * wi is the water saturation in layer i at the outlet at time t j , j = 0, ..., M.

Calculating pseudorelative permeabilities
The end point pseudorelative permeabilities for oil and water are calculated as in "Introduction", i.e., k rw S wc = 0,k ro S wor = 0 . The pseudorelative permeability for water at residual oil saturation, k rw S wor , and the pseudorelative permeability for oil at connate water saturation, k ro (S wc ) , are calculated from Eqs. 9 and 10. The pseudorelative permeabilities corresponding to the sequence S * wi , generated by calculating the average water saturation at the (11) outlet face of the grid block at the breakthrough times (or the selected times), are calculated using

Example calculations
A hypothetical reservoir with 20 distinct layers was used to calculate the pseudorelative permeability curves for (1) the constant flow rate case and (2) the constant pressure boundary case. The characteristics of the reservoir are shown in Table 1: In this example, the relative permeability of each layer is assumed derived using Corey model as below: It is assume that for this example Corey water and oil saturation exponent (a w , a o ) are equal to 2, a w = 0.3 (water rel. perm end point) and a 0 = 0.8 (oil rel. perm. end point) for all layers. Despite the same endpoints and same saturation exponents, the shape of relative permeabilities for each layer is different, due to different connate water saturation and different residual oil saturation. It should be mentioned that the aforementioned exponent and coefficients could be different and for simplicity, in this example, they assumed to be same for all layers. The relative permeabilities for each layer are shown in Fig. 2.
The results of generating the pseudorelative permeability (for all layers) by the two methods are shown in Fig. 3. This figure shows the pseudorelative permeability of water and oil versus average water saturation for both the constant  Fig. 3 is that they are different for the constant pressure boundaries and constant flow rate cases. Point A on Fig. 3 represents the last breakthrough of a layer, at this point breakthrough at all the layers have occurred. After this point all layers have a positive contribution to pseudowater relative permeability. Therefore, the slope of pseudorelative permeability of water increases and pseudorelative permeability of water reaches a maximum value of 0.3 as the coefficient of the Corey model for the water relative permeability, and the pseudorelative permeability of oil gently decreases to zero. Point S w = 0.744 . represents the infinite time when all layers produce water with outlet water saturation of S w,out i = 1 − S ro i . Point A in Fig. 3 represents the last breakthrough of a layer, at this point breakthrough at all the layers has occurred. After this point all layers have a positive contribution to the water pseudorelative permeability and the slope of the water pseudorelative permeability increases reaching a maximum value of 0.3 and the pseudorelative permeability of oil decreases to zero.

Results and discussion
Based on the above example (Table 1), two reservoir models were constructed in Schlumberger's Eclipse 100 black oil reservoir simulator. Onis a 3D cross-sectional fully layered model with 100 (5 m) * 1 (50 m) * 20 (2 m) grid blocks and each layer has its own properties as shown in Table 1, which is considered as the base case. Other PVT and Rock data, initial condition and datum depth of this model are shown in Table 2.
The second model is an areal pseudoized model with 100 (5 m) * 1 (50 m) * 1 (40 m) grid blocks, using the average porosity, thickness weighted averaged permeability and pseudorelative permeabilities generated by the four pseudoization methods discussed in "The method by Hearn/Dykstra-Parsons". These four methods include a new proposed method and Hearn method for constant flow rate boundary condition and another new method and Dykstra-Parsons method for constant pressure boundaries. The information in Table 2 is also used for this model as well.
The pseudorelative permeabilities using Hearn, Dykstra-Parsons, the new method with constant pressure boundaries, and the new method with constant flow rate are shown in Fig. 3.
Both 3D and 2D models were used for two operating conditions of constant pressure boundary and constant flow rate. The constraints for the two cases are different: the constant pressure boundary case uses BHP target of 50 bara for producing well and 150 bara for injection well, while the constant flow rate case has a reservoir producing flow rate target of 10 rm 3 /D, which by considering oil formation volume factor would be 7 Sm 3 /D.
Permeability distribution for 3D and 2D reservoir model for the discussed example is shown in Figs. 4 and 5, respectively.
After running the models, the results are compared in Figs. 6,7,8,9,10 and 11. Figures 6,7 and 8 are for the constant flow rate case and Figs. 9, 10 and 11 are for the constant pressure boundary case. Figure 6 shows that for constant flow rate case, the Hearn/ Dykstra-Parsons method results in a smooth decline in the oil production flow rate. However, the base case model shows a sharp decline in oil flow rate at year 13. The new method for generating pseudorelative permeability curves captures this trend. Figure 7 shows the same in terms of recovery factor as a function of time. Figure 8 shows the same in terms of water cut versus time (for the constant flow rate case). As shown, using the fractional flow theory to its fullest extent is superior to the use of piston-like displacements for the cases where constant flow rates in grid blocks are used to calculate pseudorelative permeabilities.  Figure 9 shows the oil flow rate for the constant pressure boundary case. From this figure, it is clear that the new method is superior to the Hearn/Dykstra-Parsons method (which is calculated based on the incorrect assumption of constant rate). Figure 10 illustrates recovery factor versus time for the constant pressure boundary case. The Hearn/Dykstra-Parsons method strongly over predicts recovery compared to the base case and the new method for generating pseudofunctions. The latter method is slightly underpredictive, particularly after water breakthrough. Figure 11 shows water cut versus time. The Dykstra-Parsons shows an early water breakthrough at year 6 and then a smooth increase in water cut. However, the base case has a sharp increase at year 13 and then stabilizes at a high value. The new method successfully captures this trend.
After year 17, all the methods have ± 2% difference in water cut with the base case (fully layered).
From these results, it can be concluded that using the new method for generating pseudorelative permeability curves for both the constant pressure boundaries and the constant flow rate cases provides results in good agreement with the results of the 3D fully layered simulation, which uses a higher number of grid blocks and is considered as the basis of comparison.
The introduced method is semi-analytical, i.e., analytically, based on fractional flow theory, pseudorelative permeabilities are generated and then are used in a reservoir model and is solved numerically using Eclipse reservoir simulator or any other commercial reservoir simulator. It should be mentioned that this new approach unlike dynamic pseudofunctions does not require fine grid simulation results.
The limitation of the presented new approach includes: 1. Supplementary program is required to implement fractional flow theory and generate pseudorelative permeability curves. 2. The model assumes no crossflow between layers, this situation in reality may occur in two cases, structurally non-communicating reservoir layers or viscous forces being dominant compared to capillary and gravity forces. As much as this assumption is achieved, the results would be more accurate. It should be mentioned that for the example mentioned in the paper, however, the reservoir layers are not structurally non-communi- cating (the vertical permeability of each layer is assumed to be half of the areal permeability of the layer), still the results achieved using the new approach are in a good agreement with full layer reservoir modeling. 3. The approach is aimed to capture the heterogeneity effects in vertical direction, and it assumes each layer has homogeneous properties within it.
A saturation position in this case can be calculated by Eq. A-1 based on fractional flow theory (Buckley and Leveret 1942): Consequently, by choosing x si = L the breakthrough time (t BT ) is calculated by the following equation: where A is the flow area, is porosity, q T is total flow rate and f ′ wi is the derivative of fractional flow at saturation S w . For t < t BT , the outlet saturation will be S wc . However, for t > t BT, we can calculate x s for some point in the interval S * w , 1 − S or where S * w is the front saturation, and obtain the saturation at the outlet face by trial and error, or a numerical method. The saturation corresponding to x s = L is considered as outlet saturation. The leading shock saturation S * w is calculated as Figure B.1 based on the Welge (1950) method.

Calculation of the water saturation at the outlet for constant pressure boundaries
A saturation position for the constant pressure boundary case is calculated by the following equation based on Johansen and James (2016): where u T (t) is the total flux velocity at time t and is given by The breakthrough time is calculated as follows: To compare water front profile of new approach with Buckley-Leverett method, a series of water saturation profile before breakthrough and after breakthrough for the two cases of constant flow rate and constant pressure boundary case was generated. The result is shown through Figs. 13a, b and 14a, b.
The results are reasonable, since using Buckley-Leverett will assume a 1D homogenous model and the saturation front will be sharp; however, the studied 3D model has vertical heterogeneity though the saturation front will smear out and our new approach despite up-scaling in vertical direction was successful to capture the heterogeneity effect in vertical direction.