Pressure Transient Behavior of Horizontal Wells Intersecting Multiple Hydraulic Fractures in Naturally Fractured Reservoirs

In this study, we present a mesh-free semi-analytical technique for modeling pressure transient behavior of continuously and discretely hydraulically and naturally fractured reservoirs for a single-phase fluid. In our model, we consider a 3D reservoir, where each fracture is explicitly modeled without any upscaling or homogenization as required for dual-porosity media. Fractures can have finite or infinite conductivities, and the formation (matrix) is assumed to have a finite permeability. Our approach is based on the boundary element method. The method has advantages such as the absence of grids and reduced dimensionality. It provides continuous rather than discrete solutions. The uniform-pressure boundary condition over the wellbore is used in our mathematical model. This is the true physical boundary condition for any type of well, whether fractured or not, provided that the friction pressure drop in the wellbore is small and the fluid is Newtonian. The method is sufficiently general to be applied to many different well geometries and reservoir geological settings, where the spatial domain may include arbitrary fracture and/or fault distribution, a number of horizontal wells with and without hydraulic fractures, and different types of outer boundaries. The model also applies to multistage hydraulically fractured horizontal wells in homogenous reservoirs. More specifically, it is applied to investigate the pressure transient behavior of horizontal wells in continuously and discretely naturally fractured reservoirs, including multistage hydraulically fractured horizontal wells. A number of solutions have been published in the literature for horizontal wells in naturally fractured reservoirs using the conventional dual-porosity models that are not applicable to many of these reservoirs that contain horizontal wells with multiple fractures. Most published solutions for fractured horizontal wells in homogenous and naturally fractured reservoirs ignore the presence of the wellbore and the contribution to flow from the formation directly into the unfractured horizontal sections of the wellbore. Therefore, some of the flow regimes from these solutions are incorrect or do not exist, such as fracture-radial flow regime. In our solutions, all or some of multistage hydraulic fractures may intersect the natural fractures, which is very important for shale gas and oil reservoir production. The number and type of fractures (hydraulic or natural) intersecting the wellbore and with each other are not limited in both homogeneous and naturally fractured reservoirs. Our solutions are compared with a number of existing solutions published in the literature. Example diagnostic derivative plots are presented for a variety of horizontal wells with multiple fractures in homogenous and naturally fractured reservoirs.

natural) intersecting the wellbore and with each other are not limited in both homogeneous and naturally fractured reservoirs. Our solutions are compared with a number of existing solutions published in the literature. Example diagnostic derivative plots are presented for a variety of horizontal wells with multiple fractures in homogenous and naturally fractured reservoirs.
Keywords Naturally fractured reservoirs · Pressure transient well testing · Pressure transient behavior of horizontal wells · Multistage hydraulic fractures · Boundary element method

Introduction
Shale gas developments, particularly in the USA, have increased applications of horizontal wells with multistage hydraulic fracturing. Drilling horizontal wells to accelerate production and increase recovery efficiency of oil and gas are now common. With the increasing number of wells, pressure transient well tests have been conducted in different geological settings, particularly in carbonate and shale reservoirs. Natural fractures are common features of these reservoirs. In many naturally fractured carbonate reservoirs, faults often have high-permeability zones and are connected to many fractures with varying conductivities. In these naturally fractured reservoirs, faults and fractures can be discrete (i.e., not a connected fracture network system). In this paper, we present a number of solutions for horizontal wells with multiple fractures in naturally fractured and homogenous nonfractured (do not contain fractures) reservoirs. These solutions are applicable to finite-and infinite-conductivity transverse and longitudinal fractures in hydraulically and naturally fractured reservoirs, where the horizontal wellbore can intersect multiple hydraulic and natural fractures. In this paper, fractures denote both fractures and faults, unless differentiation is necessary.
Most of the published solutions on horizontal wells with multiple fractures in both fractured and nonfractured homogeneous reservoirs do not include the wellbore in the mathematical model. For transverse fractures, the effect of the horizontal wellbore cannot be ignored if the fracture conductivity is finite because the wellbore is perforated in single or multiple locations for the fracture initiation, as well as providing access for the fluid flow into the wellbore. In other words, perforation channels provide hydraulic communication between the wellbore, fractures, and formation if the well is not completed barefoot.
In this paper, the formation (matrix) is assumed to be transversely isotropic (k mh = k mx = k my ) and vertically anisotropic with horizontal and vertical permeabilities, k mh and k mv , respectively, whereas fractures are assumed to be isotropic. Each fracture may have a different permeability (k f ), length (l), and aperture (b). A nominal dimensionless fracture conductivity is defined as F D = k f b k mh l c , where l c is a nominal fracture half-length that is considered to be used as a reference or characteristic length, provided it is not much smaller or bigger than the mean fracture length. The dimensionless fracture conductivity for each fracture will be F Di . Al-Kobaisi et al. (2006) provided an extensive list of references and discussions concerning the pressure transient behavior of a single horizontal well intersecting multiple transverse vertical hydraulic fractures. They also discussed a number of flow regimes that can be observed during pressure transient tests. In general, the pressure transient behavior of a single horizontal well intersecting multiple transverse vertical hydraulic or natural fractures lies between two limits. First, for high fracture conductivities (F D ), the horizontal wellbore effect is negligible at early times, and the fractures dominate the behavior of the system in homogenous reservoirs. But at late times, there will be a significant contribution to flow from the formation directly into the horizontal sections of the wellbore. Second, as F D becomes smaller (F D ∞), the unfractured horizontal sections of the wellbore start to contribute to flow. For F D < 1, the flow contributions from multiple fractures become negligible, and the horizontal well dominates the behavior of the system; it is therefore extremely important to include the wellbore and the unfractured horizontal sections of the well. As we stated previously, most of the published solutions did not include the wellbore and the unfractured horizontal sections; consequently, some of the flow regimes observed from these solutions are incorrect, or simply they do not exist.
The inner boundary condition for a horizontal well with multiple transverse fractures is considered to be of the mixed type. In a mixed boundary value problem (MBVP), we have both the Neumann and Dirichlet boundary conditions imposed on the inner boundary wellbore surface. Furthermore, as stated above, except for openhole completions, hydraulic or natural fractures communicate with the wellbore through perforation tunnels, some of which may not be 100 % efficient (i.e., they are not cleaned up properly and totally debris free).
The majority of the published solutions do not specify where the wellbore pressure, which is uniform over the inner wall of the wellbore (sandface), is evaluated. This inner boundary condition over the wellbore is called the uniform-pressure condition. In addition to the uniform-pressure condition, the wellbore pressure is approximated at the middle point of the fracture-well intersection in the literature. In general, the middle point wellbore pressure approximation is crude and often incorrect when the pressure diffusion reaches any boundaries (i.e., top, bottom, faults, fractures, etc.); see Tables 2 and 3 of Biryukov and Kuchuk (2012b). For most single-fractured wellbore, the wellbore pressure is also approximated at an appropriate equivalent-pressure point (Gringarten and Ramey 1975). For some cases, particularly, horizontal wells and fractures, the wellbore pressure is also approximated by averaging over the flowing surface (Biryukov and Kuchuk 2012b;Kuchuk 1994;Ozkan and Raghavan 1991;Streltsova 1979;Wilkinson and Hammond 1990;Yildiz and Bassiouni 1990).
The published solutions for the pressure transient behavior of a single horizontal well intersecting multiple transverse vertical hydraulic or natural fractures can be divided into four categories as follows: 1. The first category solutions for multiple transverse fractures are generated using the superposition theorem in space and the pressure distribution solution for a single vertical fracture. Chen and Raghavan (1997), Guo et al. (1994), Horne and Temeng (1995), Raghavan et al. (1994), Zhou et al. (2013) used this approach to derive their solutions, in which the same flowing wellbore pressure is used for all fractures without a wellbore, and it is assumed that the uniform-flux or infinite-conductivity rectangular fractures fully penetrate the formation. For example, the dimensionless wellbore solution ( p w D ) for a reservoir with three vertical fractures (Horne and Temeng 1995) can be written as where ( p D ) 11 is the dimensionless pressure due to the first fracture, ( p D ) 12 is the dimensionless pressure due to the interference between the first and second fractures, and so on, p w D is the dimensionless flowing wellbore pressure, the normalized flow rates are (q D ) i = q i /q t , q i is the flow rate of the ith fracture, q t is the total flow rate, and i=1, 2, and 3. The definition of the dimensionless pressures are given below. Except for very high F D values, a few flow regimes are spurious due to these types of solutions. If a horizontal well intersects transverse finite-conductivity fractures, the first flow regime will be a fracture-radial or distorted fracture-radial flow regime. It will not be a fracture-linear or formation-linear flow regime. The first flow regime will be a formation-linear flow for an infinite-conductivity fracture transversely intersecting a horizontal well. In this formulation, there is no horizontal well contribution and the fractures do not intersect with a wellbore. An imaginary wellbore, in which p w D is the same, connects the fractures. 2. The second category solutions are similar to those described in the first category, but they include the contributions of the unfractured horizontal sections of the wellbore (Kuchuk and Habashy 1994;Raghavan et al. 1994). Raghavan et al. (1994) stated that their model also allows for multiple perforations along the horizontal well (in between fractures). The solution given in Eq.
(1) is also applicable for this case with a slight modification. As in the first category, except for very high F D values, a few flow regimes at early times will be incorrect and are artifacts of these types of solutions because multiple transverse fractures do not intersect the actual wellbore. 3. The third category solutions are hybrid ones, in which piecemeal Hegre 1991, 1994) or numerical/analytical (Al-Kobaisi et al. 2006) solution approaches are used. In these solutions, the wellbore is included. The flow in the fractures converges toward the wellbore, exhibiting a fracture-radial or distorted fracture-radial flow regime; see Al-Kobaisi et al. (2006). However, these solution do not include the contributions of the unfractured horizontal sections of the wellbore. As F D becomes smaller, the solutions do not change gradually into a horizontal well solution. This condition will be discussed in detail later. 4. The first numerical solution for the pressure transient behavior of a single horizontal well intersecting multiple transverse vertical hydraulic or natural fractures in the petroleum literature was presented by van Kruijsdijk and Dullaert (1989). After this paper, a number of numerical solutions were presented for horizontal wells intersecting multiple transverse vertical hydraulic or natural fractures or faults, including Akresh et al. (2004), Al-Thawad et al. (2001. The published solutions for horizontal wells in naturally fractured reservoirs have been based on the outdated Barenblatt et al. (1960) and Warren and Root (1963) dual-porosity models. In these solutions, a naturally fractured reservoirs is transformed into an equivalent (fictitious) homogeneous medium for fractures and an equivalent homogeneous matrix medium by use of the s → s f (s) transformation, where f (s) = ω(1−ω)s+λ (1−ω)s+λ , s is the Laplace transform variable, λ is the interporosity flow coefficient, and ω is the storativity ratio. As Braester (2005) stated, naturally fractured reservoirs are heterogeneous. It is therefore unrealistic to model naturally fractured reservoirs with only two ω and λ parameters, although there may be a few reservoirs that can be modeled using the s → s f (s) transformation. In the first paper published on this subject, de Carvalho and Rosa (1988) used the pseudosteady-state flow conditions in matrix blocks with the s → s f (s) transformation to derive the pressure transient solution for a horizontal well in a naturally fractured reservoir. After the de Carvalho and Rosa (1988) solution, many more solutions have been published in the petroleum literature using the → s f (s) transformation, while varying the matrix block shape and matrix block flow condition, and adding the interporosity skin factor for the matrix (Aguilera and Ng 1991;Du and Stewart 1992;Williams and Kikani 1990), in addition to the more recent papers by Abdulal et al. (2011), Brohi et al. (2011), Guo et al. (2012, Ketineni and Ertekin (2012), Lu et al. (2009), Medeiros et al. (2007, 2008, 2010, Nie et al. (2012a, b), Torcuk et al. (2013).
As reported by Kuchuk et al. (2014), if the permeabilities of the fracture segments are a few orders of magnitude larger than those of the matrix blocks, but not k m k f , then a macroscopic representative elementary volume can be constructed (i.e., dual-porosity models can therefore be used for continuously fractured reservoirs). If k m k f , then a macroscopic representative elementary volume (REV) cannot be constructed. Thus, analytical or numerical techniques should be used to model the pressure transient behavior of fractured reservoirs without fracture-segment homogenization or upscaling (i.e., fractures have to be modeled explicitly). In general, dual-porosity models should not be used for discretely fractured reservoirs (Kuchuk et al. 2014). Kuchuk et al. (2014) also showed that the resistance interface condition used by Barenblatt et al. (1960) and Warren and Root (1963) to specify fluid transport between the equivalent matrix and fractured media, dominates the dual-porosity model pressure transient behavior, irrespective of whether it is a vertical or horizontal well. The pseudosteady-state flow in the equivalent matrix medium that we observe in the Warren and Root (1963) model is a consequence of the resistance interface condition. This is not the actual pseudosteady-state flow that takes place after the transient flow in the equivalent matrix medium. The resistance interface condition for pseudosteady-state flow and the interporosity skin with transient flow in the matrix are of a nonphysical nature. The inclusion of these two nonphysical phenomena in the dual-porosity models introduces serious nonuniqueness problems for interpretation, significantly distorts flow regimes, and creates a naturally fractured reservoir look-alike pressure behavior but without any physical meaning of the parameters.
In a naturally fractured reservoir, a horizontal well usually intersects hundreds or thousands of fractures and several conductive and/or nonconductive faults [see Fig. 1.10 of the Nurmi et al. (1995) and Fig. 1 of the Kuchuk and Biryukov (2015) papers], particularly in carbonate formations (Al-Thawad et al. 2001. It was reported by Kuchuk and Biryukov (2015) and Kuchuk et al. (2014) that wellbore-intersecting fractures dominate the pressure transient behavior of both continuously and discretely naturally fractured reservoirs.
Next, we present the mesh-free semi-analytical technique for modeling the pressure transient behavior of hydraulically and naturally fractured horizontal wells in continuously and discretely 3D fractured and nonfractured reservoirs.

Mathematical Model
In this section, we present new semi-analytical solutions in the real-time domain for horizontal wells in reservoirs bounded at the top and bottom by the no-flow horizontal boundaries. Let us consider a single-phase slightly compressible fluid flow in an infinite reservoir with a no-flow condition on the top and bottom boundaries. The reservoir is naturally fractured and consists of fractures and the formation (matrix). In this model, the formation between the fractures becomes the matrix without predetermined shapes, and it can easily deal with fracture properties that exhibit power-law distributions or other type distributions, such as log-normal. The formation thickness h, porosity φ m , total isothermal compressibility (c t ) m , horizontal and vertical permeabilities k mh and k mv , and fluid viscosity μ are assumed to be constant, and time and pressure invariant. The fluid is assumed to be Newtonian. The fracture porosity φ f , total isothermal compressibility (c t ) f , and permeability k f are also assumed to be constant, and time and pressure invariant. The model also applies to multistage hydraulically fractured horizontal wells in homogenous (unconventional oil and gas) reservoirs.
The origin of the global Cartesian coordinates {x = 0, y = 0, z = 0} is in the center of the reservoir, and the top and bottom boundary planes are at z = ± h 2 , respectively. For this system, the pressure diffusion equation can be written as where P(x, y, z, t) = p 0 − p(x, y, z, t) is reservoir pressure change, and p o and p denote the initial and reservoir pressures, respectively. Let us consider that the reservoir is producing from a horizontal well with multiple production zones having radii r whi and lengths 2l whi , with their centers positioned at z chi from the bottom no-flow boundary (Fig. 1), where the subscripts c, w, h, and i denote center, wellbore, horizontal, and ith well (production section) respectively, i = 1, . . . , N h and N h denotes the number of producing horizontal well sections. Modeling horizontal wells is a complicated problem due to vertical anisotropy that implies the use of Mathieu functions, which are computationally expensive, the problem is therefore simplified by assuming that the well has a square cross section rather than circular one, while keeping the total production area and the length of the well the same. The effect of this assumption is negligible in the pressure response compared with a circular cross-sectional horizontal well. The assump- tion of a square cross-sectional wellbore is considered to be a fair trade-off in exchange for simplifying derivations. For a square cross-sectional horizontal well, the governing partial differential equation given by Eq. (2), and the initial condition given by Eq. (3) should be completed with the following conditions; see Fig. 1: where Υ wbhi is the flow rate entering the horizontal well from the reservoir. Because we are going to use different coordinate systems, it is assumed for the simplicity of notations that all functions are being expressed in the coordinate systems corresponding to the subscripts of their arguments (i.e., we omit the transformation signs from a local coordinate system Let us also consider a number of discrete fractures in the reservoir. For the sake of simplifying derivations, we assume that the fractures are vertical and have a rectangular shape with a half-length l f i and a half height h f i . Fluid flow in such a fracture, in its own coordinate system (x f i and z f i lie in a fracture plane and y f i is an outer normal; the origin lies in the center of the fracture at the distance z c f i from the bottom of the reservoir, see ∂P where, i = 1, 2, 3, . . . , N f , and N f is the total number of fractures in the reservoir, is the pressure derivative jump across the fracture, and Υ f i is the flow rate entering the fracture from the reservoir. Note that in Eq. (10), we assume an incompressible fluid flow in fractures (i.e., we do not consider any effects of fluid compressibility in the fracture). Biryukov and Kuchuk (2012a) showed (see Fig. 10 of their paper) that a compressible fluid flow in the fractures has almost no impact on wellbore pressure response. Equations (7)-(10) are only valid for an isolated fracture that is not intersecting anything else. When the ith fracture intersects a number of i k fractures, where k = 1, . . . , K i along the lines we should replace Eqs. (7)-(10) with the following: ⎠ , P where Eqs. (15) and (16) correspond respectively to pressure and the flux density continuity conditions at the intersection line, while Eq. (14) indicates the total flux conservation in the system of intersecting fractures. The term q f ii k corresponds to the flux density exchange between the ith and i k th fractures. Note that if the intersection line lies on the edge of the fracture, then the corresponding portion of the boundary condition in Eq. (12) is "overwritten" by Eqs. (15) and (16). When the horizontal well intersects a number of fractures, Eq. (6) is no longer valid, and then the following equation should be used instead: which accounts for additional influx terms Υ hii k due to flow from fractures into the wellbore. Fracture flow equations Eqs. (7)-(10) for the ith fracture intersecting with the horizontal well section j should be changed accordingly as where α i j is the angle between the horizontal well axis and the ith fracture plane, and x f i = ξ i j and z f i = ζ i j are the coordinates of the point of intersection of the well axis and the fracture plane in the ith fracture coordinate system. For simplicity, we assume that one fracture can intersect one horizontal well section only. If that is not the case, the fracture can be split into several fractures, each intersecting with only one horizontal well section. To simplify Eqs.
(2) through (25), the following dimensionless variables in the global coordinate system are defined as and the dimensionless variables in the local coordinate system are defined as where k f and l c denote the fracture permeability and the characteristic half-length or a reference length, respectively; l i denotes the half-length of the ith fracture; k mh , φ m , and (c t ) m denote the horizontal permeability, porosity, and total compressibility of the formation (matrix); h denotes the reservoir thickness; q hi denotes the flow rate into the ith horizontal section; q denotes a reference flow rate; μ denotes fluid viscosity; and t denotes time. As discussed previously, in fractured reservoirs, since there are many fractures of different lengths; we therefore choose a nominal half-length to be a reference or a characteristic length. In our solutions, we assume that the formation matrix can be anisotropic with the horizontal and vertical permeabilities, k mh and k mv , respectively, whereas the fracture media are assumed to be isotropic and each fracture may have a different permeability. For simplicity, we will omit the subscript D in the following sections because all derivations below will be done using the dimensionless variables given in Eqs. (26) and (27). The mathematical statement of the problem can now be written for the reservoir producing from horizontal well sections as ⎠ , p p Finally, the system should be completed by imposing conditions at each horizontal wellbore section on the flow rates q hj and pressure change p wbh j that are not independent. For example, if we assume that all horizontal production zones are part of one well and are hydraulically connected and produce with a constant flow rate q, we should add the following closure equations:

Transient Solution
In this section, we will describe a solution technique that is used for solving Eqs. (28)-(43). Let us first build an expansion for a flux density distribution on fractures. We consider a grid in the form where M i and N i are the number of grid cells in the x and z directions, respectively. Now, let us assume that the flux density can be approximated as where ω denotes the boundary element in the form Although we could use Chebyshev polynomials as we did in the 2D case (Biryukov and Kuchuk 2012a), we did not find them to be efficient for approximating the flux density distribution near the point of intersection with the horizontal well. On the other hand, with boundary elements ω, we can easily vary the local grid size in order to achieve the desired precision. Therefore, we decided to increase the grid density near the edges of the fracture and near the horizontal wellbore intersection (if there is one). The question of gridding will be discussed in detail in "Appendix 1." The flux density distribution defined by Eq. (44) will induce the following pressure change in the reservoir: where ⊗ denotes the time convolution operation, and The evaluation of Λ and the time integrals given in Eq. (46) are discussed in detail in "Appendix 1." We also consider a grid along the fracture intersection lines in the form = z 2 f ii k , and assume that the interfracture flow exchange terms can be expressed as where it is assumed that z f i (z n f i k i ) = z n f ii k . For the horizontal well, we define a grid in the form −1 = x 1 hi < x 2 hi < · · · < x n hi < · · · < x N hi +1 hi = 1 along the well and −1 = y 1 hi < y 2 hi < · · · < y m hi < · · · < y M hi /4+1 hi = 1 at z hi = ±1 and −1 = z 1 hi < z 2 hi < · · · < z m hi < · · · < z M hi /4+1 hi = 1 at y hi = ±1 in the transverse direction, where N hi and M hi are the number of grid cells in the x and z directions, respectively Let us also consider the following point source distribution on the wellbore: The pressure change in the reservoir induced by this source distribution can be expressed by the following equation as where with Ξ being defined by Eq. (48), and The total flow rate over the ith producing horizontal section due to the pressure change is obtained by substituting p hi from Eq. (55) into Eq. (41) as and Next, we seek the solution of Eqs. (28)-(43) in the following form: It can be observed that Eqs. (28)-(30) are satisfied automatically. The remaining equations cannot be satisfied everywhere because we have restricted ourselves to a limited number of expansion functions in p hi and p f i . We therefore choose a number of collocation points on every fracture and well segment (must correspond to the number of expansion functions) and consider these equations only on this set. We suggest taking a point in the middle of every grid cell to construct our expansion functions as x n hi = x n hi + x n+1 hi 2 , n = 1, . . . , N hi , on the ith horizontal wellbore section.
The quantities M f i , M hi , N f i , and N hi are the numbers of collocation points and/or grid cells in corresponding directions for fractures and horizontal wells, respectively. N f ii k is the number of collocation points and grid cells along the intersection line of ith and i k fractures. Their values should be chosen based upon the desired accuracy and the computational time preferences. This procedure is discussed in "Appendix 1." In the same manner, the collocation points on the fracture intersection lines can be defined as Considering Eqs. (30)-(43) only on these collocation points, we can obtain the following system of equations for the unknown variable coefficients c mn f i (t) and c mn hi (t): p f j (t) = p wbhk (t), if the jth fracture intersects the kth horizontall well, q n j j k (t) = q n j k j (t), if the jth fracture intersects the j k th, Equations (68) and (69) correspond to the pressure conditions that are prescribed on the well segments and fractures, i.e., to Eqs. (31) and (32), respectively. Equations (70) and (71) correspond to the flow rate conditions imposed on the horizontal well sections and fractures, i.e., to Eqs. (37) and (38), respectively. Equations (72) and (73) correspond to the pressure continuity condition on the fractured horizontal well section and the fracture-fracture intersection lines, i.e., to Eqs. (35) and (39), respectively. Finally, Eq. (74) indicates the flux continuity along the intersection line of two fractures (see Eq. 40), while Eqs. (75) and (76) correspond to the closer conditions imposed on the wellbore pressure and flow rates. Note that the equations corresponding to the fracture collocations points situated in the horizontal well cross section, as well as the coefficients c mn f i corresponding to the same grid cells should be removed. Also, for simplicity of notation, we used a superindex [i] = { f i, hi}, where N o = N h + N f is the total number of "objects" in the reservoir and where Φ(x, z, x 1 , x 2 , z 1 , z 2 , ν) is defined as the solution of the following problem: E(x, z, x 1 , x 2 , z 1 , z 2 , ν, ξ, ζ, α, β) is defined as the solution of the following problem: ϕ(x, z, x 0 , z 1 , z 2 , ν) is defined as the solution of the following problem: ν, ξ, ζ, α, β) is defined as the solution of the following problem: These mixed boundary value problems are time independent and relatively simple; therefore, they can be solved by a variety of methods. In "Appendix 2," we describe one possible approach based on the analytical elements method. Now, let us return to the system defined by Eqs. (68)-(76). To solve this system, we need to discretize it over time. The system contains the convolution operation, which is only convenient for discretizing on a fixed-step-size grid. At the same time, we would like to be able to increase the step size as time grows (due to the logarithmic behavior of pressure response). So, we suggest considering the following time grid: t s l = l × 2 s , l = 0, 1, 2, 3, 4 and s = s 1 , . . . , s 2 , which effectively consists of overlapping fixed-step subgrids (with steps of 2 s ). In our computations, we use s 1 = −35 and s 2 = 8, which covers the dimensionless time range from 2 −35 ≈ 2 × 10 −11 to 2 10 ≈ 1000, which is sufficient for the most cases. Let us also define a discrete convolution operation⊗ on this time grid as dτ . Now, the discretized version of the system defined by Eqs. (68)-(76) can be obtained by simply replacing the continuous convolution with one that is discrete, and the continuous time with discrete time values from the time grid. The values for the time points that are not in the grid can be obtained by interpolation.
The dimensionless wellbore pressure ( p w D ) with the wellbore storage (C D ) and skin (S) effects can be obtained using the superposition theorem that is given as an integro-differential equation of a convolution type, namely: where p D is the dimensionless constant-rate wellbore pressure (solutions are given above) and the dimensionless wellbore storage, which is defined as The above integro-differential equation can be easily solved numerically using a variety of methods. The most straightforward method consists of replacing the integral with its sum representation and derivatives with finite differences, which reduces the problem to a set of linear equations (Kucuk 1986).

Comparison of the Results
In these examples, we investigate the pressure transient behavior of a single horizontal well intersecting a number of vertical hydraulic and natural fractures in fractured and nonfractured reservoirs. The pressure solutions given previously will be compared with some well-known analytical solutions in the literature. We compare the pressure solution for a single horizontal well given by  with our solution. Figure 3 shows a comparison of the derivative of the two solutions for a single horizontal well. As shown in Fig. 3, our results (this study) compare very well with those of the  results.
The pressure transient behavior of a single finite-conductivity fracture has been well studied, and accurate solutions are presented both in graphical and tabular forms; see Cinco-Ley et al. (1978), Lee and Brockenbrough (1986). Figure 4 shows a comparison of our solution for a single vertical finite-conductivity fracture without a horizontal well with the t D F D = 0.2 π (Cinco) F D = π (Cinco) F D = 2 π (Cinco) F D = 10 π (Cinco) F D = 20 π (Cinco) F D = 0.2 π (this study) F D = π (this study) F D = 2 π (this study) F D = 10 π (this study) F D = 20 π (this study) Fig. 4 Comparison of dimensionless pressure for a single finite-conductivity vertical fracture with varying conductivities from Cinco-Ley et al. (1978) and this work  Riley et al. (2007) presented one of the most accurate solutions available in the petroleum literature. Riley et al. (2007) obtained the dimensionless pressure for a vertical fracture with finite conductivity, where the fracture is assumed to have an elliptical cross section, the flow within the fracture to be incompressible, and the reservoir to be infinite. Their solution is a line-source solution, from which they computed the wellbore pressure by using an equivalent wellbore radius. This is a reasonable approximation for the uniform-wellbore-pressure condition if the time is not very small. The Riley et al. (2007) results, as can be observed from Table 1 of Biryukov and Kuchuk (2012a), are approximately 1 % higher than the Biryukov and Kuchuk (2012a) results. This difference might be due to the use of an equivalent wellbore radius with a linesource solution by Riley et al. (2007), whereas Biryukov and Kuchuk (2012a) used a finite (actual) wellbore radius (in the dimensionless form r Dw = 10 −3 ) with the uniform-pressure boundary condition, which is the true boundary condition at the wellbore. For transverse fractures, the effect of the horizontal wellbore cannot be ignored if the well is either entirely completed barefoot (openhole) or cased and perforated over the entire wellbore. The flow contribution from the formation directly into the unfractured horizontal sections of the wellbore is negligible for high-conductivity fractures if the wellbore is cased and perforated in clusters over small sections, and hydraulically fractured. This is a common practice in horizontal wells in shale formations. The next example will illustrate the importance of inclusion of the wellbore in the solutions. As we stated previously, most published solutions on horizontal wells with multiple fractures in fractured and nonfractured (homogenous) reservoirs do not include the wellbore effect in the solution.
For this case, we will use the example given in one of the most recent articles published by Zhou et al. (2013). For this single-rectangular-fracture case, fracture, formation, and fluid  Table 1 of Zhou et al. (2013) and in Table 1 of this paper with a few additions: a horizontal well and F D = 10,000 instead of 20. The reservoir is homogenous and infinite, the center of the well is located at {0, 0, 0} (in the middle of the formation from the top and bottom), and the transverse vertical fracture intersects the 1000-ft horizontal at {0, 0, 0} and penetrates the entire formation. We will investigate the behavior of three cases: (1) a single vertical fracture (denoted without well in Fig. 5), (2) a single vertical fracture with the horizontal wellbore (denoted with wellbore) without the flow contribution from the formation into the horizontal section of the wellbore, and (3) a single vertical fracture with the wellbore and the flow contribution from the formation into the horizontal section of the wellbore (denoted with horizontal well).
As shown in Fig. 5, the single vertical fracture derivatives with and without wellbore are exactly the same because as discussed previously for high fracture conductivities (F D ), the horizontal wellbore effect without the noncontributing horizontal section is negligible at early times, and the fractures dominate the behavior of the system in homogenous reservoirs. As expected, the derivatives exhibit a formation-linear flow regime (m = 1/2) until about t D of 30, after which an infinite-acting pseudoradial flow regime is observed. The same figure also presents the derivative of the horizontal well without the fracture in a homogenous reservoir as a bench mark. The horizontal well derivative exhibits a first radial flow regime around the well and then an intermediate horizontal well linear flow regime from t D of 0.02 to 0.4 before a horizontal well infinite-acting pseudoradial flow regime. The derivative (denoted with horizontal well in Fig. 5) of the single vertical fracture with the contributing horizontal well section over its entire length also exhibits a formation-linear flow regime (m = 1/2) until about t D of 3 × 10 −5 , after which it deviates from the formation-linear flow regime until about t D of 0.02, after which the horizontal well dominates the derivative behavior, as shown by the "only horizontal well' derivative in Fig. 5. Figure 6 presents the dimensional derivatives to give an idea about the observability of these flow regimes. The formation-linear flow regime of the derivatives with and without wellbore last about 3 h before the start of an infinite-acting pseudoradial flow regime. The formation-linear flow regime of the derivative of the single vertical fracture with the contributing horizontal well lasts about 0.008 h (less than 1 min). Even without skin and wellbore storage effects, this flow regime would not be observable. Note that the derivative for this case exhibits a trilinear flow regime from 0.1 to 10 h, perhaps a look-alike one. It exhibits an intermediate horizontal well linear flow regime from 10 to 20 h. The infinite-acting pseudo-  Fig. 7 Comparison of pressure changes in a single vertical fracture with and without a wellbore and a horizontal well radial flow regime starts for all derivatives about 5000 h, which is unlikely to be achieved for a reasonable well test duration. Figure 7 presents the pressure changes for the three fracture cases and a single horizontal well. A very large pressure change difference is observed between the fracture with and without wellbore cases and the fracture with the horizontal well contribution case. Figure 8 presents the percentage differences for pressure changes and derivatives for the fracture without wellbore cases and the fracture with the horizontal well contribution case. As shown in these plots, more than one hundred percent difference in pressure changes and derivatives are observed. Note also that the difference in derivatives goes down after 10 h, and eventually will be small when all derivatives reach the infinite-acting pseudoradial flow regime at about 5000 h.
As shown here, the solutions (discussed above) without wellbore or the horizontal well contribution may have serious consequences for model identification and history matching, which may lead to serious inaccuracies in estimating well, fracture, and formation parameters from transient well test data. This large difference in the pressure changes particularly for long times, as shown in Fig. 8, will have a very large impact on the rate behavior of the reservoir.
We compare our solution with the one given by Al-Kobaisi et al. (2006) shown in Fig. 9. This horizontal well example is given in Fig. 16 of the (Ozkan 2014) paper. As we discussed unfractured horizontal sections of the wellbore, but, unlike most other solutions, it does include the wellbore in the fracture plane. Therefore, as can be observed in Fig. 9, the derivative curves first exhibit a fracture-radial flow regime (m = 0) and then a formation-linear flow regime (m = 1/2). Our results compare remarkably well with those of Al-Kobaisi et al. (2006) if we exclude the flow contributions from the unfractured horizontal sections of the wellbore. A slight difference can be observed between the derivatives of the two solutions due to our exact treatment of the uniform-wellbore-pressure condition. For low-finite-conductivity fractures, the duration of the fracture-radial flow regime becomes longer, and the formation-linear flow regime becomes shorter. The derivatives for F D = 1 do not exhibit a formation-linear flow regime. The derivatives for F D = 100 exhibit a very short fracture-radial flow regime, which cannot be observed in real pressure transient well tests. As we said above, our solutions normally include the flow contributions from the unfractured horizontal sections of the wellbore, but for this case we did not include them. Figure 10 compares our solution that includes the flow contributions from the unfractured horizontal sections of the wellbore with the Al- Kobaisi et al. (2006) solution given in Fig. 16 of their paper. We assume that the horizontal well length in this case is 200 ft, which is the same as the fracture length (see Table 1 of Al-Kobaisi et al. 2006) and is a very short well. As can  (2006) and our solutions diverge very rapidly for small F D values because the derivative behavior is dominated by the horizontal section of the wellbore. As we stated above, for large F D values, the derivatives compare with each other perfectly because the wellbore-intersecting fracture(s) dominate(s) the transient behavior of the system due to a very short well length. In other words, it is very clear from Fig. 10 that the unfractured horizontal sections dominate the transient behavior of the well for F D values less than 100. For F D = 1, we observe a horizontal well radial flow regime around the wellbore, but do not observe a formation-linear flow regime (m = 1/2). As can also observed from the figure, the difference between the derivatives with and without contributions from the unfractured horizontal sections of the wellbore for F D = 10, 100, and 1000 after t D of 0.01 are almost the same because of a very short well. As shown in Fig. 11, when the horizontal well length is increased from 200 to 1000 ft, the difference between the derivatives are clearly observable as in the example given in Fig. 5.
This example was also presented by Zhou et al. (2013) (see Fig. 14 of their paper) for the transient behavior of a horizontal wells intersected by six multiple hydraulic fractures. The fracture, formation, and fluid properties are given in Table 2 of Zhou et al. (2013) as: φ m = 0.1, (c t ) m = 3 × 10 −6 psi −1 , μ = 0.6 cp, k mh = k mv = 0.1 mD, h = 50 ft, r w = 0.354 ft, q = 63.65 RB/D, l f (x f ) = 300 ft, and F D = 30. Zhou et al. (2013) stated: "A case of a horizontal well with six stages of uniform fractures is simulated, as shown in Fig. 13. The thick red line represents the horizontal wellbore, and the blue diamond is the heel of the horizontal section." The horizontal well length obtained from Fig. 13 of their paper is 500 ft. They have not specified whether the wellbore is cased and perforated in clusters over the fractured sections, the entire horizontal wellbore is cased and perforated, or completed barefoot. We will investigate the behavior of two cases: (1) the wellbore is cased and perforated in clusters over the fractured sections and (2) the wellbore is cased and perforated over the entire horizontal well. Figure 12 shows a comparison of derivatives of six fractures that transversely intersect a horizontal well in an infinite reservoir for both cases: perforated in clusters over the fractured sections case is denoted "cluster" and the cased and perforated over the entire horizontal well case is denoted "entire wellbore." As expected, the derivatives of F D = 30, 100, 1000, and 10,000 exhibit a formation-linear (m = 1/2) flow regimes, after which the distorted formation pseudosteady-state flow (m ≈ 1) is observed because the formation volume between the fractures reaches a pseudosteady-state flow regime, which occurs after the formation-linear flow. In this case, we have included the derivative of F D = 30 because this was used by Zhou et al. (2013) for the six hydraulic fracture case. Finally, all derivatives approach the formation infinite-acting pseudoradial flow regime. The derivatives of F D = 1, and 10 exhibit a fracture-radial flow regime at early times, whereas the derivative of F D = 10 exhibit also a bilinear flow regime. As can be seen from this figure, the differences in the derivatives are not large because the fracture length (600 ft) is longer than the horizontal well length (500 ft). On the other hand, as shown in Fig. 13, which presents the percentage differences in pressure changes between two cases, the differences are large for small F D values; it is about 70 % for F D = 1 and more than 1 % for F D = 100. As we stated previously, the large differences have serious consequences for model identification and history matching. Figure 14 compares derivatives of six fracture solutions from Zhou et al. (2013) and this work for F D = 1, 10, and 30. The cases are shown in this figure: (1) the pressure and derivative data provided by Zhou (2014) (an imaginary wellbore, in which a common wellbore pressure is assumed to connect the fractures) are denoted "without-SPE 157367"; (2) the data from this work without wellbore are denoted "without-this study"; and (3) the wellbore is cased and perforated over the entire horizontal well is denoted "entire wellbore." As can be seen in Fig. 14, as in the single fracture case given above, the Zhou (2014) derivatives do not compare well with the derivates from this study. We observe a bilinear flow regime (m = 1/4) for F D = 1 instead of a fracture-radial flow regime. The differences in the derivates for F D = 10 is small but considerable for F D = 30, and the formation-linear flow regime is quite distorted. As in the single fracture case (Fig. 8) and in Fig. 13, the differences in pressure changes (not shown here) are considerably large. Let us consider a 1000-ft horizontal well in the x direction completed in the reservoir, where the center of the well is located at {0, 0, 0}, and the well is hydraulically fractured, as shown in Figs. 15 and 16. We will investigate the behavior of these three cases: (1) a horizontal well with only hydraulic fractures in a homogenous (nonfractured) reservoir; (2) a horizontal well in a naturally fractured reservoir without any hydraulic fractures; and (3) a horizontal well in a naturally fractured reservoir with 40 hydraulic fractures (Fig. 15) for F D = 0.01, 0.1, 1, 10, 100, 1000, and 10,0000. For the example of the horizontal well in a naturally fractured reservoir with hydraulic fracture, formation, and fluid properties are given in Table 2. In this system, each transverse hydraulic fracture intersects the horizontal well every 20 ft and extends 400 ft (the half fracture length l f ) in the y direction. The total number of hydraulic fractures in the x direction is 40. For each run, the dimensionless conductivity of 40 fractures remains the same, but is varied from 10 −2 to 10,0000 to compute dimensionless pressures and derivatives. For the first case,  the reservoir model is a horizontal well only, with hydraulic fractures in a homogenous (nonfractured) reservoir (Fig. 16). Figure 17 presents the dimensionless pressure derivatives for the first case for F D = 0.01, 0.1, 1, 10, 100, 1000 and 10,000. As can be seen from this figure, only F D = 1000 and 10,000 derivatives exhibit a short formation-linear (m = 1/2) flow regimes. The formation pseudosteady-state flow (m = 1) then proceeds. The formation pseudosteady-state flow regime takes place if the ratio of d/(2l f ) is small. The ratio for this case is 0.02. As shown in Fig. 17, the formation pseudosteady-state flow regime takes place after the formation-linear flow before the infinite-acting pseudoradial flow regime due to the whole formation. As F D becomes smaller, the effect of fractures on the derivatives becomes negligible, and the derivatives basically exhibit a horizontal well behavior without fractures. As shown in Fig. 17, the derivative for F D = 0.01 exhibits horizontal well flow regimes. If the entire horizontal wellbore is cased and perforated or completed barefoot, the very low-conductivity fracture derivatives will not exhibit fracture-radial flow regimes, but will exhibit horizontal well flow regimes. For instance, as shown in Fig. 17, we observe the first horizontal well radial flow regime . If the horizontal wellbore is cased and perforated in clusters over the fractured sections, the low-and moderate-conductivity fracture derivatives will exhibit fracture-radial flow regimes, but it will not exhibit horizontal well flow regimes. Figure 18 presents the dimensionless pressure derivatives for the second case (a horizontal well in a naturally fractured reservoir without hydraulic fractures) for F D = 0.01, 0.1, 1, 10, 100, 1000, and 10,000. As can be seen from this figure, F D = 10, 100, Dimensionless pressure derivatives of a horizontal well that intersects 40 hydraulic fractures in a naturally fractured reservoir for various F D values 1000, and 10,000 derivatives exhibit a very long formation-linear (m = 1/2) flow regime compared to the linear flow regime observed in Fig. 17. For this case, eight natural fractures with different lengths intersect the horizontal well and extend from the well center into the formation that is 4.5 times the well length. Therefore, the formation-linear flow regime is very long. The derivative F D = 10 exhibits a 1/4-slope bilinear flow regime at early times most likely after a fracture-radial flow regime. This radial flow regime cannot be the formationradial flow regime around the horizontal well because the derivative value (slope) is small. As with the first case, the effect of natural fractures on the derivatives becomes negligible, and the derivatives basically exhibit a horizontal well behavior without fractures for small F D values. It should be noted that a formation pseudosteady-state flow regime is absent for this case. Figure 19 presents the dimensionless pressure derivatives for the third case (a horizontal well in a naturally fractured reservoir with hydraulic fractures, as shown in Fig. 15) for F D = 0.01, 0.1, 1, 10, 100, 1000, and 10,000. As can be seen from this figure, the behavior in this case is very similar to that of the first case as shown in Fig. 17 because the wellboreintersecting 40 hydraulic fractures and eight natural fractures dominate the behavior of the system.

Conclusions
In this paper, we have presented solutions for the pressure transient behavior of horizontal wells intersected by multiple hydraulic and/or natural fractures in homogenous (nonfractured) and naturally fractured reservoirs. A mesh-free semi-analytical technique was used for obtaining pressure transient solutions to multistage fractured horizontal wells in continuously and discretely naturally fractured and homogenous reservoirs. Our solution technique is based on the boundary element method, which has advantages such as the absence of grids, reduced dimensionality, and the ability to obtain continuous rather than discrete solutions. This technique is sufficiently general and can be applied to many different reservoir geological settings and well geometries.
In our solutions, naturally fractured reservoirs can contain periodically or arbitrarily distributed finite-and/or infinite-conductivity fractures with different lengths and orientations. There are many factors, including fracture conductivities, lengths, and distributions as well as whether or not fractures intersect the wellbore, that dominate the pressure transient behavior of horizontal wells intersected by multiple hydraulic fractures in naturally fractured reservoirs.
Our solutions are compared with a number of existing solutions published in the literature. It has been shown that most of the published solutions ignore the contribution to flow from the formation directly into the unfractured horizontal sections of the wellbore. With one exception, the published solutions did not include the wellbore in the fracture plane. Some of the flow regimes from these solutions are therefore incorrect and superfluous. Our solution includes the effects of the wellbore in the fracture plane and the unfractured horizontal sections of the wellbore. In addition, our solutions are also given for horizontal wells with multistage hydraulic fractures in continuously and discretely naturally fractured reservoirs without using the no-fracture containing Warren and Root dual-porosity type models. In our solutions, all or some of the multistage hydraulic fractures may intersect the natural fractures, which is very important for shale gas and oil reservoir production.
It is shown that the solutions without the wellbore and the unfractured sections cannot capture the true early-time response, such as fracture-radial flow. More than one hundred percent difference in pressure changes and derivatives are observed among these solutions with and without the unfractured sections. This large difference will have serious consequences for model identification and history matching, which may lead to serious inaccuracies in estimating well, fracture, and formation parameters from transient well test data.
The exact treatment of the uniform-wellbore-pressure condition and the inclusion of the wellbore and the unfractured sections of the horizontal well result in the identification of new flow regimes that were not apparent from most of the published solutions.
Diagnostic derivative plots have been presented for a variety of horizontal wells with multiple fractures in homogenous and naturally fractured reservoirs. It has been shown that these reservoirs exhibit many different flow regimes.
where r wxi = r whj sin α i j l f i and r wzi = r whj h f i . Four cells contained inside of the well cross section and the corresponding coefficients c mn f i should be discarded from the systems in Eqs. (68)-(76). q 1 (x) = N n=1 a n T 2n Note that we used only even polynomials and the same expansion coefficients on opposite sides due to the symmetry of the problem. Now, Φ can be written as where S 2n (x, z, r ) = 1 2π π 0 ln r 2 (x − cos θ) 2 + y 2 cos(2nθ)dθ, F(x, y) = 1 2π ln((x − ξ) 2 + (z − ζ ) 2 )dξ dζ = 1 2π x y(ln(x 2 + y 2 ) − 3) + x 2 atan(y/x) + y 2 atan(x/y) .
Finally, similar technique can be used to find the solution of Eqs. (89)-(92) in the following form: