On the injectivity estimation in foam EOR

In this work, we study injectivity issues caused by the use of the Peaceman equation in the numerical simulation of chemical enhanced oil recovery (EOR) processes aimed at reducing fluid mobility, such as foam injection, on coarse grids. Employing analytical solutions, we demonstrate that the Peaceman equation, commonly applied to mathematical modeling of the injectivity in commercial simulators, leads to errors of more than two orders of magnitude in the injection well pressure drop when the foam flow effects near the well are considered for assuming homogeneous mobility in a coarse-well block. To circumvent this issue, we investigate numerical treatments focused on the grid of well blocks through local grid partitioning strategies (radial and Cartesian) to improve the injection well bottom-hole pressure (BHP) estimation. This methodology does not change the input data nor the injectivity model characteristics of the commercial simulator. It does not significantly affect the computational cost of the simulation, since the grid treatment occurs only in the blocks containing the wells. Thus, the radial and Cartesian grid partitioning for the well block are compared using the STARS simulator. Our results show the clear capability of the methodology to reduce the well BHP overestimation, mitigating the errors caused by the Peaceman equation. Indeed, in some simulated scenarios, the BHP overestimation was reduced a hundredfold after applying the partitioning technique. Furthermore, we discuss the choice of simulation parameters leading to more accurate and numerically stable results.

The ratio between the first radial block radius and the well radius

Introduction
The injection of foam in porous media has been identified as a promising enhanced oil recovery technique (Gassara et al. 2017). During foam injection, better conformation control is expected than when no-foamed gas (either in a single slug or in a WAG scheme) is injected, because foam is shown to increase apparent viscosity of the injected fluid, leading to reduced gravity override and viscous fingering (Smith 1988;Hirasaki 1989a, b). However, when adopting the Peaceman well model (Peaceman 1978) [the most used in commercial software, such as STARS (Computer Modeling Group 2019)] in foam simulations, the well injectivity can be significantly underestimated when combined with a coarse grid (Rossen 2003(Rossen , 2013Leeftink et al. 2015;Li et al. 2017;Gong et al. 2019Gong et al. , 2020Soulat et al. 2020). The problem of injectivity loss in foam simulations occurs, because the Peaceman equation assumes uniform water saturation near the well. That issue can be exacerbated in numerical simulations due to the reservoir spatial discretization, in which inevitably the dimensions of any grid block containing a well are much larger than the wellbore radius itself, and the wells are mathematically represented as sources (injection well) or sinks (production well) within the coarse grid block (Peaceman 1978). This well modeling approach is usually not adequate when dealing with viscous fluids, such as heavy oil (Computer Modeling Group 2019), polymer solution (Li et al. 2017), or foam (Gong et al. 2020;Rossen 2013). The reason is that it cannot adequately capture near-wellbore effects that can have a significant influence on overall reservoir simulation results (Computer Modeling Group 2019).
Some techniques have been proposed in the literature to overcome the injectivity issues related to the foam EOR technique. Leeftink et al. (2015) present a methodology based on the computation of the correct well injectivity from the analytical solution for the water saturation around the well region using a radial flow model. However, their solution is over-restricted by simplifying hypotheses, and the methodology lacks a procedure applicable in commercial simulators. On the other hand, Soulat et al. (2020) focused on developing a methodology that could be embedded into commercial simulators. They propose a modification to the well index that takes into account only the foam quality. Nonetheless, the well index alteration is oversimplified and cannot be used when only gas flows through the well (e.g., during a SAG process). To reduce the overestimation of using the Peaceman model in the SAG foam process, Gong et al. (2019) proposed a negative constant for the skin factor based on a prediction obtained by a simplified radial flow model with fitted data from core flood experiments. Li et al. (2017) studied the difficulties in using the Peaceman model combined with an EOR by polymer flooding. In this case, since the Peaceman equation assumes an average velocity in the well block, the polymer apparent viscosity is incorrectly estimated, and so is the injectivity. The authors then propose a modification of the skin factor of the classical well model that can be readily input into STARS simulator.
In this context, this paper first investigates the issues encountered when the Peaceman model is used to compute the injectivity of foam using analytical solutions, as presented by Leeftink et al. (2015). This study demonstrates that the association between the assumptions of the Peaceman model with an excessively coarse-well block generates wrong injectivity estimates when the foam is used to reduce the mobility of the gas phase. Our goal is to present a simple tool that engineers can immediately use to circumvent the injectivity estimate issue. That is why, we focus on the implementation in STARS simulator. Thus, two types of grid partitioning are studied: Cartesian and hybrid. The Cartesian is characterized by increasing the number of Cartesian grid blocks around the wells (therefore, using smaller blocks). The hybrid grid method, readily available in STARS, is based on the definition of a cylindrical grid in the near-well region that can better capture the radial flow nature of the well injection/production flow. The hybrid grid technique depends on parameters associated with the number of grid partitions and the radius of the innermost cylinder. Notice that changing the grid partition using Peaceman model changes the physics of the model; it is not equivalent to 1 3 grid refinement in the sense of standard Numerical Analysis. One of the goals of this work is to obtain the proper choice of these parameters generating better estimates for the bottom-hole pressure (BHP), and also controlling numerical instabilities.
Differently from several other methods proposed in the literature (Leeftink et al. 2015;Soulat et al. 2020;Gong et al. 2019;Li et al. 2017), none of the well block partitioning approaches changes the input data or the commercial reservoir simulator injectivity model (including formulas or values for the well index or the skin factor, for instance). The results obtained with these strategies show a better estimation of the BHP (or equivalently, the well injectivity) as the number of well block partitions is increased, especially when the innermost block is closer to the wellbore. Moreover, the injectivity issues due to the use of the Peaceman equation are mitigated, and the deviations of the injectivity in real foam flooding applications from the simulations are greatly reduced. This leads to more reliable foam simulations that can be used to design field applications with more accuracy.
Local grid partitioning has been successfully applied to petroleum reservoir simulations for decades (Pedrosa Jr. and Aziz 1986;von Rosenberg 1982;Aziz 1993;Ding et al. 1998). This type of partitioning depends on several parameters, which are not directly determined. Up to our knowledge, a systematic parameter study is lacking in the context of the injectivity misestimation during the injection of high viscous fluids, such as foamed gas. We demonstrate the influence of the parameters (number of grid partitions and the radius of the innermost cylinder) of the hybrid grid technique in the estimation of BHP, emphasizing the values that present the best results and those resulting in numerical instabilities.
This work is organized as follows: The second section presents the mathematical formulation of the problem; the third section highlights the well injectivity issues through an analytical solution by the Method of Characteristics; next, in the forth section, we present grid treatments to the well block and numerical results using STARS; Finally, the fifth and sixth sections present discussions and conclusions.

Radial flow model including foam effects
We consider the problem of determining the pressure drop due to the injection of N 2 gas at a fixed volumetric rate Q through a well placed in the center of a reservoir block initially saturated with an aqueous solution of surfactant. For the matter of simplifying the presentation, let us assume both the rock and fluids are incompressible, the flow is isothermal, the gravity effects can be neglected, and, in the neighborhood of the injection well, the rock is homogeneous, the fluids are immiscible, and the oil phase is not present. Under these hypotheses, the flow is radial with the superficial velocity (Chen et al. 2006) where r is the radius, and H is the block thickness. To formulate the problem, we need some standard fractional flow theory relations (Chen et al. 2006;Buckley and Leverett 1942), presented next. The superficial velocity of the water phase is defined by where f w denotes the fractional flow of water (the ratio between the mobility of water phase, w , and the total fluid mobility, rt = w + g ): f w = w ∕ rt . The mobility of water ( w ) and gas ( g ) phases are defined by w = k rw ∕ w , g = k f rg ∕ g , where w and g are the viscosities of water and gas and k rw and k f rg are the Corey-type relative permeabilities of water and (possibly foamed) gas phases given by Notice that k 0 rw and k 0 rg represent the end-point relative permeabilities; S wc and S gr denote the connate water and residual gas saturations. The reduction of gas mobility due to foam injection (Computer Modeling Group 2019) assumes local equilibrium between the dynamics of creating and coalescing bubbles and follows an implicit texture model that can incorporate the effects of surfactant in water (F 1 ), the presence of oil (F 2 and F 5 ), the balance between viscous forces and surface tension forces (i.e., capillary number) (F 3 and F 4 ), the salinity of the brine (F 6 ), and the dry-out effects (FDRY). FMMOB defines the maximum gas mobility reduction. This work is focused on the suitability of the Peaceman equation to compute injectivity during foam injection, and therefore, we shall use a simplified foam model. Thus, we take F i = 1, i = 1, … , 6 . The only remaining function, FDRY, is represented by default in STARS (Computer Modeling Group 2019) as where SFBET indicates the abruptness of the dry-out effect, and SFDRY is a critical water saturation (below which foam collapses). The parameter SFDRY is considered constant in agreement with the hypotheses of this work: the MRF does not depend on surfactant concentration, oil saturation, salinity, or flow rate. These expressions are functions of the water saturation S w , which can be simplified by changing the variable to The remaining parameters are listed in Table 1. The relative permeability and foam parameters were defined in Kapetas et al. (2017) and Rossen and Boeije (2014) by adjusting the models to coreflood experiments on Bentheimer sandstone. A plot of the relative permeability curves (Eqs. (3) through (5)) is shown in Fig. 1. An illustration of the impact of foam on total mobility can be seen in Fig. 2.
The water saturation along the radius is the solution of the mass conservation equation (Chen et al. 2006) where is the porosity of the grid block containing the injection well. We define dimensionless variables x D and t D to describe cylindrical flow with radial symmetry (Leeftink et al. 2015) (7) where r w is the wellbore radius, and r b = h∕2 is the block radius. Substituting Eqs. (1), (2), (8), and (10) into Eq. (9)  yields To simulate the drainage of the well block, we define the left (injection) state J = S − = 0 and right (initial) state I = S + = 1.

Estimating injectivity in a grid block
Let us describe how to estimate the injectivity, why the wrong injectivity estimates appear in simulations, and finally, how to circumvent this issue.  Table 1 Parameter values for the foam flow in porous media based on data from Kapetas et al. (2017) and Rossen and Boeije (2014) for Bentheimer sandstone Peaceman injectivity in a grid block Following Peaceman (1978) and Chen and Zhang (2009), the difference between the pressure in a well and in its block can be computed as where P w is the well pressure, P re is the pressure at the equivalent radius (which will also be the pressure in the well block, because the pressure is assumed piecewise constant), and r eq is the equivalent radius, defined as The pressure drop estimate (Eq. 12) is widely used in reservoir simulations. As it is based on the average saturation in the well block (and therefore, on a uniform rt ), its application on coarse grids is limited to scenarios where the mobility of the fluid mixture does not change much near the well. Unfortunately, during foam injection, the total fluid mobility changes considerably (see Fig. 2), causing the pressure drop estimate given by Eq. (12) to be useless, as it often underestimates the injectivity of the well (Leeftink et al. 2015). An alternative to the Peaceman model was presented by Leeftink et al. (2015) providing a better injectivity estimate for the foam EOR processes. The idea is to drop the piecewise-constant hypothesis for fluid phases' saturations (at least in the well block) and compute the pressure drop by integrating the ODE from r w to r eq . The water saturation profile is given by the solution of Eq. (11) using the Method of Characteristics (see Appendix A) as plotted in Fig. 3. The natural limitation of this approach is that it requires an analytical solution, which is not available for compressible foam flow in porous medium.

How wrong injectivity estimates appear in numerical simulators
Let us compare the pressure drop computed using Eq. (14) with the one given by Eq. (12), where rt is dependent on the water saturation of a cylindrical well block, which is computed via the integration of the mass conservation equation by a combination of an implicit Euler scheme and Newton's method. Figure 3 shows these pressure drops and the one , computed using STARS software for compressible fluids using parameter values listed in Table 1. As one can see, the large spatial variation in water saturation (therefore, in the total mobility) causes the hypothesis of piecewise-constant mobility in the Peaceman model to be far from valid when foam alters the gas phase mobility. The pressure drop computed assuming uniform mobility increases until the water saturation in the well block reaches the SFDRY value, at which point the maximum pressure drop is attained. As the gas injection continues, the pressure drop rapidly decays. Although the physics of this problem requires a higher pressure drop initially to maintain the flow rate (due to the mobility reduction caused by foam) and a subsequent decrease in the pressure drop (due to the destruction of the bubbles caused by the water saturation reaching the SFDRY value), the analytical solution developed using the Method of Characteristics indicates that the Peaceman approach underestimates the injectivity most of the time, as can be seen in Fig. 3.
Notice that the simplifying hypothesis of incompressible phases (adopted to obtain the analytical solution) leads to a faster decrease in the water saturation of the well block, which in turn leads to an earlier increase in the pressure drop (after the injection of around 0.4 PV of gas) when compared to the pressure drop simulated in STARS (which reaches its maximum value after the injection of around 0.7 PV of gas). However, the pressure drop peak region and the overestimated behavior obtained by the Peaceman model in STARS remain.
In Fig. 4  phenomenon also observed in other works, as commented by Leeftink et al. (2015). As shown in Fig. 4a, the increase in the FMMOB parameter implies more significant differences between the pressure drops obtained using Peaceman and MOC, indicating that the error caused by the Peaceman equation is accentuated with the increase in the mobility reduction factor. On the other hand, the increase of SFDRY (presented in Fig. 4b) causes a reduction in the maximum pressure drop given by the Peaceman equation; this is due to the critical water saturation being reached earlier, causing the destruction of the bubbles and, consequently, an abrupt decrease in the pressure drop required to maintain the injection flow rate.

Circumventing injectivity issues in STARS simulator
As observed previously, the use of the classical Peaceman well model (Peaceman 1978) for reservoir simulation can lead to improper handling of the total mobility near the well and, therefore, resulting in an overestimation of the pressure drop (or equivalently the well BHP). This behavior is especially exacerbated when the well radius is much smaller than the well block itself, resulting in a poor estimate of the water saturation (and, therefore, of the phase mobility) in the well block, a region that presents high saturation variation (Leeftink et al. 2015;Gong et al. 2020). For an injection well located within a coarse grid block, if the BHP is fixed, the injectivity is excessively low, and if the injection flow rate is imposed, the BHP can get extremely high (Computer Modeling Group 2019). Similarly, for producer wells, the breakthrough time prediction cannot be adequately handled (Letkeman and Ridings 1970;MacDonald 1970;Computer Modeling Group 2019). A local Cartesian well block grid partitioning technique was proposed in Pedrosa Jr. and Aziz (1986) and von Rosenberg (1982) improving well modeling for two-phase flow. This method defines a grid partitioning of the well region only, avoiding an excessive global grid refinement and reducing computational costs. A schematics of the local Cartesian grid partitioning can be seen in Fig. 5a. The effect of the additional blocks on the computational costs, however, should not always be neglected. Also, a Cartesian grid partitioning does not match the radial geometry of the nearwell flow (Pedrosa Jr. and Aziz 1986).
Knowing that the accurate approximation of near-wellbore flow phenomena has a strong influence on the overall reservoir simulation, the hybrid grid concept was developed by Pedrosa Jr. and Aziz (1986) to match the need for an accurate and simple way to represent the two-phase flow in near-wellbore regions in reservoir simulators. In the hybrid grid approach, the entire reservoir is divided into well and reservoir regions, as shown in Fig. 5b: • Well regions A cylindrical grid system is used in the region around the well within the Cartesian well block, allowing for a more rigorous treatment of the near-well radial flow effects and mitigating the overestimation problem in the Peaceman well model. • Reservoir region Far away from the wells, the flow geometry may be considered to be approximately linear. Therefore, a Cartesian grid is appropriate for this region.
We follow both Cartesian and hybrid methodologies to study how to circumvent the injectivity misestimate for the foam EOR.
In STARS, the hybrid grid feature has several parameters; among them, we chose to analyze the following: the number of radial blocks, n rd , the ratio between the first radial block radius (denoted as r 1 ) and the well radius (Computer Modeling Group 2019) as shown in Fig. 5b. In the next section, we show how to improve the well injectivity estimate by adjusting these parameters. Notice that changing the grid partitioning for the foam flow simulations changes the physics through Peaceman equation. That is why, the effects of the grid refinement (see the next section) are not equivalent to the standard Numerical Analysis, where the physics is fixed.

Parameter investigation through numerical simulations
This section reports computational results based on a fivespot pattern obtained using the STARS simulator, version 2019.10. The production wells were located far enough apart, with a BHP of 1,000 kPa, so that no gas breakthrough would occur during the simulated period. The remaining (16) r 1 = i r w , parameters are the ones from Table 1. More details and the input file are in the supplementary material.
The Cartesian and hybrid strategies are used to produce two sequences of partitioned grids. The resulting injection well BHP for these grids is shown in Fig. 6. When a local (hybrid or Cartesian) grid partitioning strategy is applied, the average BHP required to maintain a fixed injection flow rate reduces significantly (Fig. 6a, b). Although even the most partitioned grids still present peaks in pressure (see the zoom in the peak regions in Fig. 6c, d), the maximum BHP clearly diminishes as the well block gets more partitioned, as shown in Fig. 7. Hence, it is expected that these peaks will vanish in an asymptotic successive grid partitioning.
In the hybrid grids, the radius of the first block is, by default, computed by STARS using relation (Computer Modeling Group 2019) In the most refined case ( n rd = 25 ), r 1 is approximately 0.781 m and, from Eq. (16), i ≈ 7.81.
We now show the results of a second study, in which we analyze the influence of i in the injectivity computation. In Fig. 8, the reference result ( i ≈ 7.81 ) is compared with i = 2.0, 3.25, 4.5, 5.75, 7.0 , which generates the values r 1 = 0.2 m, 0.325 m, 0.45 m, 0.575 m, 0.7 m for the radius of the first block. The results show that the well BHP peaks get lower and narrower as i value reduces, √ .

Fig. 5
The two block grid treatment strategies studied in this work. a well block region with a Cartesian grid partitioning; b well block region with a hybrid grid partitioning 1 3 i.e., as the grid partitioning becomes more concentrated near the well (see Fig. 9). The minimum value for the injection well BHP peak is 1,975 kPa, when n rd = 25 and i = 2.0.

Discussion
The approach studied in this work has advantages over the analytical methodology presented by Leeftink et al. (2015): (1) the analytical technique can compute only the well pressure drop, not being able to provide the well BHP (that is usually the desired quantity for well operation). Up to our knowledge, this limitation is due to the difficulties to obtain an analytical solution for the pressure equation; (2) STARS is based on a mathematical model that includes complex physical phenomena of multi-phase flow in porous media, such as compressibility, three-phase flow, and various EOR techniques; when the well block grid partitioning is applied in STARS, these phenomena are automatically taken into account. Analytical solutions are not available for these advanced model characteristics due to inherent difficulties that analytical methods, such as MOC, have when handling the equations. On the other hand, the hybrid grid partitioning technique can have practical limitations: • By manipulating the parameter , that controls the ratio between the outer radius and the inner radius for the remaining annular blocks (i.e., = r 2 ∕r 1 = r 3 ∕r 2 = … ), it is possible to increase the number of block grid partitions beyond n rd = 25 . However, the number of grid partitions has a limit. It is not possible to indiscriminately partition, because the computational cost increases, and numerical instabilities appear after a grid partitioning threshold. • Although the maximum value of the BHP decreases as r 1 → r w , numerical instabilities arise in the interval 1 < i < 2 . These instabilities increase the computational cost as STARS is forced to use smaller time steps in the simulation. The BHP peaks can also be affected by these instabilities, presenting higher values than when i ≥ 2 , which is not consistent with our expectations.  In the results presented in this work, the parameter combination with the best output was n rd = 25 , i = 2.0 and computed by STARS using the relation (Computer Modeling Group 2019) It is worth noticing that the outcome of using a specific parameter combination can strongly depend on the grid characteristics and the other simulation parameters. It can be inferred that the choice for moderate values of i and n rd , combined with the default calculation, leads to more promising results and reduced numerical instabilities. Even though the numerical experiments using well block grid partitioning still present high BHP peaks, these peaks are expected to vanish in an asymptotic behavior. However, these high peaks can be used to estimate the deviations of the simulations results from the reality of field applications.

Conclusions
• The analysis by the Method of Characteristics shows the pressure drop computed using the Peaceman equation (implemented in STARS) can greatly misestimate the injectivity during foam injection. The leading causes of this issue are the simplifying Peaceman model assumptions, such as a uniform water saturation in the well block, associated with a coarse grid. • The well block grid treatment adopted in STARS circumvents the injectivity problem during foam injection in porous media; • Differently from other methodologies, the employed block well partitioning technique (Cartesian or hybrid) does not change the input data nor the characteristics of (18) = h 2 r 2 1 1 2(n rd − 1) .
the injectivity model of the STARS simulator. It applies a local grid partitioning only in the blocks that contain wells, without significant computational impact. • Our results show that the well block grid partitioning, either Cartesian or hybrid, significantly improves the estimation of the injection well BHP needed to maintain a constant injection flow rate. • The hybrid grid partitioning can better represent the radial flow characteristics near the well region and requires fewer grid partitions when compared to the Cartesian grid. • The increase in well block grid partitioning, following both strategies, significantly reduces the BHP peaks. This reduction is especially accentuated with the hybrid grid associated with a thin first radial block. • The proper choice of parameters for the hybrid grid technique can improve the estimation of BHP. On the other hand, inadequate choices can increase computational costs and generate inaccurate estimates due to numerical instabilities.   where S M is obtained through and is the shock front's velocity defined by Rankine-Hugoniot equation The solution (Eq. 19) has a simple interpretation in the context of the method of characteristics as a rarefaction wave followed by a shock wave. Here, rarefaction wave is the spreading continuous part of the solution and the shock wave is the discontinuous one; for more details, see Smoller (1994). According to the theory, the shock curve must satisfy Rankine-Hugoniot equation (Eq. 21) and a tangent condition graphically represented in Fig. 10. Riemann problem solution profile and the corresponding characteristic curves are plotted in Fig. 11a, b.