Boltzmann transformation of radial two-phase black oil model for tight oil reservoirs

Tight oil accumulates in impermeable reservoir rocks, often shale or tight sandstones. The flow behaviour of tight oil in unconventional reservoirs is described by peculiar complexities such as the typical low permeability to viscosity ratio and the dissolution of some gases in the oil phase. Reservoir simulations that consider these complexities negligible stand the potential of poorly characterizing the reservoir flow dynamics. The adoption of similarity transformation effectively reduces the complexities associated with the flow equations through spatial variable (r) and temporal variable (t). The Boltzmann variable ξ=rt\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left(\xi =\dfrac{r}{\sqrt{t}}\right)$$\end{document} is introduced to facilitate the reformulation of transient two-phase flow phenomenon in a radial geometry. The technique converts the two-phase Black oil model (thus highly nonlinear partial differential equations (PDEs)) to ordinary differential equations (ODEs). The resulting ODEs present a reduced form on the flow model which is solved by finite difference approximations (the Implicit-Pressure-Explicit-Saturation (IMPES)) scheme. No loss of vital flow characteristics was observed between the Black oil model and the similarity transform flow model. Furthermore, the similarity approach facilitated the determination of pressure and saturation equations as unique functions of the Boltzmann variable. This derivation is applied to an infinitely acting reservoir where the Boltzmann variable tends to infinity (ξ→∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi \rightarrow \infty$$\end{document}). Finally, this case study’s analytical solution formulated critical relations for fluid flow rate and cumulative production, which are useful for single-phase flow analysis.


Introduction
The decline in conventional hydrocarbon resources coupled with the increase in energy demand has encouraged the development of unconventional resources. The production of oil from conventional resources has peaked and is currently on a terminal, long-run global decline (Gordon 2012). According to U.S. Energy Information Administration (2013); Kuppe et al. (2012), tight oil systems are examples of unconventional resources that provide a significant amount of petroleum for the world's energy needs. Literature highlights several definitions for unconventional resources based on varying factors such as reservoirs with permeability threshold less than 0.1mD (Kazemi 1982;Belhaj et al. 2003;Alkuwaiti et al. 2021); nature of basin and trap mechanisms (Meckel and Thomasson 2008); and fluid types (maturity, density and viscosity) (Thakur and Rajput 2011). In this paper, the graphical definition proposed by Cander (2012) that incorporates both rock property (permeability) 1 3 and fluid property (viscosity) is adopted. Therefore, unconventional reservoirs define reservoirs whose permeabilityviscosity ratios require appropriate techniques to alter the rock permeability or the in situ fluid viscosity to achieve commercially viable production rates.
Production of tight oil comes from very low permeability rocks that must be stimulated using hydraulic fracturing mechanisms to create sufficient permeability for matured oil and natural gases to flow at economical rates (Firincioglu et al. 2012;Orangi et al. 2011). The low permeability of tight oil reservoirs requires production with large pressure drawdowns. The pressure drawdown is often large enough to engender the drop of flowing pressures below the bubble point of in situ liquids hence, causing the evolution of dissolved gases (Aziz and Settari 1979;Chen and Ewing 1997).
Several mathematical models Peaceman (1977); Aziz and Settari (1979); Dake (1978); Carlson (2006); Chen and Ewing (1997); Yu-Shu (2016) have been used to analyze the flow dynamics of conventional reservoirs. However, the development of models for unconventional reservoirs presents peculiar complexities, as discussed by Aldhuhoori et al. (2021aAldhuhoori et al. ( , 2021b. The black oil model, for example, provides a traditional approach for studying multiphase and multicomponent flow systems with a potential mass exchange between phases in a porous medium (Peaceman 1977). The nonlinearities associated with the models describing multiphase flow in tight oil reservoirs present a challenging task to engineers in finding immediate solutions (Tabatabaie and Pooladi-Darvish 2016). Furthermore, the analysis and understanding of the factors that affect the flow performance of these unconventional reservoirs are critical for their efficient exploitation.
These inherent complexities necessitate the adoption of modelling and simulating approaches that appropriately describe such problems without losing vital information about the flow phenomenon (Marshall 2009). The radial diffusivity equation is widely used in the petroleum industry to mimic hydrocarbon flow with several applications in well testing (Ertekin et al. 2001;Ahmed 2018). Unfortunately, the work of Tabatabaie and Pooladi-Darvish (2016) highlights a transient linear flow model applied in the study of boundary dominated flow. The implementation of similarity variable theory in flow analysis of tight oil reservoirs is one of several approaches to finding solutions to reservoir flow problems. The similarity theory is fundamentally based on a combination of variables peculiar to the equations under consideration. One of such variable combinations is defined by the Boltzmann variable ( ). The Fickian diffusion equation is pivotal in studying concentration-dependent mass diffusion in several fields. It has been established that finding closed-form solutions to problems with concentrationdependent diffusion coefficients were complicated (Kass and O'Keeffe 1966). Boltzmann (1894) proposed a novel analytical approach for evaluating concentration-dependent diffusion coefficients. The study transformed the Fickian diffusion equation from a partial differential equation to an ordinary differential equation, commonly referred to as the B-equation. The B-equation was subsequently used to determine the diffusion coefficients (Okino et al. 2012;Ahmed et al. 2015).
Carslaw and Jaeger (1959) adopted the Boltzmann transformation approach to determine solutions of heat conduction problems. Similarly, Bird et al. (1962) investigated several problems associated with heat and mass transfer and fluid flow. Ayala and Kouassi (2007) also applied the Boltzmann transformation in the study of gas condensate reservoirs where an analytical solution for fluid pressure in terms of saturation was developed. Other investigations have also established varying forms of the saturation-pressure relationship, such as those proposed by Raghavan (1976) using the producing gas-oil ratio (GOR). Behmanesh (2016) developed an alternative solution to establish the saturationpressure relationship by adopting a numerical approach.
However, this paper explores the use of the Boltzmann variable to transform the radial two-phase Black oil flow model into reduced forms. The reduced form (ODEs) and the PDEs were numerically solved using an Implicit-Pressure-Explicit-Saturation (IMPES) finite difference technique. This approach also facilitates the derivation of an analytical solution under the assumption that the reservoir is infinitely acting (a limiting case of short producing time at constant pressure production) for → ∞.

Mathematical formulation
In this paper, the radial diffusivity of two-phase (oil and gas) in porous media is described by a similar form of the Black oil model, henceforth named as such (Peaceman 1977;Aziz and Settari 1979;Yu-Shu 2016). The formulation presents a pair of coupled nonlinear PDEs resulting from the combination of the following three governing equations: (a) the law of mass conservation (the continuity equation), (b) Darcy's empirical law, (c) equation of state (Aziz and Settari 1979;Chen and Ewing 1997;Peaceman 1977;Bear and Bachmat 2012). Based on a set of assumptions, the two-phase radial diffusivity model for oil and gas, respectively, is given as: Equation (2) specifically accounts for both gases dissolved in solution and those that evolve out of the solution when in situ pressures fall below bubble point pressure. The evolution of gases out of solution is described using the oil-gas solubility ratio (R s ).

Similarity transformation of flow model
Prior to the similarity transformation, the parameters defined in Table (1) are adopted to reformulate Eqs. (1) and (2) as follows: A step-wise transformation of Eqs. (3) and (4) into ODEs is implemented by introducing the Boltzmann variable ( = r √ t ) as shown in Appendix (8) and (9). The transformation process yields ODEs for oil and gas, respectively, as follows: Equations (5) and (6) present a coupled system of ODEs that is solved by an Implicit-Pressure-Explicit-Saturation (IMPES) scheme. The adoption of the IMPES scheme necessitate the reformulation and decoupling of Eqs. (5) and (6) into pressure and saturation equations. The derivation of a pressure equation is achieved by first expanding the derivatives of all parameters that are functions of pressure, p and oil saturation S o . Considering the ODE for oil flow given by Eq. (5), the expansion of in terms of p and S o yields: Notably, Eq. (7) has two unknowns, pressure and saturation. If saturation were known, Eq. (5) could be solved directly to find pressure as a function of the Boltzmann variable. However, saturation profile is not known a priori. In order to complete the system of equations and unknowns, Eq. (6), representing the flow of gas, is employed to obtain a saturation equation. Again, the derivation of the saturation equation require the expansion of the derivatives of the parameters R, b and in Eq. (6) to obtain: The saturation equation for the flow problem is finally derived through further re-arrangement of Eq. (8) as:

Saturation and pressure profile of flow model
A similar decoupling process is adopted to reformulate the radial diffusivity of two-phase Black oil model defined by Eqs. (1) and (2). The decoupling of the nonlinear PDEs result to pressure and saturation equations that are analogous to the similarity Eqs. (7) and (9). First, Eq. (1) is multiplied by the binomial term (B o − R s B g ) , which represents the critical formation volume factor of oil in the medium. Equation (2) on the other hand is multiplied by the gas formation volume factor, (B g ) . The addition of the resulting equations yield the right hand side (RHS) given by: whereas the resulting LHS expression is given by: The adoption of the chain rule facilitates the expansion of all the time derivative terms in the RHS Expression (10). Nevertheless, a fundamental assumption on the state of the flow problem requires that the total saturation of oil and gas is 1, ( S o + S g = 1 ). This primary assumption, therefore, causes all the time derivatives of saturation to resolve out of the RHS Expression (10) to yield: From a mathematical perspective, the decoupling approach contributes to the derivation of the following compressibility terms that influence the nature of the flow of the in situ fluids.
(a) Gas compressibility: Additionally assuming that the formation compressibility c f is negligible, the right hand side (RHS) Expression (12) simplifies to: The combination of the derived left-hand side (LHS) Expressions (11) and right-hand side (RHS) Expression (13) as well as the introduction of the parameters given in Table  (1) yields the pressure equation as: In conclusion, Eq. (3) describing the radial flow of oil in the derivation process is used to determine the saturation profile once the pressure profile is obtained from Eq. (14). Similar expansions of Eq. (3) by the chain rule yields a saturation equation given by:

Implicit-pressure-explicit-saturation (IMPES) scheme
The equations associated with the two-phase Black oil model are highly non-linear. Therefore, a numerical approximation technique is implemented to offer a comparative basis to the similarity solution. Literature proposes several schemes like the Fully Implicit (FI), Implicit-Pressure-Explicit-Saturation (IMPES), Adaptive Implicit Method (AIM) and other schemes. The FI is unconditionally stable when solving for the unknown pressure and saturation. However, the implicit nature presents a highly expensive computational approach per timestep compared to the IMPES approach (Marcondes et al. 2009;Chen et al. 2006). The IMPES is considered one of several practical approaches in reservoir simulation studies (Coats 1982;Fanchi 2001;Chen et al. 2006). The procedure facilitates the decoupling of the flow model into a pair of pressure Eq. (14) and saturation Eq. (15). The formulation of a unique pressure equation that is independent of the ordering of variables or ordering of equations is one advantage of implementing IMPES. The pressure equation is solved implicitly while lagging all saturation dependencies to the old time level. After obtaining the pressure solution at a new time level, the saturation equation is solved explicitly (Maciasa et al. 2013). As pressure is the only unknown variable determined by a linear system of equations, the approach requires small computational effort per timestep to implement. However, the explicit part must conform to the Courant-Friedrichs-Lewy (CFL) condition (CFL number < 1) to avoid instability in the saturation solution. For this reason, the IMPES approach is limited to small timesteps (Marcondes et al. 2009). Thomas and Thurnau (1983) proposed the Adaptive Implicit Method (AIM), which combines the inherent advantages of the IMPES and FI schemes to achieve better stability, accuracy, and efficiency. The AIM approach implements the FI method in regions where the IMPES method exhibit instabilities while implementing the IMPES method in the rest of the reservoir (Fagin and Stewart 1966;Fung et al. 1989;Watts and Shaw 2005). Unknown pressure and saturation are primarily determined implicitly in complicated subdomains where the CFL number>1, while only pressure is implicitly determined in less heterogeneous subdomains. A switching criterion based on error estimation at the beginning of each time step is used to control the change of the computational schemes. However, this approach introduces a load imbalance problem which is impractical in parallel computing environments (Lu and Wheeler 2009).

Analytical solution for an infinite acting system ( → ∞)
A reservoir that exhibits no apparent outer boundary effects on fluid flow is termed infinitely acting. In this paper, the condition for infinite-acting reservoir systems corresponds to the Boltzmann variable approaching infinity ( → ∞ ). From a practical point of view, the implications include small values of time or large values of distance since is a function of r and t. The derivation of analytical solutions for pressure and saturation under the condition of large values is developed as follows.
At large values of , the term p → 0 . Hence Eq. (7) reduces to: However, since the saturation profile is not known a priori. The flow equation of gas in Eq. (9) is employed to complete the system of equations needed to determine the unknowns (pressure and saturation).
The infinite-acting flow problem ( → ∞ ), defined by the system of equations (Eqs. (16) and (9)), is further constrained by the following initial and boundary conditions; The constrains given by Eqs. (17), (18) and (19) are applicable for the investigation of an infinite acting reservoir, with constant pressure production at the wellbore. Eq. (17) indicates that, the pressure is constant and equal to the pressure at the producing face or wellbore (p w ) . Equation (18) and (19) representing the reservoir pressure (p e ) and oil saturation (S o ) , respectively, are uniform (thus remain unchanged) at initial time (that is, at t = 0 ) to the far boundary of the reservoir ( r → ∞).

Pressure solution
The development of the analytical solution for pressure for the infinite acting flow problem ( → ∞ ) involves the substituting Eq. (9) into Eq. (16) to obtain the following equation expressed in terms of pressure only: In order to reduce the nonlinearity associated with Eq. (20), the two-phase pseudo-pressure (m) is introduced in the reservoir flow behaviour (Fetkovich 1973;Raghavan 1976;Behmanesh 2016). In this paper, a similar approach is used, and the two-phase pseudo-pressure is expressed as: Let, the first derivative of the two phase pseudo-pressure be defined as dm = i dp which resolves Equation (20) to: Under the unique condition of → ∞ Eq. (22) is reformulated as: where: Replacing the full expressions for and into Eq. (24) gives Eq. (25) is similar to the well known hydraulic diffusivity of a single-phase reservoir, adjusted to reflect the effect of two-phase flow and the evolution of gas out of the oil when the pressure drops below the bubble point pressure. Equation (26) and (27), on the other hand, account for the total compressibility of the fluids and fractional flow of oil, respectively. Several physical interpretations are imbibed by Eqs. (26) and (27) despite the rigorous mathematical approach adopted in their derivation. The different terms of Eq. (26), more specifically, are explained as follows: dR s dp represents the amount (positive quantity) of gas released per unit pore volume at reservoir conditions during the pressure drop of dp. (b) The second term S g B g dB g dp represents the effect of gas compressibility on the flow. (c) The third term The solution of Eq. (23) constrained by the two conditions for pressure (Eqs. (17) and (18)  where m i and m w are two phase pseudo-pressures evaluated at the initial and flowing pressure respectively.
In the event where the variation of ∞ with m is negligible, ∞ is assumed to be a weak function of m. Mathematically, the assumption implies that ∞ = ∞ = constant which redefines Eq. (30) as: Additionally, assuming that all terms in ∞ are constant and equal to their initial values, ∞i , Eq. (31) resolves to:

Saturation solution
The analytical solution for saturation of the infinite acting problem ( → ∞ ) is determined by introducing dm = i dp to reformulate Eq. (9) into: The incorporation of the full expressions of the parameters, , , R and b in Table 1 into Eq. (33) yields:: where; Since → ∞ , the evaluation of the coefficients dm d in Eq.
(34) at their initial values yields: where; The introduction of into Eq. (38) and subsequent integration of the resultant gives: Considering the special value Erfc(∞) = 0 and evaluation S o ( ) at = ∞ gives: Hence, the saturation solution is obtained as;

Production parameters
The derived pressure and saturation solutions are useful for deriving oil and cumulative production equations. The cumulative oil production from hydrocarbon-bearing reservoirs is generally expressed as: where the oil flow rate ( q o ), according to Darcy's equation, is given by Equation (45) is reformulated in terms of the two phase pseudo pressure as: Hence, the substitution of the derivative of the pseudo pressure (Eq. (32)) yields the oil production equation as: Finally, the cumulative oil production given by Eq. (44) can also be expressed:

Results and discussion
The behaviour of tight oil reservoirs during transient radial flow under constant pressure production is studied. The reservoir simulation is conducted for the similarity and two-phase Black oil formulations. The analysis of the flow dynamics is based on the data adopted in Tabatabaie and Pooladi-Darvish (2016). The unique descriptions of the various known parameters necessary for the simulation are presented in Figs. (1, 2, 3, 4, 5 and 6).
In the analysis of the formulations, a cylindrical reservoir of radius(r), 800m and pay thickness(H), 50m, is considered. The porosity and the initial permeability of the flow problem are taken to be 0.1 and 0.01md, respectively. The flow model is initially saturated with oil at a saturation pressure of 50000kPa and produced at a constant flowing pressure of 10000kPa. The fluid properties of the reservoir are presented in Fig. (2, 3, 4). In this paper, the Corey-type relative permeability functions defined in Eqs. (49) and (50) are employed to relate the variation of relative permeabilities to the saturation and illustrated by Fig. (1).
In addition, the endpoint relative permeabilities are taken as 1 regardless of the fluid phase. Similarly, the gas and oil relative permeability exponents are each considered to be 2. Figure (2) through to (6) describe the Pressure-Volume-Temperature (PVT) parameters adopted in this work. Figure  (2) illustrates a typical behaviour of the oil formation volume factor. Below the bubble point pressure, the oil formation volume factor increases with pressure. The volumetric factors (B o ) and (B g ) readily relate the volume of fluid that are obtained at the surface (stock tank) to the volume that the fluid occupy when compressed in the reservoir. Fig. (4) defines the relationship between the solution gas-oil ratio to pressure. For this paper, R s increases approximately linearly with pressure, which is a function of the oil and gas composition. Tight oils contain high amounts of dissolved gas; hence, the oil-gas solubility ratio increases with pressure as observed in Fig. (4) until the bubble point pressure is reached. Upon attaining bubble point, the variation becomes constant, and the oil becomes undersaturated.  Under constant pressure production at the wellbore, reservoir fluids propagate towards the wellbore from the undisturbed regions in the reservoir over time. The constant pressure production of 10000kPa causes, instantaneously, a two-phase flow of oil and gas through the porous media into the wellbore. The two-phase phenomenon ensues since the wellbore pressure is below the bubble point pressure of 50000kPa, hence, causes the evolution of gases out of the solution. The nature of the pressure drop over the domain consequently affects the saturation distribution of fluids (oil and gas). (Figs. 7,8 and 9) During the flow period, it is observed that the pressure profile for each time step does not significantly differ. This phenomenon may be attributed to the assumption imposed on the proportion of oil to gas in the system. The assumption of a slightly compressible flow system causes flow that is synonymous with a single-phase flow system. The transformation of the radial two-phase Black oil model (PDEs) to ODEs was dependent on the condition that all independent variables (r and t) can be combined to a single form defined by the Boltzmann variable ( ). Figures  (10), (11) and (12) compare the solutions of the ODEs and the nonlinear PDEs and showed no vital information lose. The pressure and saturation profiles in real-time domain collapse unto a single curve when plotted versus the Boltzmann variable, as observed in Figs. (10), (11) and (12). All the real-time solutions (p(r, t) or S o (r, t) ) can be readily calculated by taking any point on the plot of p( ) or S o ( ) and assigning their values to a corresponding distance found from = r √ t at any particular time.
In developing the analytical solution for the limiting case of → ∞ , the infinite-acting condition subjects the radial flow regime to take linear flow patterns for larger reservoir extents. The subsequent derivations for oil production rate Despite the development of the analytical solution, the practicality remains stalled. From Eqs. (31) through to (48), the initial gas saturation is zero due to the infinite boundary condition imposed on the flow problem. Hence, evaluating ∞ at initial conditions removes the effect of gas mobility from the analytical solution (since, k rg g i = 0 ). However, it accounts for the effect of gas on oil flow, by changing k ro , and the effect of gas evolution on compressibility, given as, dR s dp i . Naturally, reservoirs with higher gas mobility tend to lose energy support faster than reservoirs with lower gas mobility and produce less. It is, therefore, necessary to take into account the effect of gas mobility in Eqs. (47) and (48) by determining a correction factor. It is worth noting that no assumptions regarding the variation of porosity and absolute permeability with pressure were considered in the radial flow model. The saturationpressure relationship used in this paper are independent of absolute permeability. As discussed by Ramey (1964), the total system compressibility is beneficial for multiphase pressure transient analysis rather than a single-phase compressibility.

Conclusion
The study of radial two-phase flow in tight oil reservoirs is appropriately characterized in this paper by adopting a system of nonlinear partial differential equations that account for the occurrence of both dissolved and evolving gases when in situ pressures fall below bubble point pressure. The evolution of gases out of solution is described using the oilgas solubility ratio (Rs). Specific to this paper, the Boltzmann variable ( = r √ t ) is introduced to facilitate the conversion of the nonlinear PDEs into a system of ODEs. The adoption of this technique enabled a straightforward estimation of unknown reservoir parameters (pressure and saturation). The similarity transforms presented in this paper enabled the estimation of critical reservoir profiles and provided a reasonable alternative to solving reservoir flow problems. The similarity approximation delivered a good mapping following comparison to the radial two-phase Black oil model (nonlinear PDEs) based on the data simulation results. This establishes that the similarity approach is satisfactorily sufficient in solving flow problems of this nature.
The subsequent extension of this paper to understand infinite-acting systems ( → ∞ ) yielded analytical solutions that conform with results obtained by Tabatabaie and Pooladi-Darvish (2016) on linear flow in tight oil reservoirs. Under the limiting condition, similar formulations for oil rates and cumulative oil productions were derived. These were determined by adopting a two-phase pseudo-pressure (m) that linearized the flow problem.
Under several simplifying assumptions to the flow, it was necessary to consider the total reservoir diversities to understand the viability of similarity transformation in reservoir engineering techniques. The data and concepts employed in this paper provides an alternative theoretical representation of flow in tight oil reservoirs.

3
The implementation of the chain rule on Eq. (B15) introduces the derivatives of the Boltzmann variable as given by: Further, the oil equation given by Eq. (A12)) can be expanded to: Substituting Eqs. (A2) and (32) into Equation (A10) gives: Then, substituting Eq. (B24) into Eq. (B23) gives: Finally, expanding and simplifying terms in Eq. (B25) yields a diffusivity equation in terms of the Boltzmann variable. These procedures establish the similarity transform of the gas equation as: Since the PDE Eq. (B26) is only dependent on the Boltzmann variable, it is re-written as an ODE in the form: Consent to participate Not applicable.

Consent for publication 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/.