Analytical Modeling of Multistage Hydraulically Fractured Horizontal Wells Producing in Multilayered Reservoirs with Inter-Layer Pure-Planar Crossflow Using Source/Sink Function Method

Abstract

This study presents a comprehensive analytical modeling technology to model transient behaviors of multilayered reservoirs with inter-layer pure-planar crossflow induced by multi-stage hydraulically fractured horizontal well (MHFHW). The objective of this study is to develop an analytical model for multilayered reservoirs in conjunction with complex MHFHW and to achieve not only accurate and efficient computation, but also well-organized solutions expressed in a systematically integrated manner. The consideration of inter-layer crossflow across adjacent layers sets up the foundation for successful modeling of multilayered reservoirs. Source/sink function method (SSFM) is applied to describe fluid flow. Unsteady-state pressure or production rate solutions of MHFHW with the advantages of fast computation, accurate, and stable solutions are achieved. Comparative and consistent outcomes generated by this work and widely applied industry software have largely enhanced our technical confidence. More importantly, innovatively defined modified dimensionless terms that integrate systematic well-reservoir geometry information, as well as rock/fluid properties of each layer, have been newly applied to regulate the new modified dimensionless rate decline curve. This new technique sheds light on the reservoir characterization practice for complicated reservoir systems. Theoretical results in terms of transient pressure and rate were generated by the proposed multilayered model (SSFM-ML) for five scenarios of general concern, under various reservoir and well parameters, which were examined and discussed to demonstrate technical robustness. Not only does this study give solutions to the targeted multiple layered reservoirs, but it also provides insights into modeling three-dimensional fluid flow in heterogeneous reservoir with complex well configurations. It is recommended that future research should be conducted for more complicated two- and three-dimensional reservoirs, using the similar strategy of developing new type curves through adopting other new forms of modified dimensionless rate and time terms.

1. Introduction

Among the proven hydrocarbon resources, unconventional oil/gas stored in tight sandstone or shale comprises the biggest portion. Tight reservoirs with complicated stratified layers with relatively larger vertically spanned thicknesses, but with poor reservoir permeability and low porosity, have become the central focus of unconventional reservoir engineering. Generally, fluid flow in multilayered reservoir formations is extremely complex in nature due to the diversity of reservoir geological and depositional environment variation over geological history. Recently, horizontal wells with multi-stage hydraulically fractured transverse fractures along wellbore trajectory have become a widely accepted and increasingly implemented technology for unconventional reservoir development, which further complicated the wellbore and reservoir system. Generally speaking, the forming of multilayered reservoirs is attributed to the deposition of rock materials over a long geological history. A typical setting of such a well-reservoir system is schematically illustrated by Figure 1, where the artificial fractures, hydraulically initiated and supported by proppant, are emitting from the horizontal wellbore as shown.

Lee et al. [1] pointed out that one of the main challenges reservoir engineers face is producing oil and gas from multiple vertical layers with different reservoir properties. They concluded that a simplistic reservoir description for layered reservoir always results in an overestimated production potential. Frantz et al. [2] also concluded that a proper description for multi-layered reservoir offers a better evaluation for post-fracture performance in comparison with conventional single-layer reservoir description for a complex tight formation. Establishing a physical-based multilayered reservoir model produced by multi-stage hydraulically fractured horizontal well (MHFHW) in order to provide a practical and reliable tool for reservoir engineers to deal with challenging testing or production data is necessary to help efficiently and economically develop tight oil/gas.

Oil/gas flow problems in porous media are mainly addressed by numerical or analytical methods. However, they both have their own limitations. Finite difference is the most commonly used numerical method to solve flow problems in reservoir engineering by discretizing the reservoir into many grids and its capability to handle complex heterogeneous reservoir model has gained popularity. However, because of hydraulic fracturing applied in a tight formation, the existence of artificial fractures along wellbore causes many difficulties for numerical solutions. Due to the large difference between fracture width and length, the gridding system assigned on fractures has to be extremely narrow but long when extended (fracture width no wider than 2 cm normally but at least 10 m long in grid length required). As a result, extensive computation, slow convergence, and solution stability become common problems confronting reservoir simulation engineers. On the other hand, analytical solutions, although accurate, stable, fast computing, and able to deal with fractures comparatively much easier, are mainly derived under the assumption of homogeneous reservoir and regular-shaped reservoir geometry. Modeling complex heterogeneous reservoirs analytically in a straight manner is nearly impossible.

For homogeneous reservoir models under various types of initial and boundary conditions, many researchers have established fundamental analytical solutions to describe fluid flow behaviors. Carslaw and Jaeger [3] applied SSFM to address heat conduction problems. Gringarten and Ramey [4] first provided instantaneous source functions corresponding to various boundary conditions describing transient pressure behaviors in reservoir engineering. Gringarten et al. [5] extended the point source solution to transient pressure behaviors of a producing/injecting fractured well with a single vertical fracture possessing infinite conductivity by integrating point source solution along fracture plane with respect to time. Ozkan and Raghavan [6] derived point source solution in Laplace domain with various well configurations. Rbeawi and Tiab [7] formulated an analytical model describing pressure behaviors of MHFHW with partially-penetrating hydraulic fractures. In addition to pressure behavior, many studies have examined the production rate response of fractured well in recent years. Mederios et al. [8] developed a semi-analytical model incorporating key features of hydraulic fractures, naturally fractures (dual-porosity system) and wellbore flow to explore production-decline characteristics of fractured well in terms of transient-productivity index. Nobakht et al. [9] provided a simple method of forecasting production in tight/shale gas reservoirs, which can be applied to MHFHW based on single fracture linear flow solution. The solution is simple and easy to generate but may not be able to handle the interference among multiple hydraulic fractures reliably. Additionally, for heterogeneous reservoirs, many studies have also tried to apply analytical methods with newly proposed innovative modeling strategies to achieve reliable and accurate analytical modeling outcomes. Zhao and Thompson [10] developed semi-analytical solutions using SSFM based on imaging principle to model pressure responses and flow characteristics in reservoirs with complex geometries, such as T-shaped, splay, and linear composite reservoirs. Medeiros et al. [8] applied Green function methods to develop a linear heterogeneous model with a horizontal well, intercepting a number of blocks, and examined transient pressure behaviors of the horizontal well. Although the well-reservoir system is simplified without the consideration of hydraulic fractures emitting from wellbore, their pioneering exploration of modeling heterogeneous reservoir is encouraging due to Green function’s analytical nature in handling the source and boundary. Zhao [11] applied the SSFM in developing the technology of modeling reservoir heterogeneity with complicated wellbore and fracture system by enhancing the application of the concept of source/sink and imaging source/sink extensively. The patented technology greatly improved the technical capacity of reservoir engineers in dealing with reservoirs with complicated heterogeneity and geometry. Zhao [12] developed a stimulated reservoir volume (SRV) model for MHFHW, based on the patented modeling methodology [11], which confines a horizontal well with multiple stages of hydraulic fractures within SRV region, using SSFM with solutions implemented in the Laplace domain, and established type curves of pressure and newly modified dimensionless type curves of production rate to help diagnose fracture length in a systematic manner.

It is believed that accurate evaluation of crossflow flux in multilayered reservoir systems plays a central role in modeling such kind of complicated reservoir systems. In literature, a great deal of research on multilayered reservoir modeling with the consideration of crossflow had been conducted with various methods. A brief summary of the most related studies is summarized in

Table 1. Literature regarding multilayered reservoir modeling with crossflow consideration.

Author(s)	Year	Well Type	Reservoir Type	Documented Solution(s)
Russell and Prats [14]	1962(a)	Vertical	Two layers	P~t, q~t
Russell and Prats [15]	1962(b)	Vertical	Two layers	P~t, q~t
Neuman and Witherspoon [16]	1969	Vertical	Two layers	PD~tD
Bourdet [17]	1985	Vertical	Two layers	PD~tD
Prijambodo et al. [18]	1985	Vertical	Two layers	PD~tD
Gao [19]	1987	Vertical	Multiple layers	PD~tD, qD~tD
Wijesinghe and Kececioglu [20]	1988	Vertical	Multiple layers	PD~tD, qD~tD
Gao and Dean [21]	1988	Vertical	Multiple layers	PD~tD, qD~tD
Onur and Reynolds [22]	1989	Vertical	Two layers	PD~tD
Ehlig-Economides [23]	1993	Vertical	Multiple layers	ΔP~t
Kuchuk and Habashy [24]	1996	Horizontal	Multiple layers	PD~tD
Sun et al. [25]	2003	Vertical	Two layers	qjD~tD, j = 1, 2
Al-Ajmi et al. [26]	2003	Vertical	Multiple layers	PD~tD
Wang et al. [27]	2005	Fracture, infinite conductivity	Three layers	PD~tD
Villanueva-Triana and Civan [28]	2013	Vertical	Multiple layers	P~t. P~r
Sun et al. [29]	2017	Fracture, finite conductivity	Three layers	PD~tDA under PSS
Lu et al. [13]	2019	Vertical	Two layers	PD~tD
This work	2022	MHFHW, infinite conductivity	Multiple layers	PD~tD, qD~tD, qDM~tDM

. Due to the comprehensive nature of the studies in this technical domain, readers are suggested to refer to Table 1 for further details. Most recently, Lu et al. [13] studied vertical well producing from an infinite two-layer reservoir analytically through using Laplace transformation, double Fourier transformation, and Green’s function method, which is a promising integrated solution system for complicated well-reservoir systems. The above-mentioned analytical models for heterogeneous reservoirs have advanced a significant step in applying SSFM or Green function to describe fluid flow in oil/gas reservoir, however, clear and organized solutions for MHFHW within multilayered reservoir with inter-layer crossflow under three-dimensional flow have not yet been proposed.

The objective of this study is to develop an analytical model for multilayered reservoirs in conjunction with complex MHFHW and to achieve not only accurate and efficient computation, but also well-organized solutions expressed in a systematically integrated manner. The initiation of this study is motivated by the powerful capacity of SSFM to describe fluid flow induced by fractures in reservoir and by the innovative strategy of modeling heterogeneous reservoirs analytically [10], and the creative strategy that has successfully re-organized the general dimensionless rate and dimensionless time type curve using the uniquely new modified dimensionless rate and modified dimensionless time [30]. This new technique sheds light on the reservoir characterization practice for these kinds of complicated reservoir systems.

2. Methodology

Illustrated in Figure 2, to facilitate demonstrating the modeling process, a simplified layered reservoir consisting of only two vertically stacked layers is exemplified and produced by a MHFHW with different vertical and horizontal fracture penetration in each layer with the lower layer wider and extended further. A layered reservoir with more layers can also be modeled accordingly. The two layers are in hydraulic communication with each other along their interface. The lower and upper layers are individually homogeneous porous mediums with constant vertical thickness $h_{e1}$ and $h_{e2}$ and have rock/fluid property ${\left(k,\varphi ,\mu ,c_t\right)}_1$ and ${\left(k,\varphi ,\mu ,c_t\right)}_2$, respectively. A horizontal well was placed within one of the layers, with N_stage number of hydraulic fractures transversely distributed along the wellbore. The fractures with infinite conductivity are assumed to penetrate into both layers. Each fracture is assumed to have the same fracture length: $2L_{f1}$ and $2L_{f2}$ in the lower and upper layer, respectively. Note that non-uniform fracture length can be assumed if needed. The fracture heights penetrating into each layer are denoted as $h_{f1}$ in Layer 1 and $h_{f2}$ in Layer 2. In this study, the production contribution of horizontal wellbore trajectory is neglected due to the fact that fluid production from fractures dominates.

A general multilayered reservoir assumes that multiple rectangular-shaped blocks with different geometries, stacked on top of one another within each being internally homogenous, and a no-flow outer boundary constitute the physical reservoir model. The modeling method consists of decomposing the original layered reservoir into a set of single-layer homogeneous sub-reservoirs that interact with each other over their hydraulically contacted interface(s) by fluid transfer (inter-layer crossflow) and pressure continuity. Hydraulic fractures are assumed to have infinite conductivity (without further consideration of the effect of fracture width) due to the significant contrast of systematical properties between hydraulic fracture and formation. SSFM was applied to describe the pressure drop at any point within each layer by integrating appropriate point source function spatiotemporally.

In the solution process, pressure drop solution written in real-time domain was converted to Laplace domain using numerical Laplace transformation to help establish the solution algorithm effectively. The methodology of performing numerical Laplace transformation for source/sink functions [10] plays a critical role in this work for the reason that pressure drop caused by plane source in real-time domain involves products of source functions in three dimensions (x, y, z) based on Newman product method [31], and the analytical Laplace transform for such a complicated procedure is not readily available. Each reservoir layer is coupled with its adjacent one(s) over the contact interface(s) based on pressure and flux continuity. The desired pressure response at any point within reservoir, flux along hydraulic fractures, and inter-layer crossflow are then solved systematically in the Laplace domain. Finally, pressure and rate solutions are converted to real time domain numerically.

2.1. Decomposing of the Layered Reservoir

It is necessary to identify simple homogeneous reservoir components that constitute this layered heterogeneous reservoir. Firstly, considering the N_stage transverse fractures stimulated inside the lower Layer 1 only, it behaves as if one well and an inter-layer plane source exist: the original fractured well with hydraulic fractures that reach the top of this layer composes only a portion of the original fractures’ vertical penetration, that is, $h_{f1}$, and an injecting planar interface at the conjunction between the two layers. The planar interface is defined as a pure-plane source because the plane is implemented to account for hydraulic fluid transfer only and is not associated with any fluid storage-related issue. Therefore, the inter-layer fluid crossflow from the upper layer to the lower layer, as illustrated by Figure 3a, is defined according to inter-layer pure-planar crossflow. Similarly, considering the upper Layer 2 only, it also behaves as if there were one well and an inter-layer pure-plane source: the original fractured well with hydraulic fractures with vertical penetration $h_{f2}$ starting from the bottom of layer 2, and a producing planar interface at the conjunction of the two layers, also treated as a pure-plane source, illustrated by Figure 3b.

Accounting the inter-layer pure-planar crossflow fluid transfer rate between the two layers as $q^{cf}$, the original reservoir system can then be decomposed into two homogeneous sub-reservoirs. Decomposing the reservoir components from one another at their junction is achieved by applying “no-flow boundary” concept. For simple reservoir systems that contain wells, the no-flow boundary is achieved by applying well-principled images. For the situation being considered, Layer 1 and Layer 2 are each individually modeled using this strategy, and the processes are graphically presented in Figure 4a,b. The decomposed Layer 1 and Layer 2 systems are subsequently coupled using pressure and rate continuity at their contact interface through the application of the pure-plane concept proposed. This original methodology modeling heterogeneous reservoir was devised by Zhao and Thompson [10].

2.2. Solutions to Layer 1 and Layer 2

For Layer 1 as illustrated by Figure 4a, where $\left(x,y,z\right)$ is defined as the sub-coordinate system, pressure drop at any location, $\Delta p_1$, can be formulated by applying superposition principle. Appendix A gives the mathematical derivation of a plane source solution in detail, and the plane source solution is expressed as Equation (A7). With the same discretizing process of plane source as shown in Figure A2, each original hydraulic fracture plane source is discretized into $m\times n$ sub-plane sources, and the injecting inter-layer pure-plane source is divided into $M\times N$ sub-plane sources (Figure 5), respectively. Therefore, the pressure drops in Layer 1 can be written as

$$\begin{array}{c}\Delta p_1\left(x,x_{e1},y,y_{e1},z,2h_{e1},t\right)=\\ \frac{1}{{\left(\varphi c_t\right)}_1}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}\frac{q_{k,i,j}^{f1}\left(\tau \right)d\tau }{2L_{f1}{h}_{f1}/mn}\left[\begin{array}{c}I{p}_{sx}\left(x_{f,i-1},x_{f,i},x,x_{e1},\tau ,t\right)\\ ·I{p}_{sz}\left(z_{f,j-1},z_{f,j},z,2h_{e1},\tau ,t\right)\\ ·p_{sy}\left(y_k,y,y_{e1},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_1}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}\frac{q_{k,i,j}^{f1}\left(\tau \right)d\tau }{2L_{f1}{h}_{f1}/mn}\left[\begin{array}{c}I{p}_{sx}\left(x_{f,i-1},x_{f,i},x,x_{e1},\tau ,t\right) \\ ·I{p}_{sz}\left(-z_{f,j},-z_{f,j-1},z,2h_{e1},\tau ,t\right)\\ ·p_{sy}\left(y_k,y,y_{e1},\tau ,t\right) \end{array}\right]-\\ \frac{1}{{\left(\varphi c_t\right)}_1}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^N}{\displaystyle {\displaystyle \sum }_{j=1}^M}\frac{2q_{i,j}^{cf}\left(\tau \right)d\tau }{\frac{x_{e1}{y}_{e1}}{MN}}\left[\begin{array}{c}I{p}_{sx}\left(x_{l,i-1},x_{l,i},x,x_{e2},\tau ,t\right)\\ ·I{p}_{sy}\left(y_{l,j-1},y_{l,j},y,y_{e2},\tau ,t\right)\\ ·p_{sz}\left(0,z,2h_{e1},\tau ,t\right)\end{array}\right].\end{array}$$

(1)

The first term on the right-hand side of Equation (1) accounts for the original fractured well in the lower layer, the second term accounts for its image well, and the third term accounts for the injecting plane source. Note that in this study, producing is recognized as “+” while injecting is recognized as “-”. In Equation (1), $q_{k,i,j}^{f1}$ represents production rate of a sub-plane from the k-th hydraulic transverse fracture along horizontal wellbore in Layer 1, and $q_{i,j}^{cf}$ denotes the production rate of a sub-plane from injecting pure-plane source.

Similarly for Layer 2, as illustrated in Figure 4b, where $\left(\xi ,\nu ,\zeta \right)$ is defined as the sub-coordinate system, the pressure drop in Layer 2 can also be generated using the superposition principle. Each original hydraulic fracture plane source is also divided into $m\times n$ sub-plane sources for the purpose of demonstrating the modeling strategy (the sub-plane source number can be different) and the inter-layer pure-plane source is the same as in Layer 1. It is important to note that the inter-layer pure-plane in Layer 1 acts as an injecting plane source, whereas the inter-layer pure-plane in Layer 2 acts as a producing plane source. The pressure drop in Layer 2 can be written as

$$\begin{array}{c}\Delta p_2\left(\xi ,x_{e2},\nu ,y_{e2},\zeta ,2h_{e2},t\right)=\\ \frac{1}{{\left(\varphi c_t\right)}_2}{\displaystyle {\displaystyle \sum }_{k=1}^{Nstage}}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}\frac{q_{k,i,j}^{f2}\left(\tau \right)d\tau }{2L_{f2}{h}_{f2}/mn}\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{f,i-1},{\xi }_{f,i},\xi ,x_{e2},\tau ,t\right)\\ ·I{p}_{sz}\left({\zeta }_{f,j-1},{\zeta }_{f,j},\zeta ,2h_{e2},\tau ,t\right)\\ ·p_{sy}\left({\nu }_k,\nu ,y_{e2},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_2}{\displaystyle {\displaystyle \sum }_{k=1}^{Nstage}}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}\frac{q_{k,i,j}^{f2}\left(\tau \right)d\tau }{2L_{f2}{h}_{f2}/mn}\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{f,i-1},{\xi }_{f,i},\xi ,x_{e2},\tau ,t\right) \\ ·I{p}_{sz}\left(-{\zeta }_{f,j},-{\zeta }_{f,j-1},\zeta ,2h_{e2},\tau ,t\right)\\ ·p_{sy}\left({\nu }_k,\nu ,y_{e2},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_2}{{\displaystyle \int }}_0^t{\displaystyle {\displaystyle \sum }_{i=1}^N}{\displaystyle {\displaystyle \sum }_{j=1}^M}\frac{2q_{i,j}^{cf}\left(\tau \right)d\tau }{x_{e2}{y}_{e2}/MN}\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{l,i-1},{\xi }_{l,i},\xi ,x_{e2},\tau ,t\right)\\ ·I{p}_{sy}\left({\nu }_{l,j-1},{\nu }_{l,j},\nu ,y_{e2},\tau ,t\right)\\ ·p_{sz}\left(0,\zeta ,2h_{e2},\tau ,t\right)\end{array}\right],\end{array}$$

(2)

where $q_{k,i,j}^{f2}$ represents the production rate of a sub-plane from the k-th hydraulic fracture.

Laplace transforms, with respect to time, were taken for Equations (1) and (2). Let L[] denote Laplace operator, according to convolution theorem, and the pressure responses in Layer 1 and Layer 2 in Laplace domain can be written as

$$\begin{array}{c}L\left[\Delta p_1\left(x,x_{e1},y,y_{e1},z,2h_{e1},t\right)\right]=\\ \frac{1}{{\left(\varphi c_t\right)}_1}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}\frac{mn}{2L_{f1}{h}_{f1}}L\left[q_{k,i,j}^{f1}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left(x_{f,i-1},x_{f,i},x,x_{e1},\tau ,t\right)\\ ·I{p}_{sz}\left(z_{f,j-1},z_{f,j},z,2h_{e1},\tau ,t\right)\\ ·p_{sy}\left(y_k,y,y_{e1},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_1}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}\frac{mn}{2L_{f1}{h}_{f1}}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}L\left[q_{k,i,j}^{f1}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left(x_{f,i-1},x_{f,i},x,x_{e1},\tau ,t\right) \\ ·I{p}_{sz}\left(-z_{f,j},-z_{f,j-1},z,2h_{e1},\tau ,t\right)\\ ·p_{sy}\left(y_k,y,y_{e1},\tau ,t\right) \end{array}\right]-\\ \frac{1}{{\left(\varphi c_t\right)}_1}\frac{MN}{x_{e1}{y}_{e1}}{\displaystyle {\displaystyle \sum }_{i=1}^N}{\displaystyle {\displaystyle \sum }_{j=1}^M}L\left[2q_{i,j}^{cf}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left(x_{l,i-1},x_{l,i},x,x_{e2},\tau ,t\right)\\ ·I{p}_{sy}\left(y_{l,j-1},y_{l,j},y,y_{e2},\tau ,t\right)\\ ·p_{sz}\left(0,z,2h_{e1},\tau ,t\right) \end{array}\right] ,\end{array}$$

(3)

and

$$\begin{array}{c}L\left[\Delta p_2\left(\xi ,x_{e2},\nu ,y_{e2},\zeta ,2h_{e2},t\right)\right]=\\ \frac{1}{{\left(\varphi c_t\right)}_2}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}\frac{mn}{2L_f{h}_{f2}}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}L\left[q_{k,i,j}^{f2}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{f,i-1},{\xi }_{f,i},\xi ,x_{e2},\tau ,t\right)\\ ·I{p}_{sz}\left({\zeta }_{f,j-1},{\zeta }_{f,j},\zeta ,2h_{e2},\tau ,t\right)\\ ·p_{sy}\left({\nu }_k,\nu ,y_{e2},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_2}{\displaystyle {\displaystyle \sum }_{k=1}^{N_{stage}}}\frac{mn}{2L_f{h}_{f2}}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^m}L\left[q_{k,i,j}^{f2}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{f,i-1},{\xi }_{f,i},\xi ,x_{e2},\tau ,t\right) \\ ·I{p}_{sz}\left(-{\zeta }_{f,j},-{\zeta }_{f,j-1},\zeta ,2h_{e2},\tau ,t\right)\\ ·p_{sy}\left({\nu }_k,\nu ,y_{e2},\tau ,t\right) \end{array}\right]+\\ \frac{1}{{\left(\varphi c_t\right)}_2}\frac{MN}{x_{e2}{y}_{e2}}{\displaystyle {\displaystyle \sum }_{i=1}^N}{\displaystyle {\displaystyle \sum }_{j=1}^M}L\left[2q_{i,j}^{cf}\left(t\right)\right]L\left[\begin{array}{c}I{p}_{sx}\left({\xi }_{l,i-1},{\xi }_{l,i},\xi ,x_{e2},\tau ,t\right)\\ ·I{p}_{sy}\left({\nu }_{l,j-1},{\nu }_{l,j},\nu ,y_{e2},\tau ,t\right)\\ ·p_{sz}\left(0,\zeta ,2h_{e2},\tau ,t\right)\end{array}\right].\end{array}$$

(4)

Zhao and Thompson [10] devised a nontrivial fast and highly accurate method of generating numerical Laplace transformation, and it is implemented in this work. The complicated Laplace transforms of Equation (3) or Equation (4), thereby, can be obtained.

2.3. Coupling of the Multilayered Reservoir System

Equations (3) and (4) were evaluated at the midpoint of each sub-plane of hydraulic fractures, which gives $\left(2·m·n·N_{stage}\right)$ pressure expressions in Laplace domain. Since production from horizontal wellbore is ignored in this study, the pressure along wellbore was assumed to be uniform, and thereby pressures at interceptions of wellbore and fractures are equal to one another. By equating the pressure at the midpoint of each sub-plane source of hydraulic fractures and equating pressures at interceptions of fractures and wellbore, it yields $\left(2·m·n·Nstage-1\right)$ linear equations. Layer 1 and Layer 2 interact with each other by fluid transfer through pure-plane (inter-layer pure-planar crossflow) and pressure continuity. Pressure continuity is achieved by equating pressure at midpoint of each sub-plane of the planar injection/production source from Equations (3) and (4), which yields $\left(M\times N\right)$ linear equations in Laplace domain. In addition to pressure continuity, material balance must be satisfied as well. The sum of rates from every sub-plane of hydraulic fractures amounts to the well production rate, which yields the last linear equation to be solved.

The above coupling conditions of the 2-layer reservoir-well system yield $\left(2·m·n·Nstage+M·N\right)$ linear equations in the unknown sub-plane rates of hydraulic fractures and of inter-layer pure-planar injection/production source. Each sub-plane rate can be obtained in the Laplace domain and substituted back into Equation (3) or Equation (4) so that the pressure at any desired location can be calculated, including pressure at any horizontal wellbore position. The Duhamel principle can then be readily applied to generate a rate solution in the Laplace domain [12]. Rates or pressures in the Laplace domain can then be inverted to the real-time domain by applying the Stehfest inversion algorithm [32]. For commingling well (without inter-layer pure-planar crossflow), it can be easily modeled by setting the inter-layer pure-plane source rate equaling zero.

3. Model Validity and Result Discussion

To validate the multilayered heterogeneous reservoir model by the proposed SSFM-ML model and examine transient behaviors of such a complex well-reservoir system, the solutions generated by the proposed model in the form of type curves in terms of pressure under constant well rate production or well rate under constant pressure production were checked against various scenarios. Dimensionless groups are defined in Appendix B. Five scenarios are examined in this study. The first four scenarios address 2-layer reservoir under various well-reservoir systems and the fifth scenario addresses a 3-layer reservoir. Solutions generated from commercial well-testing software Kappa (V5.20) and numerical simulation software CMG under some cases are also presented for comparison.

The purposes of the five scenarios are as follows. The first scenario focuses on testing solution stability and serving as a limiting case by setting the fracture length to be extremely short, which makes the fractured well ultimately approach a vertical well, so that the generated solution should also approach that of a corresponding vertical well. The second scenario investigates the effect of mobility ratio of the two layers on pressure and rate behaviors of fractured well and inter-layer crossflow response. Outcomes from comparable two-layer reservoir models of CMG, Kappa’s analytical, Kappa’s numerical, and our work using SSFM-ML were compared to enhance our technical confidence. The crossflow flux profile across the inter-layer pure-plane, generated by SSFM-ML model, was fully investigated. The third scenario conducts sensitivity analysis of fracture length for fractured well response under a 2-layer reservoir model. In light of fostering the field application of the technique, this was performed to demonstrate how the re-organized modified dimensionless rate decline curves appear graphically. The fourth scenario is to examine the behaviors of fractured well under two-layer reservoir model without the effect of inter-layer crossflow by setting inter-layer crossflow rate equaling zero in the model. This simulates the general commingled flow of two-layer reservoirs as a special case. The fifth scenario addresses a three-layer reservoir with a sensitivity analysis of storativity among the three layers to showcase the applicability of the generalized modified dimensionless rate decline type curve of q_DM vs. t_DM.

Scenario 1: Limiting case of approaching vertical well solutions for 2-layer model with inter-layer pure-planar crossflow. Assuming that there is only one transverse fracture along the horizontal wellbore with a very short fracture length, the fracture height is comparatively much greater. The result generated by the model must then approach the behavior of a partially penetrating vertical well or a horizontal well, because the fracture geometry is closely approaching a line source due to the extremely short fracture length.

Figure 6 shows the transient pressure derivative of a one-stage fractured well with various dimensionless fracture half-lengths of $L_{f1D}=L_{f2D}=1, 2.5, 5$ ($\mathrm{with} r_{wD}=1$), and both layers are identical in terms of fluid/rock properties. Additionally presented is a standard solution of a vertical well with limited entry in a homogeneous reservoir from Kappa. When the dimensionless fracture half-length equals 1, the one-stage fractured well behaves closely similar to a vertical well with the same level of vertical penetration as its pressure derivative curve almost overlaps that of the vertical well and the flow regimes are clearly exhibited: early in time, the pseudo-radial flow is presented with a flat pressure derivative; after pressure-transient reaches horizontal boundaries (x-axis), the linear flow (channel flow) with pressure derivative displaying half a slope is evolved and finally the pseudo-steady state (PSS) flow is shown after all boundaries are felt. When dimensionless fracture half-length reaches 5, an early formation linear flow with pressure derivative displaying half a slope in reflecting fracture production behavior is successfully exhibited but lasts only a short period of time. Fracture length is still considered very short, and the pressure derivative starts to show a pseudo-radial flow behavior.

A vertical well-producing reservoir with multiple layers or a horizontal well producing reservoir with multiple linear blocks with distinct properties characterized from well logs is commonly seen in field practice. Figure 7 provides typical pressure derivative responses accounting for a vertical well with limited entry in a reservoir with 2 layers, through using a fracture half-length of $L_{f1D}=L_{f2D}=r_{wD}=1,$ produced from a 2-layer reservoir under 3 cases: (1) wellbore fully penetrates only the upper layer/block; (2) wellbore goes into the lower layer/block with 50% partial penetration; (3) wellbore fully penetrate both layers/blocks. Solutions generated from the 2-layer model of a vertical well by Kappa are also presented in Figure 7. The two sets of solutions show very good agreement, which, to a large degree, demonstrates solution accuracy. One may realize that this scenario also approaches the work of Medeiros et al. [8] for horizontal well intercepting a few blocks. Testing data of vertical or horizontal wells in field practice often show these types of pressure derivative responses. With better understanding of the pressure and production performance features of the existing numerous multilayered wells, more accurate reservoir characterization with less uncertainty could be achieved by applying the analytical strategy of this study.

Scenario 2: Effect of heterogeneity under various${\mathit{M}}_{\mathit{2}}/{\mathit{M}}_{\mathit{1}}$ratios for 2-layer reservoir with inter-layer pure-planar crossflow by MHFHW. The influence of mobility ratio between the two layers is examined in this scenario, with 10-stage fractures penetrating both layers. Figure 8 shows pressure derivative behaviors of the well under various levels of heterogeneities as $M_2/M_1=100, 10, 1, 0.1,$ and the solutions generated by CMG are presented as well. As illustrated in Figure 8, the early formation linear flow, with half a slope pressure derivative, marks the first flow regime, implying that each fracture is produced within its own individual linear-flow-expanding drainage area. Pressure derivative is followed by a sharp rise, with the slope almost reaching units, which is caused by the intense interference among fractures, and the unit slope accounts for the SRV boundary influence. The work of Zhao [12] has indicated the strong multiple fracture interference right after the early time period, and the phenomena of the sharp rise of derivative has also been noticed by Chen and Raghavan [33].

In this scenario, where the fracture vertical penetration is comparatively much smaller than the whole reservoir vertical thickness, the pressure-transient before hitting the nearest boundaries is likely to travel far into the reservoir based on the effective system flowing diffusivity. Therefore, the system behaves as if the well was producing from an overall single-layer homogeneous reservoir. Figure 8 exhibits the flat pressure derivative values during the pseudo-radial flow regime. Zhao and Thompson [10] and Spivey and Lee [34] have pointed out that the stabilizing pressure derivative during pseudo-radial flow in such two-region heterogeneous reservoir should approach the value given by the following formula,

$$\frac{d{P}_D}{dln{t}_D}\to \frac{1}{1+\frac{M_2}{M_1}}$$

(5)

Additionally, to compare solution accuracy from the SSFM-ML model against Kappa’s numerical solutions of the pressure and its derivative for one case of $M_2/M_1=1$, i.e., a homogeneous reservoir, the comparison is provided in Figure 9. However, it turns out that Kappa’s numerical provides a solution, with reduced accuracy for cases of MHFHW produced in a 2-layer reservoir with inter-layer pure-planar crossflow, as indicated by the red lines. Note that the solution from Kappa’s numerical solution requires more vertically stacking layers to simulate the designed reservoir case properly, with limited vertical fracture penetration. For the solutions presented as the red lines, 20 sub-layers were integrated in Kappa’ numerical model to generate the reliable solution, which costs much longer computing time than SSFM-ML model does. The solutions from Kappa’s analytical model of homogeneous reservoir are also presented for reference. One can conclude that the solutions generated by these three models are consistent, and thus enhance our technical confidence in applying them. However, generally speaking, the advantages of SSFM-ML modeling are that the MHFHW can be defined and, most importantly, computed accurately in a three-dimensional multilayered reservoir. Because the crossflow flux can be accounted for more accurately on the inter-layer crossflow pure-plane in the SSFM-ML model, it offers solutions with higher accuracy.

Application of the modified dimensionless terms. In addition to studying the pressure responses, the rate responses of the fractured well producing at constant bottom pressure are generated and provided in Figure 10a. It can be seen that, in general practice, a full comparison and understanding of the rate outcomes is hard to reach, because of the fact that MHFHW configuration can be complex and reservoir heterogeneity is heavily involved. Zhao et al. [30] devised a technique to systematically analyze the sequence of flow regimes of 2-dimensional homogeneous reservoirs in MHFHW rate response by defining new terms of the modified dimensionless production rate $q_{DM}$ and the modified dimensionless time $t_{DM}$, which, in a general physical sense, integrate the geometry of fracture and reservoir system for data analysis in a systematic manner. In this study, further improvement has been achieved by extending this novel concept to 3-dimensional multilayered reservoir and integrating heterogeneous reservoir characteristics, including mobility ($k/\mu$) and storativity ($\varphi c_t$) of each layer, into the definition of $q_{DM}$ defined in Equation (A20) and $t_{DM}$ defined in Equation (A21). Therefore, a multilayered system can be evaluated systematically with fracture geometry, reservoir geometry, and rock/fluid properties information integrated.

The production rate responses presented in Figure 10a are plotted again in Figure 10b by applying the newly defined modified dimensionless terms of $q_{DM}$ and $t_{DM}$. It can be seen clearly in Figure 10b that all curves merge together in early linear and late boundary-dominated flow (BDF) periods. This outcome greatly facilitates the comparison, understanding, and application of the MHFHW rate solutions and offers an obvious advantage in rate decline analysis. This technique can be applied in oil and gas reservoir engineering worldwide to help improve the understanding of MHFHW practice effectively and globally, because the solutions under various scenarios become systematically comparable. The mathematical descriptions of $q_{DM} vs. t_{DM}$ under early linear and late time BDF conditions are given by Equations (6) and (7) below, which are also plotted in Figure 10b as Line 1 and Line 2 accordingly. One needs to realize the fact that the amalgamation of these two lines forms a consolidated reference base for systematical rate decline analysis and universally applicable reservoir characterization strategy.

$$q_{DM, Linear}=\frac{1}{\sqrt{0.5{\pi }^3{t}_{DM}}},$$

(6)

$$q_{DM,BDF}=\frac{1}{2\pi t_{DM}}.$$

(7)

Figure 11 presents the production rate from each layer (layered rates) and an inter-layer pure-planar crossflow rate (from the less permeable layer underneath to the more permeable upper layer) when the well is producing under constant rate conditions. General observations from this figure indicate that the layered rates that display flow characteristics can be categorized into four phases that are strongly influenced by dynamic flow regime evolvement, as shown in Figure 10, and that the overall production rate from the more permeable upper layer is always higher than that of the less permeable layer underneath, i.e., $q_D^{f2}>q_D^{f1}$. It is also observable that the evolving flow regimes around the fractures affect the inter-layer pure-planar crossflow, $q_D^{cf},$ along with time, causing it to steadily increase until a full PSS flow regime is reached. Further scrutinization of Figure 11a reveals that the mathematical expressions of $q_D^{f1}$ and $q_D^{f2}$ during early time Phase 1 can be formulated as Equations (8) and (9). The red dash lines of 1 and 2 are plotted according to Equations (8) and (9) with M₂/M₁ = 10, respectively, and similarly, the red dash lines of 3 and 4 are plotted with M₂/M₁ = 100.

$$q_D^{f1}\left(t_D\right)=\frac{1}{1+\frac{L_{f2D}{h}_{f2D}}{L_{f1D}{h}_{f1D}}\sqrt{\frac{C_{s2}{M}_2}{C_{s1}{M}_1}}}$$

(8)

$$q_D^{f2}\left(t_D\right)=\frac{1}{1+\frac{L_{f1D}{h}_{f1D}}{L_{f2D}{h}_{f2D}}\sqrt{\frac{C_{s1}{M}_1}{C_{s2}{M}_2}}}$$

(9)

Scenario 3: Effect of various levels of horizontal fracture extension to the width of 2-layer reservoir with inter-layer pure-planar crossflow for MHFHW. This scenario explores wellbore pressure and rate responses that correspond to various ratios of $\frac{2L_{f1D}}{x_{e1D}}=\frac{2L_{f2D}}{x_{e2D}}=0.05, 0.2, 0.4, 0.8 and 1$, with the fixed fracture height partially penetrating both layers under $M_2/M_1=10$.

Figure 12a provides pressure derivative behaviors under the above-mentioned fracture length to reservoir width ratios. For the one with the shortest fracture length ($\frac{2L_{f1D}}{x_{e1D}}=\frac{2L_{f2D}}{x_{e2D}}=0.05$), it shows that spherical flow is clearly evolved with −1/2 slope when 0.1 < t_D < 100, because the fracture vertical penetration and reservoir vertical thickness ratio is as small as 0.02. For production rate decline analysis, Figure 12b shows the modified dimensionless rate $q_{DM}$ versus modified dimensionless time $t_{DM}$. One of the advantages of presenting the rate profile in the manner of $q_{DM} vs t_{DM}$ is that by plotting it this way, the early- and late-time rate responses corresponding to various reservoir and fracture geometry can be systematically compared as shown in this figure, which is very useful and meaningful when matching field data.

Scenario 4: Effect of reservoir heterogeneity under various M₂/M₁ and C_S2/C_S1 for 2-layer commingled fractured vertical well without inter-layer pure-planar crossflow. In this scenario, the lower layer has larger reservoir length, width, and thickness, and the fracture extends horizontally twice longer in the lower layer than in the upper layer. Figure 13a presents pressure derivative behaviors of the fractured vertical well, producing two vertically stacked non-communicating layers under various mobility and storativity ratios from the proposed SSFM-ML model without inter-layer pure-planar crossflow. The solutions under the same 2-layer reservoir and well model by Kappa are provided in Figure 13a for comparison. The solutions from these two methods are consistent. It is noticeable that derivative behaviors for M₂/M₁ = 10 and 50 increase sharply after early time when the ratio of C_S2/C_S1 = 1, which is not seen for homogeneous case, and in particular a unit slope is clearly evolved for the case of M₂/M₁ = 50. This is because production is mainly provided by the upper layer (with better flowing properties) and the lower layer cannot transfer fluid to the upper layer due to the non-existence of inter-layer pure-planar crossflow. The commingled well acts as if it was effectively produced in a bounded reservoir with the same size as the upper layer. When the pressure drops near wellbore it becomes greater and greater, the lower layer starts to increase its production share, and pressure derivative behavior starts to deviate away from the PSS-like feature. In field practice, when testing data of wells produced from two layers with great distinction in flowing properties showing PSS-like features, careful analysis is required to confirm whether the response reflects the volumetric flow of both layers or if it mainly reflects the layer with better flowing properties.

Figure 13b shows the production behaviors of the 2-layer commingled fractured well in terms of $q_{DM} vs. t_{DM}$. Early and late time responses for all cases also hold together on their respective stems, similar to the previous scenario, in which inter-layer pure-planar crossflow exists. The time duration between the end of the linear flow stem and the start of fully developed PSS flow stem reflects the complex flow regime dynamics. More complex geometries of the reservoirs and the fracture systems cause variation during this period.

Scenario 5: Effect of storativity ($\mathit{\varphi }{\mathit{c}}_{\mathit{t}}$) among 3-layer reservoir with inter-layer pure-planar crossflow for MHFHW. Figure 14 presents the pressure derivatives of the MHFHW with six stages of fractures for the three situations with different storativity ratio combinations among the three layers, while mobility is assumed to be uniform throughout all three layers. Overall, Figure 14 shows that smaller storativity ratios in the upper two layers, relative to the bottom layer, yield a greater pressure drop along the entire life of the well. It also shows that early linear flow ends later for reservoir with greater storativity. It is because of the fact that a greater storativity leads to a smaller diffusivity, causing the pressure-transient to travel more slowly, thus also leading to the later arrival of the PSS flow. As described by Zhao and Thompson [10] and Spivey and Lee [34], pressure derivatives should approach 0.5 once the flow regime enters pseudo-radial flow regime regardless of storativity distribution, as long as mobility throughout the reservoir is homogeneous.

Figure 15a presents the corresponding rate profiles of the well and illustrates that a greater storativity supports a better production rate for the entire well life. By applying material balance time, the production rates from reservoirs with different storativity ratio combinations of the three layers all decline following a unit slope trend during BDF. Replotting the production behaviors of Figure 15a in the manner of $q_{DM} vs t_{DM}$ on Figure 15b, all curves again merge together during the early time linear flow regime and late time PSS flow regime, respectively. It also can be seen that $q_{DM} vs. t_{DM}$ curves almost overlap each other during the entire production, offering a valuable characteristic for the analysis of this type of reservoir setting. In fact, it is expected that if field data are accurate enough the changing trend of nearby layered reservoir properties can be further deduced, which may become necessary if quality data from a comparative set of field cases with multilayered reservoir structure are collected.

4. Conclusions and Recommendation

Based on the study, the following conclusions can be drawn.

SSFM was applied to analyze multilayered reservoirs produced by complex well configuration, and solutions with a high level of accuracy were achieved. The proposed analytical modeling strategy holds advantages in many aspects compared to commercial software—accurate solutions and computing efficiency, for example. The powerful modeling technology allows reservoir flow models to closely conform to reservoir characterization based on field information interpreted from well logging, sedimentology, and stratigraphy sources.

The new concept of pure-plane source defined in this work serves the purpose of accounting for the hydraulic fluid transfer for flow in porous media, thus helping to establish a consolidated general theoretical framework for the modeling methodology. The inter-layer fluid crossflow between adjacent layers is defined accordingly as inter-layer pure-planar crossflow to support a more rigorous description for the commonly encountered crossflow phenomena in multilayered reservoirs.

Theoretical solutions from the proposed SSFM-ML multilayered model show that vertical/horizontal well (line source) solution can be closely approached by pure-plane source with smaller vertical penetration and/or horizontal extension, being properly stretched to approximate wellbore diameter, and that vertical well or horizontal well intercepting multiple blocks can also be addressed by the proposed analytical multilayered model.

The newly proposed dimensionless modified rate and time, q_DM vs. t_DM, advances the innovative concept originally proposed by Zhao et al. for homogeneous 2-dimensional reservoir to 3-dimensional multilayered reservoir successfully. This methodology offers a unique way to reduce the uncertainty in matching system parameters. The application of this technique greatly helps organize standard production decline solutions and assist analyzing field production data for MHFHW producing from multi-layered reservoir in a systematically integrative manner.

Because of the unique features of the modified dimensionless type curves of q_DM vs. t_DM and its potential for analyzing and characterizing real reservoirs in field practice, it is highly recommended that future research should be conducted for more complicated two- and three-dimensional reservoirs using a similar strategy of adopting other new forms of modified dimensionless rate and time terms.