Analytical solution to pressure transient analysis in non-Newtonian/Newtonian fluids in naturally fractured composite reservoirs

During polymer flooding into the reservoir using shear-thinning non-Newtonian fluids, an interface exists between the in situ Newtonian crude oil in the reservoir and the injected polymer solution. This paper examines the application of non-Newtonian well test analysis techniques to develop an analytical solution to non-Newtonian/Newtonian fluids with different flow index composite reservoirs for estimating the interface boundary conditions in such reservoirs. Mathematical models were presented and solved analytically using the Laplace transform. The solution was inverted to real domain with the Stehfest algorithm (Stehfest Commun ACM, 1970). Solutions were obtained for three boundary conditions that are infinite acting, constant pressure and no-flow boundaries. New pressure and pressure derivative type curves are developed for naturally fractured reservoirs. A general solution was obtained, which is adaptable to the case of non-Newtonian/Newtonian fluid interface. Numerical examples are used to estimate the radius of investigation, reservoir and power-law flow parameters using type curve matching and Tiab’s direct synthesis technique.


List of symbols B
Formation volume factor, rm 3 /stm 3 C D Dimensionless wellbore storage c t Total compressibility Pa −1 h Formation thickness, m H Consistency index (power-law parameter) Pa s n k Formation permeability, m 2 m, n Flow behavior index (power-law parameter), dimensionless P D ′ Dimensionless wellbore pressure derivative P D Dimensionless wellbore pressure drop P i Initial pressure, Pa P wD Dimensionless wellbore flowing pressure P D Dimensionless pressure drop in Laplace space q Volumetric flow rate, m 3 /s r Radial distance, m r D Dimensionless radial distance at any point, (r/r w ) R eD Dimensionless radial distance at boundary conditions (external)

Introduction
Enhanced oil recovery (EOR) encompasses processes of improving the recovery of hydrocarbons from oil reservoirs after primary production has declined. Among chemical-based enhanced oil recovery processes, polymer flooding is a highly practiced process in conventional oil reservoirs. Polymer flooding utilizes the injection of polymer solutions into oil reservoirs (Manrique et al. 2017).
The presence of a polymer in water increases the viscosity of the injected fluid, which upon injection reduces the water-to-oil mobility ratio and the permeability of the porous media, thereby improving oil recovery. The net result is the increase in oil displacement, reservoir sweep efficiencies, and pressure gradient, especially in heterogeneous oil reservoirs (Manzoor 2020). Well test analysis is often carried out before, during, and after polymer flooding to obtain an accurate reservoir description and to evaluate fluid mobilities within the project area (Huh and Snow 1985). Many field applications of polymer flooding have been reported as technically and economically successful. However, Non-Newtonian power-law behavior of polymer solutions such as lack of a simple proportionality relation between shear stress and rate of shear, non-uniform polymer distribution within the reservoir, and the decrease of apparent viscosity with increasing shear rates (Tariq Fariss Al-Fariss 1984) has led to various setbacks during polymer flooding. As a result, different works have been done to investigate how to maximize the oil recovery during polymer flooding. Manzoor (2020) studied the effect of injection pressure and polymer concentration on flooding performance.
During polymer flooding, the reservoir is often made up of two or more regions and, as such, is regarded as a composite reservoir. Most times, the two regions, which might be oil and polymer solution/water regions, have different rock and fluid properties. Ikoku and Ramey (1979) derived a linearized partial differential equation for describing the radial flow of a slightly compressible non-Newtonian power-law fluid through a homogeneous porous medium. They obtained approximate analytical solutions to the equation by considering constant rate injection in an infinitely large radial flow system. The solutions provided a new method of well test analysis for non-Newtonian fluids injected into reservoirs. However, they did not consider the discontinuity of fluid in the reservoir at the interface. Lund and Ikoku (1981) addressed this neglect through numerical method by extending the flow of non-Newtonian fluids through a porous media onto composite reservoirs containing non-Newtonian and Newtonian fluids. Ertekin et al. (1987) studied the pressure transient behavior of non-Newtonian/Newtonian composite fluid flow in a porous media. Their work employed numerical method and primarily focused on a porous media with a finite-conductivity vertical fracture. Gencer and Ikoku (1984) analyzed injectivity and falloff test data during two-phase flow of non-Newtonian (power-law) and Newtonian fluids to determine the saturation gradients that can develop. They considered simultaneous (multiphase) flow of oil and water, and effects of relative permeability characteristics of the porous media, finite radius of the skin region, skin factor, injection time, viscosity of the Newtonian fluid, and rock compressibility on the pressure response.
Katime-Meindl and Tiab (2001) applied Tiab's direct synthesis (TDS) technique to the case of flow of non-Newtonian fluid through porous media. Pressure and pressure derivative data of a well test data in a homogeneous reservoir are interpreted without the use of type curve matching or regression analysis. The analysis considered the effect of skin factor and wellbore storage, and equations for calculating mobility and radius from wellbore to the nearest boundary were obtained. Escobar et al. (2010) utilized the pressure derivative curve in interpreting well test data in reservoirs with non-Newtonian power-law fluids. They obtained equations to estimate the permeability, skin factor and non-Newtonian bank radius. The characteristic features found on the pressure and pressure derivative curves enabled the application of Tiab's direct synthesis (TDS) technique to their work. No analytical mathematical model developed to describe the pressure behavior in composite reservoir with non-Newtonian/Newtonian fluids.
In this study, we developed a general analytical solution to describe pressure transient behavior in a reservoir with two immiscible different non-Newtonian fluids having different flow indexes. Figure 1 represents the physical description of the composite reservoir. TDS technique is used for the analysis.

Mathematical model
The linearized partial differential equation for naturally fractured reservoir (NFR) as presented by Omosebi and Igbokoyi (2015) is: All the dimensionless variable and derivation are in "Appendix". (1) (2) f (s) = (1 − )s + (1 − )s + The assumptions used in developing this model are: (1) fluid flow in the reservoir is radial, (2) laminar and Darcy's flow is valid, (3) isothermal, single phase and slightly compressible fluid with constant properties, (4) homogenous and isotropic system, (5) pseudoplastic (shear-thinning) fluid obeys Ostwald-de Waele power-law relationship and (6) steady-state effective viscosity. The mathematical model also assumed there is no transition region.
The approximate long-time solution at the wellbore when r D is unity is: The solutions to Eq. 1 in Laplace space with various boundary conditions are presented in "Appendix". Solution with the dimensionless wellbore storage C D and skin S is:

Type curve
Apart from the flow index, interporosity transfer and storativity ratio, two key parameters are used to characterize the type curve. They are alpha and theta (see Eqs. 38 and 59 in "Appendix"). and While alpha is a function of diffusivity ratio D, theta is a function of mobility ratio M. Therefore, The effect of mobility and diffusivity ratios is demonstrated in Figs. 2 and 3 as in the work of Ambastha (1988). Figure 4 is the curve generated for the three outer boundary conditions. Figures 5 and 6 show the condition when the inner region is Newtonian and the outer region is non-Newtonian. This is to replicate a case when the pressure measurement is taking place at the producer of the Newtonian fluids; for example, crude oil that exhibits Newtonian behavior. The early pressure derivative is 0.5, while the late pressure derivative exhibits non-Newtonian behavior. Figure 7 is a reversed behavior of Fig. 5. However, the non-Newtonian flow regime in the inner region has a strong influence in the outer region with Newtonian flow. Odeh and Yang (1979) derived the expression for the radius of investigation in non-Newtonian fluid flow in porous media. This expression is simplified to: Fig. 2 Effect of mobility ratio: ω = 0.01, λ = 10 -5 , n = 1.0 (inner region) and m = 1.0 (outer region), infinite system

Fig. 3
Effect of diffusivity ratio: ω = 0.01, λ = 10 -5 , n = 1.0 (inner region) and m = 1.0 (outer region), infinite system Ikoku and Ramey (1979) also used a similar approach to arrive at their formula for radius of investigation as: Equations 14 and 15 are used in this study to compute the polymer front. Equation 14 gives radius of investigation in centimeter, while Eq. 15 gives it in meter. Ikoku and Ramey (1979) also presented the formula for estimating the effective viscosity of the polymer fluid as: where the effective viscosity is not given, Eq. 16 is used to estimate the effective viscosity.
The main purpose of this work is to get a general solution for a composite reservoir with two regions having  Fig. 8. This type of solution is particularly useful for heavy crude that exhibits non-Newtonian behavior and being displaced by polymer fluid for increased recovery.

Applications to numerical examples
The two numerical examples used to verify the mathematical models are from Lund and Ikoku's (1981) work.
The numerical examples are for homogeneous system with storativity ratio ω equal to 1.0.

Injection test
The match in Fig. 9 was obtained with a flow index n of 0.9 (Table 1).
We estimated the effective viscosity with Eq. 16. The value of the effective viscosity with n of 0.6 is 0.0000402 Pas, while n equals 0.9 giving 0.00424 Pas. The effective viscosity of the polymer cannot be less than that of the Newtonian fluid in the formation; hence, we have confidence in the match with 0.9 as the value of n, and therefore, the value of the polymer effective viscosity used in this work is 0.00424 Pas.
We used pressure derivative plot to estimate the time at the interface. Lund and Ikoku (1981) used the plot of ΔP against t v to estimate the time at the interface. The pressure derivatives give better estimate of point of discontinuity. From Fig. 10, we estimated the interface to be at the time  (Table 2).
The alpha (α) value for both match is 0.003 and theta (θ) is 0.00135.
The effective compressibility in the polymer zone, which is region I, is estimated from the value of alpha and time match.
From the time match: From the value of alpha:

3
Falloff test from the same numerical model.
The Agarwal equivalent time is used for all the plots in the falloff test (Figs. 11, 12).
The same values of alpha and theta were used (Table 3).

3
From the time match: From the value of alpha:

Tiab's direct synthesis method
Interporosity flow parameter is obtained by substituting the co-ordinates of the minimum point of the trough, t min and (t x ΔP′) min according to Tiab et al (2006). At that point, For storage capacity ratio, the method of calculation is obtained using the pressure derivative co-ordinate of the minimum point and the radial flow regime line, according to Tiab et al (2006). It is given as: At the infinite-acting radial flow regime of non-Newtonian fluid, the log-log plot of the pressure derivative versus time will give a straight line with slope, v. The value of n can then be determined from the slope. By substituting the dimensionless pressure equations into the pressure derivative equation, the mobility ratio (effective permeability-viscosity ratio) is calculated by Igbokoyi and Tiab (2007) For the skin factor, the following equation from Omosebi and Igbokoyi (2015) is obtained: Wellbore storage is computed with the usual equation.
Permeability in region II and the skin factor can be computed with the following equations presented by Tiab (1995

3
Skin factor S using Newtonian formula:

Discussion of the results
The new analytical model developed was used to analyze simulated data from the work of Lund and Ikoku. The major differences between Lund and Ikoku, and this study results are the value of flow index n in region I and the polymer front. The analytical method for estimating radius of investigation developed by Ikoku and Ramey could not confirm the polymer front estimated by Lund and Ikoku. While Lund and Ikoku used the plot of ΔP versus t v to estimate the polymer front, we used pressure derivative curve to estimate the polymer front, which should be more accurate. Moreover, we used Odeh and Yang's method to estimate the polymer front, and the results are in good agreement with this study's. Permeability estimated from both type curve and TDS methods for injection and falloff tests is consistent. Therefore, the new analytical method developed is considered valid.

Conclusions
A general solution for a composite reservoir with two regions in a radial system is developed. The solution is applicable to both non-Newtonian and Newtonian composite reservoirs in two-region radial system. The model was applied to solve simulated numerical examples of injection and falloff tests. Type curve and Tiab's direct synthesis methods were used to analyze the numerical data, and consistent results were obtained.
The radius of investigation for non-Newtonian derived by Ikoku and Ramey, and Odeh and Yang is based on the same assumptions. However, Odeh and Yang's method gives more accurate results.

Appendix: Mathematical derivations
The general partial differential equation (Ikoku & Ramey 1979) Lund and Ikoku (1981) to non-Newtonian/Newtonian fluid flow, at interface R D (i.e., rates are equal at the interface): Region 1: Non -Newtonian 1 ≤ r D ≤ R D Laplace solution for the dimensionless pressure drop in infinite (non-Newtonian) flow for region II is: From (60) and (68), At interface,r D = R D , P D1 = P D2 Therefore, Set: where: Combining Eqs. 69 and 70, we have: At wellbore condition, r D = 1 ∶ Combine Eqs. 77 and 79, we have: Laplace solution for the dimensionless pressure drops for constant pressure boundary system for region II is Therefore, Also,
Code availability Available.

Conflict of interest Not applicable.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.