A 2-D analytical solution to infinite conductivity fracture injection pressure transient analysis in non-Newtonian/Newtonian composite infinite reservoirs

Two fluid regions are usually created during either water or polymer injection into an oil reservoir for enhanced oil recovery. The two fluid regions create a composite reservoir system that requires analytical solution that can be used to evaluate the injection performance. In this study, we present a 2-D solution for non-Newtonian fluid in the inner zone coupled with infinite boundary in the outer zone using power law model. The outer zone fluid can either be Newtonian or non-Newtonian. Moreover, the mathematical model is formulated such that either zone can be with non-Newtonian power law fluid or Newtonian fluid. The infinite conductivity fractured solution obtained was compared with the ones obtained using the radial line source solution that had already been developed in the literatures, and a perfect match was obtained. An application was made to three examples that have been presented in a previous work, and satisfactory results were obtained. The dimensionless pressure drops during the early linear flow regime for all values of flow indices are the same, but with deviation during radial flow regime.


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 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 Pi Initial pressure, Pa P wD Dimensionless wellbore flowing pressure P D Dimensionless pressure drop in Laplace space P WD Dimensionless pressure drop in Laplace space q Volumetric flow rate,m 3 /sec g w G Radial distance at wellbore, m R D Dimensionless radial distance at any point, (r/x f ) S Skin factor, dimensionless t Time, s t*ΔP' Pressure derivative, Pa t D Dimensionless time x, y Length in x and y directions, respectively x f Fracture half length l Matrix characteristic length j Number of normal sets of fractures; j = 1, 2, 3. Γ(x) Gamma (factorial) function ϕ Porosity, fraction of bulk volume λ Interporosity flow parameter, dimensionless μ Viscosity, pa-s θ Effective mobility ratio for a two-region composite reservoir ω Storage capacity ratio, dimensionless π 3.142

Introduction
Injecting foreign immiscible fluid into a reservoir or porous media will create a composite feature within the reservoir. In the oil exploitation, polymer, water or waste slurry can be injected into a reservoir for the purpose of enhanced oil recovery or waste disposal. Monitoring the injected fluid's front is therefore necessary. Lund and Ikoku (1981) presented a numerical radial model for a composite reservoir with non-Newtonian fluid in the inner region (Zone 1) and Newtonian fluid in the outer region (Zone 2). Martha and Ertekin (1983) presented a Numerical Simulation of a Power-Law Fluid Flow in a Vertically Fractured Reservoir with finite conductivity. Ertekin et al. (1987) also presented another numerical method with similar composite reservoir as in Lund and Ikoku, but used finite conductivity hydraulic fracture. Prior to Ertekin et al.'s work, Ikoku and Ramey (1979), and Odeh and Yang (1979) have presented radial analytical line source solutions for a non-Newtonian fluid flow in porous media using power law model. The two solutions are similar and essentially the same in analytical approach. Vongvuthipornchai and Raghavan (1987) presented a numerical method for evaluating infinite conductivity hydraulic fracture with non-Newtonian power law fluid flow in porous media. They observed that the dimensionless pressure drops during early linear flow regime are the same for all flow behavior indices. Raghavan and Chen (2017) also presented a new analytical line source solution for non-Newtonian fluid, which was used to compute the dimensionless pressure drop for infinite conductivity hydraulic fracture. The early linear flow is not present in the work of Raghavan and Chen (2017) as observed in the numerical work of Vongvuthipornchai and Raghavan (1987). In this study, we developed a 2-D analytical solution for infinite conductivity hydraulic fracture in a composite reservoir with non-Newtonian fluid in the inner region (Zone 1) and Newtonian fluid in the outer region (Zone 2) as shown in Fig. 1. It is a general solution which can also be used for non-Newtonian fluid with different flow behavior index in Zone 1 and Zone 2. To the best of our knowledge, this analytical solution for a composite reservoir is appearing in the literature for the first time. Previous analytical works (such as Ambastha (1988) and, Abbaszadeh and Kamal (1989)) are on composite reservoir with Newtonian fluids in both Zone1 and Zone 2. van den Hoek et al. (2012) provided a method of analyzing pressure transient behavior in fractured polymer injection well. van den Hoek's method of evaluation is to identify the non-Newtonian radial flow regime in the inner region (Zone 1) and fixed a straight line of the radial flow regime to estimate the flood front. van den Hoek et al. did not provide a matching in their field examples to capture the entire behavior of the pressure transient data.
We compared our 2-D analytical solution with the results obtained using Ikoku and Ramey (1979) and Odeh and Yang line source solutions to obtain infinite conductivity hydraulic fracture solution; we obtained a good match. However, our result is at variant with that of Raghavan and Chen (2017) when the flow behavior index is less than 1.0. We also used the 2-D solution we obtained, applying Cinco-Ley and Meng's (1988) method to compute the dimensionless pressure drop for finite conductivity fracture and compared the results with that of Ertekin et al. (1987); we obtained a fairly perfect match that gave us confidence in our solution. van den Hoek et al. (2012) presented three field examples which practically demonstrated the problem we are addressing in this study. We used our 2-D model to match the pressure transient data and a good match was obtained. Thus, our 2-D model presented a new analytical method for analyzing fractured polymer or waste/drilling cuttings slurry injection well where two fluid zones are created due to the mobility and diffusivity ratios, and different flow behavior index in the two zones. Essentially, the major contribution of this study is the presentation of analytical solution that can be used to evaluate pressure behavior in a composite reservoir with two zones of different flow behavior index. It could be either non-Newtonian/Newtonian or non-Newtonian/non-Newtonian as the case may be. In Well Test Analysis, methods of estimating formation parameters are usually derived from the analytical solution.
The analytical solution in Laplace space is numerically converted to real-time space using Stehfest (1970) Algorithm 368.

Mathematical model
The linearized partial differential equation for naturally fractured reservoir (NFR) is All the dimensionless variables and derivations are in "Appendix A".
The assumptions used in developing this model are: 1. Fluid flow in the inner region (Zone 1) is governed by the infinite conductivity hydraulic fracture model in Eq. 1. 2. The reservoir is infinite. 3. Both regions could exhibit non-Newtonian flow. 4. Laminar and Darcy's flow is valid. 5. Isothermal, single-phase and slightly compressible fluid with constant properties. 6. Homogenous and isotropic system. 7. Pseudoplastic (shear-thinning) fluid obeys Ostwald de Waele power-law relationship. 8. Steady state effective viscosity. 9. The mathematical model also assumed there is no transition region. 10. Warren and Root (1963) matrix pseudo-steady state transfer function is used in the mathematical model.
The solution to Eq. 1, applying to the inner region (Zone 1), after using a g-transformation as shown in "Appendix A" is x Di and y Di is the coordinate at the flood front. When n = 1 in Eq. 6, a Newtonian radial coordinate is obtained. (1) We have assumed homogeneous isotropic formation; hence, the relative permeability in both directions is taken to be unity.
The solution for an infinite system (Carslaw and Jaeger, 1959) is: The governing equation for Zone 2 in an infinite system is Wellbore storage C D and skin S can be incorporated in the wellbore solutions in the Laplace space as shown in Eq. 9.

Type curve and validation
We used the radial line source solutions developed by Ikoku and Ramey (1979), and Odeh and Yang (1979) to develop the solution for an infinite conductivity hydraulic fracture and compared the results with our 2-D analytical solution. We also used Eq. 34 in the work of Rghavan and Chen (2017) which is the analytical solution developed in their work for an infinite conductivity hydraulic fracture to obtain similar curve for flow behavior index n of 0.2, 0.6 and 1.0. The results obtained are shown in Fig. 2. As shown in Fig. 2a, b, Raghavan and Chen (2017) solution did not match other solutions at flow behavior index n of 0.2 and 0.6 except for 1.0. The reason for this discrepancy could be due to the line source solution used by Raghavan and Chen (2017) or in the conversion of their line source solution to that of an infinite conductivity hydraulic fracture.
We also compared our results with that of Ertekin et al (1987) for a finite dimensionless conductivity value of π. The result is shown in Fig. 3. Ertekin et al (1987) used numerical method which probably may have contributed to the observed variation. However, the pattern remained the same.
The dimensionless pressure drop for various values of n computed with the infinite reservoir described in Eq. 7 for both homogeneous and naturally fractured reservoir shown in Fig. 4a, b, respectively.
The plot of the dimensionless pressure drop of non-Newtonian fluid with different flow behavior index (Fig. 4) reveals that the early pressure behavior is practically the same and exhibits linear flow except in the late radial flow regime. Having validated our solution for a non-composite reservoir, we proceed to develop the solution for a composite reservoir.

Composite reservoir
The two equations used in the modeling of the composite system with an infinite conductivity hydraulic fracture in the inner zone are Eqs. 5 and 8. The details of the derivation are in "Appendix A" under composite system. The final wellbore dimensionless pressure drop equation in Laplace space is Eq. 56 in "Appendix A". Figures 6 and 7 are generated with R D = 2.5, alpha = 0.001 and theta = 0.1. Figure 5 demonstrates the effect of mobility and diffusivity ratios as in the work of Ambastha (1988) and Abbaszadeh and Kamal (1989). The value of n and m in Fig. 5 is 1.0 and the value of the pressure derivative must be 0.5 at the longtime radial flow when theta value is unity (1.0). As defined in "Appendix", alpha and theta represent, respectively, the diffusivity and mobility ratios of Zone 1 to Zone 2 as shown in Eqs. 50 and 51. Since Eq. 56 is a general solution, it must replicate the case of Newtonian fluid in both zones and this is shown in Fig. 5.  Figure 6 is for homogeneous reservoir and naturally fractured reservoir. The pressure derivatives at the late radial flow in Zone 2 is less than 0.5 because of the value of theta (mobility ratio) used which is less than 1.0. As shown in Fig. 5, when the mobility ratio is less than 1.0, the dimensionless pressure derivative of the late radial flow regime will be less than 0.5. When the mobility ratio is greater than 1.0, the dimensionless pressure derivative of the late radial flow regime will be greater than 0.5. As shown in the work of Abbaszadeh and Kamal (1989), the mobility ratio has a stronger influence on the pressure derivative behavior. There are many non-Newtonian viscosity models (Savins (1969)). Our mathematical model is only developed for OstWaal-de-Waele power law model. However, it is applicable to Newtonian model by setting the flow behavior index to unity. 2. The model assumed no transition region. 3. Warren and Root (1963) matrix pseudo-steady state flow into the fracture is assumed. This may not be applicable to matrix region of high natural fracture intensity where complex models like trilinear and triple-porosity models exist. Van den Hoek et al. (2012) presented three practical examples that we used to practically validate our solution. Table 1 shows the parameters used in Eqs. 56 and 9 to generate the type curves in obtaining the match for the three examples. The injected polymer is in Zone 1 with flow behavior index n, while Zone 2 contains the in situ fluid with flow behavior index m equal to 1.0. This is a case of non-Newtonian and Newtonian fluids in Zone 1 and Zone 2, respectively. If R D is the radial distance of the flood front, assuming a square model, then,

Application to practical example
Example #1 Fig. 9 shows the type curve match for the falloff data in a polymer injection well after few months of injection. During radial flow, the dimensionless wellbore pressure drop and its pressure derivative for a non-composite reservoir are given by Eqs. 11 and 12, respectively.
On a log-log scale, a straight line can be fixed on the pressure and pressure derivative data. The slope of such straight line on the pressure derivative data is v from which the value of n can be determined Using the dimensionless definition in Eqs. 21 and 22 in "Appendix", Eq. 12 becomes: Equations 11-14 are derived from Eq. 43. The factorization and inversion of Eq. 56 for a longtime case is presented in "Appendix". As shown in Eq. 64 in the "Appendix",  Equations 64-66 show that A is a function of the position of the flood front, mobility and diffusivity ratios and flow behavior index n and m. We fixed Eq. 15 on the dimensionless pressure derivative curve by computing the dimensionless pressure derivative from the values of t D, using Eq. 12. The numerical value of A is determined by trial and error until the log-log straight line is perfectly on the non-Newtonian radial flow region of the dimensionless pressure derivatives. For example, Fig. 8 is the type curve for example #1. The straight line with the red color was drawn with Eq. 12, while the straight line with the green color was drawn with Eq. 15 with the A value of 0.8. Figure 9 shows the complete type curve match, though without the red line. The reason for Fig. 8  Refer to Table 1 for the parameters used to generate the type curve. In this example, both the type curve match and van den Hoek et al. (2012) Fig. 10 shows the type curve match for this example. The pressure data are from different injection well, and the polymer injection has been on for about nine months Example #3 This injector is the same as in example #1. The type curve match is shown in Fig. 11. As indicated by van den Hoek et al. (2012), using the workflow for examples #1 and #2 is not possible because the large wellbore storage has eliminated the non-Newtonian radial flow in Zone 1 (polymer zone). van den Hoek et al. (2012) assumed a flow behavior index of 0.6 to draw the straight line. This might lead to uncertainty in the estimation of the flood front. Using a type curve match method should give a better result with lower degree of uncertainty. The flow behavior index we used is 0.8. The reason for not being able to match the early part of the pressure data could be due to the drop in the liquid level during the falloff test, a phenomenon usually experienced in the phase segregation case.

Example #2
The point of the intersection of the straight line with the pressure derivative curve at best can only be the transition zone. The flood front radial dimensionless distance used in the type curves is provided in Table 1.

Discussion of the results
The examples #1 and #3 are from the same well. Therefore, we do not change the values of alpha and theta. We, however, maintained the same values for example #2.
Having developed the analytical 2-D mathematical model, our goal is to use the model in evaluating polymer flood front. van den Hoek et al.'s (2012) examples provide appropriate case for practical validation. Formation parameters can be evaluated from both zones if injection and necessary formation data are available. Using Tiab Direct Synthesis (TDS) for evaluation has been discussed extensively by Katime-Meindl and Tiab (2001), Igbokoyiand Tiab (2007) and, Omosebi and Igbokoyi(2015). However, as indicated in Eq. 64 in "Appendix", a correction factor is needed for correcting the pressure derivative equation of non-composite system to that of a composite system. A major highlight from this study is the fact that the pressure derivative data in a composite system during non-Newtonian radial flow are strictly numerically not equal to its equivalent in the non-composite system. As earlier indicated by van den Hoek et al. (2012), a correction factor is needed to adjust the pressure derivative log-log straight line equation for non-composite system to that of a composite system before Eq. 12 can be used in the analysis or application of the TDS formulae from non-composite system.
From the three examples, we observed that as n tends to 1.0, the value of A in Eqs. 15 and 16 tends to 1.0. We have used the available field data to prove the validity of our mathematical model for the case of reservoir flooding that result in two flow behavior indices in the reservoir. Even though the flow behavior index we used in example #2 and #3 is different from that of van den Hoek et al.'s (2012), the point of intersection of the two straight lines with the pressure derivative curve on the match is practically the same in each of the examples #2 and #3. As earlier mentioned, both indicated the same transition zone in each case.

Conclusions
A 2-D general solution for evaluating fractured injection well performance with non-Newtonian fluid in a composite reservoir is developed. The computed dimensionless pressure drop is favorably compared with the wellbore infinite conductivity dimensionless fractured pressure developed from the line source solutions of Ikoku and Ramey (1979), and Odeh and Yang (1979). We also used the work of Ertekin et al (1987) to validate our solution for a finite conductivity case. The composite mathematical model in an infinite system was applied to solve the examples presented by van den Hoek et al. (2012). Although the match obtained in examples #2 and #3 is with higher values of flow behavior index n, we have confidence in the overall results because of the consistency in the flood front that agreed with the duration of injection.
The 2-D solution is applicable to naturally fractured reservoirs and also the case of non-Newtonian/non-Newtonian composite system as in polymer flooding of a heavy crude reservoir where the heavy crude exhibits non-Newtonian behavior. As indicated by van den Hoek et al. (2012) and the longtime inversion of the analytical solution for the composite system in this study, a correction factor is needed to apply the TDS formulae developed for non-Newtonian radial flow from the non-composite system.

I: Infinite conductivity hydraulic fracture
The general 2-D non-Newtonian fluid flow equation for the proposed model is The dimensionless form of Eq. 18 is given as: From Warren and Root (1963), Using the same derivative transformation for y D , apply to Eq. 28 with relative permeability of unity gave: Define v as: Solution to Eq. 33 in an infinite system is: For infinite conductivity hydraulic fracture: By making At the early linear flow, sfs is large and is approximately sω. By using the expression in Raghavan and Chen (2017) in Eq. 40 of their paper, we can show that At early time, the expression in Eq. 38 for all values of n 0 ≤ n ≤ 1 gave: For homogeneous formation when ω = 1;

II: Infinite Conductivity Hydraulic Fracture in a Composite Reservoir
The two equations for the composite system in an infinite system are: for the inner region with infinite conductivity hydraulic fracture and, for the outer region. and (39a) where n and m are the flow behavior indices in the inner and outer region, respectively. By considering g spatial variable as an equivalent radial distance (Vongvuthipornchai and Raghavan 1987;Equation B-9).
At the interface P D1 = P D2 .Therefore, Combining Eqs. 48 and 53, at the interface; Combining Eqs. 52 and 54. Note that when m = 1.0, v 2 is zero and A 1 vanishes. A is a factor to be determined by trial and error as described in the three examples.