Shale Reservoir Drainage Visualized for a Wolfcamp Well (Midland Basin, West Texas, USA)

Abstract

Closed-form solution-methods were applied to visualize the flow near hydraulic fractures at high resolution. The results reveal that most fluid moves into the tips of the fractures. Stranded oil may occur between the fractures in stagnant flow zones (dead zones), which remain outside the drainage reach of the hydraulic main fractures, over the economic life of the typical well (30–40 years). Highly conductive micro-cracks would further improve recovery factors. The visualization of flow near hypothetical micro-cracks normal to the main fractures in a Wolfcamp well shows such micro-cracks support the recovery of hydrocarbons from deeper in the matrix, but still leave matrix portions un-drained with the concurrent fracture spacing of 60 ft (~18 m). Our study also suggests that the traditional way of studying reservoir depletion by mainly looking at pressure plots should, for hydraulically fractured shale reservoirs, be complemented with high resolution plots of the drainage pattern based on time integration of the velocity field.

1. Introduction

The shale revolution responsible for North America’s oil and gas renaissance has seen three waves. The first wave was the rapid expansion, between 2005 and 2010, of hydraulically fractured shale gas wells in the Barnett and other gassy shale plays (Haynesville, Fayetteville, Woodford, Eagle Ford, and Marcellus) with major growth currently confined to the Marcellus shale [1]. The second wave was a shift by North-American operators, starting in 2008, from gas to liquid-rich shale plays like the Eagle Ford and the Bakken [2,3]. The third wave has only just begun with the hydraulic fracturing of multi-stage wells now being successfully applied to re-develop thick shale formations in the Permian Basin. Our study focuses on the production performance of a multistage hydraulically fractures well in the Midland Basin, which comprises the eastern part of the Permian Basin.

With oil and gas prices in the epoch 2014–present being significantly lower than when the North American shale bonanza started in the previous decade [4], operators need innovative insights to improve recovery factors to gain more production income without cost increases. Currently, recovery factors are still a fraction of what industry achieved in conventional hydrocarbon reservoirs [5]. Reservoir and completion engineers are challenged to improve recovery factors of hydrocarbon resources in shale formations. Currently, such recovery factors are typically quite low, perhaps 20% in gas-rich shale reservoirs and less than 10% in liquid-rich plays, which compares poorly to those achieved in conventional reservoirs where ultimate recovery factors range between 40% for oil and 80% for gas [5]. The well stimulation process by hydraulic fracture treatment is a key tool to improve the recovery factors of shale gas and liquids.

Companies operating the unconventional shale plays, have made advances by adopting a manufacturing model, drilling multiple horizontal wells from a single well pad, completed with 100s of closely spaced hydraulic fractures per well. Next, developing new flow visualization methods to pinpoint where oil moves into the well via the hydraulic fractures has become very important to find ways to further improve the recovery of oil from shale. We developed a grid-less and meshless model tool based on Complex Analysis Methods (CAM), which allows rapid visualization of flow attributes around individual hydraulic fractures, using: (1) pressure field contour plots, (2) velocity field contour plots, (3) streamlines, and (4) drainage contours. The accuracy of the CAM tool has been benchmarked against independent reservoir simulators, which showed no discernable difference in results [6,7,8], except for the higher resolution of the CAM results. The method is suitable for modeling the drainage pattern near multiple hydraulic fractures and natural fractures. Figure 1a,b show an example of how natural fractures can completely redistribute the flow in a reservoir, dependent on the fracture characteristics [9].

Our present study highlights two new insights. The first new insight is that using pressure plots as a proxy for reservoir depletion zones can be misleading. Pressure plots traditionally used to show the depletion of hydrocarbon reservoirs generally suggest that large areas are drained by hydraulically fractured wells [10,11,12]. Our present study shows that velocity plots and drained area plots reveal much more detail about the hydrocarbon migration in shale reservoirs than pressure plots alone. Although pressure plots and the associated velocity field, which is linear to the pressure gradient, are controlled by the same reservoir parameters, velocity front tracking allows the construction of drainage contours, based on the time-of-flight concept, around the stimulation zone. Flow velocities change rapidly with spatial changes in the pressure gradients that drive the fluid flow and only a fraction of the moving fluid reaches the well, due to very slow flow in ultra-low permeability reservoirs (see velocity quantification later in this study).

The second new insight is that micro-cracks may affect the flow field and drain fluid from the matrix deeper than when such cracks would not exist. However, fracture diagnostics precludes the confirmation of the existence of such cracks. Due to the lack of high-resolution fracture diagnostics [13], the extent of the fracture network often is incompletely disclosed and the fracture conductivity of known fractures remains shrouded in uncertainty. Our flow visualization method can show both the matrix drainage by the main fractures and any deeper reach into the matrix when micro-cracks occur as offshoots from the main fractures. This study documents our results and the methodology used is succinctly outlined.

2. Prior Work

Traditionally, the visualization of the fluid migration into the hydraulic fractures of shale reservoirs makes use of general purpose reservoir simulators, primarily developed to model heterogeneities in conventional reservoirs represented by spatial changes in porosity and permeability. Such grid-based solution methods have been modified to account for flow through matrix and fractures of hydraulically fractured wells [14]. Extremely narrow grid sizing is required near fractures to resolve the flow, which makes finite-element solutions time-consuming and costly to produce. A brief review of common methods for modeling flow in hydraulically fractured reservoirs is given below.

2.1. Dual Porosity Models and Limitations

Modeling flow in fractured porous media, such as rocks and unconsolidated sands, are among the most challenging problems both in hydrogeology [15,16] and in reservoir engineering [17,18]. One approach prevalent in simulations of fractured reservoirs uses the dual-porosity model [19,20] and its numerous expansions, which consider a description for a continuously fractured (pseudo-fractured) reservoir, as opposed to descriptions of discretely fractured reservoirs. The reservoir is modeled by separating domains of connected fractures with 100% porosity and a certain assumed permeability with matrix pores acting as storage space. The flow from the matrix to the wellbore is clearly affected by the presence of fractures, which the dual-porosity model attempts to capture by a transfer function using a shape factor. Among the several shortcomings of the dual porosity model, most importantly is the lack of any explicit parameter accounting for fracture density. Transfer functions and shape factors were introduced to control the exchange of fluid between the matrix and the fractures [21,22,23]. Numerous shape factors have been suggested [24,25,26,27,28] after the early work [19,20]. For example, mass transfer from matrix to wellbore would be 25% more effective using a shape factor for a rhombic fracture system rather than a cubic fracture system [28]. Later shape factors included time-dependency of the pressure transfer function. However, all proposed modifications of the shape factor affect the diffusion equation such that a larger shape factor (in the denominator of Equation (2) in Warren and Root [20]) will result in steeper pressure gradients.

2.2. Discrete Fracture Network Models and Limitations

The benefits of high-resolution modeling of the detailed fracture patterns in shale by numerical methods include [29]: (1) accurate estimation of reserves and recoverable resources, (2) validation of analytical models, and (3) optimization of fracture and completion designs. The effect of natural fractures (positive/negative) on the stimulation of unconventional reservoirs has been well documented [30,31,32]. The introduction of constrained discrete fracture networks (DFN) [33] where a secondary variable was used to control the randomness of the DFN models did not solve multiple problems inherent to a discrete description of the natural fractures. Alfred et al. [34] describe additional issues of a DFN approach, which includes the major issue of upscaling [35] when transferring the fracture information to the continuum world commonly used in reservoir simulation and in our case, in geomechanical modeling. To honor structural geology concepts rather than incomplete statistics, the Continuous Fracture Modeling (CFM) approach was proposed [36,37,38].

Aimene and Nairn [39] and Aimene and Ouenes [40] introduced the Material Point Method (MPM) to facilitate the representation of the natural fractures in a multi-scale continuum problem. MPM combines Lagrangian and Eulerian descriptions of the material and in doing so overcomes the inherent drawbacks of methods relying on using only one perspective. This is achieved by discretizing the material into Lagrangian material points. These material points reside in an Eulerian mesh that is used to solve the equations of motion during a time step. An interpolation from the material points to the mesh nodes allows for this dual perspective to be used. The MPM is an example of a meshless method [41] as the mesh is not required for a full description of the material, it is only used for calculation purposes. The proposed MPM workflow [42,43] reduces uncertainties in fracture design and optimization through the use of half lengths derived from the validated geomechanical strain.

The limitations of local grid refinements (LGR) to capture the pressure field near fractures in DFN are well known. One drawback is that LGR is computational intensive, especially when using structured Cartesian grids [44]. Another shortcoming is that such Cartesian grids cannot effectively represent complex and non-planar fracture networks. The latter can be better captured using unstructured grids to define Voronoi cells [29,45,46]. Some gains in computational time can be achieved using coarse scale parameters from numerical homogenization to upscale the non-transient regions in the reservoir model [47]. However, simulations with gridded and meshed finite elements remain computationally intensive, especially for multi-stage wells, which may have several hundred perforation zones in a single well. Further, reduction of computation time is typically achieved by applying so-called stencils to represent a multiple fractured well system by collocation of repetitive flow cells around a single or a small subset of the hydraulic fractures. The justification of such simplification typically were based on older completion techniques with well spacing of 120 m (394 ft), which reduced the effect of frac interference [29,44]. However, with today’s tight frac spacing reduced to 1/10th (or even 1/20th in 2018 fracture treatments) of the spacing used in the cited studies, frac interference effects occur early in the well-life and preclude the use of stencils. We agree with prior authors [48,49] that representing a multiple-fractured horizontal-well system by repetitive elements of a single fracture is not warranted.

An additional limitation of numerical codes for multi-stage fractured wells with long computation times is that achieving global optimization remains out of reach when multiple realizations of all potential fracture designs and well spacing cannot be computed due to practical time constraints. Numerous attempts exist to improve optimization routines, including, but not limited to, Bayesian optimization [50], parametric sensitivity analysis [51,52], genetic algorithms for derivative-free optimization [53,54], and particle swarms approach [55,56].

2.3. Analytical Models

Early analytical models of flow through rock fractures are mainly limited to cubic law formulations for fluid flow transport through just a single fracture [57], expanded in later work with factors accounting for wall roughness [58,59,60]. The transient effect of stress on flow in a single fracture due to dilation effects has also been quantified [61,62,63], with further efforts to account for inertia effects due to turbulent fluid at high rates of injection [64]. Recent studies account for solute transport [65,66], and chemical interaction [67,68], all confined to a single fracture model.

2.4. Drained Rock Volume

Prior studies have emphasized that the drainage region where reservoir fluid is moving due to the pressure gradient differs from the drained region, defined by the much smaller volume of fluid that has moved into the well after a certain time [69,70]. Closer study reveals the important distinction between stimulated rock volume (SRV) and fluid that actually reaches the well within practical field operation times, which we refer to as the drained rock volume (DRV). Pressure plots alone convey reservoir regions affected by flow, but are inconclusive to establish the DRV. Pressure plots do not reveal directly where fluid is stagnant. Pressure is a scalar quantity and velocity is a vector quantity. The difference is relevant, because the fact that fluid drainage rates are highest near the tips of hydraulic fractures is not quickly seen from pressure plots. Contouring of the velocities reveals that new insight immediately (see Section 4). Additionally, we can use velocities to track the drained volume using time-of-flight contours. The so-called depth of investigation increases as the pressure wave propagates from the well system deeper into the reservoir and closely corresponds to the outline of the SRV [71,72,73]. However, the areal extent of fluid drained may be much smaller than suggested by using pressure plots only. Using pressure depletion plots in combination with velocity field contour plots and time-of-flight based drainage contours, using the transient velocity field vectors to track individual fluid elements, is very useful for hydraulically fractured wells and arguably should be used as a new industry standard.

Exactly how much drainage occurs, and where the fluid produced in fractured wells comes from, can be visualized with our method at a high resolution [8]. Our present approach applies calculus from complex analysis to provide closed-form solutions for single phase flow involving an unlimited number of interacting fractures [9].

3. Methodology

Our modeling approach using CAM, side-steps the drawbacks of LGR and other gridded solutions (outlined in Section 2.2), by using a grid-less and meshless solution strategy, while still being able to model transient flow by time-stepping streamline integrals with very small time intervals to achieve smooth particle paths [74]. Single-phase flow in porous media can be modeled with analytical descriptions as has been advocated in many classical studies (e.g., Muskat [75,76]). The closed-form solution is commonly used to map curvilinear line-integrals to identify streamlines [77,78]. Our novel approach using complex analysis methods has been applied to fluid flooding and sweeping of hydrocarbon reservoirs with finite boundaries [6] and infinite boundaries [79,80]. Investigations of how hydraulic fractures drain the matrix in shale reservoirs are also possible using our method for streamline simulation and fluid time-of-flight tracing based on complex analysis methods. Applications include the visualization of hydrocarbon migration toward hydraulically fractured wells [9], and fluid drainage near multi-fractured horizontal wells with fracture hits [8,74,81]. The key algorithms and key workflow steps used in our model method are explained below.

3.1. Complex Potential and Generic Velocity Potential

The complex potential Ω(z) links the stream function (ψ) and potential function ($\varphi$):

$$\mathsf{\Omega }(z)=\varphi (x,y)+i\psi (x,y) ,\hspace{1em}{(\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(1)

The velocity components (v_x and v_y) can be obtained using $v_x=\partial \varphi /\partial x=\partial \psi /\partial y$ and $v_y=\partial \varphi /\partial y=-\partial \psi /\partial x$. The full velocity field V(z) is given by the velocity potential:

$$V(z)=\frac{d\mathsf{\Omega }(z)}{dz}=v_x-i{v}_y,\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1}).$$

(2)

The fluid flow is assumed irrotational, incompressible, immiscible and isothermal, such that the viscosity and density of the reservoir fluid remains constant. Capillary pressure, wettability and relative permeability effects are neglected. Details on complex analysis and flow applications can be found elsewhere [82,83,84,85,86,87].

3.2. Single Point Source

The interval source and dipole elements used in our study are based on the superposition of point sources and point sinks. The complex potential centered of a source/sink flow centered at the origin with time-dependent strength m(t), is:

$$\mathsf{\Omega }(z,t)=\frac{m(t)}{2\pi }\log (z),\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1}).$$

(3)

The flow in a porous medium due to a pressure gradient near a vertical injector wellbore in the planar (x,y) field is modeled by a point source (Figure 2a). Likewise, flow toward a vertical producer wellbore is represented by a point sink (Figure 2b). The corresponding velocity field is given by [79,80]:

$$V(z,t)=\frac{m(t)}{2\pi (z-z_c)},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1}),$$

(4)

where m(t) is the well strength (producers m(t) > 0; injectors m(t) < 0) in wells located at z_c.

Strength m(t) for a well with productivity q(t) (m³·s⁻¹), reservoir height h (m), and reservoir porosity n is:

$$m(t)=\frac{Bq(t)}{hn},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1}),$$

(5)

Factor B is the volume factor defined by the ratio of any oil volumes under reservoir conditions divided by that oil volume in the stock tank at surface conditions.

Figure 3 gives streamlines and time-of-flight contours for a homogenous horizontal reservoir due to fluid injection from a single vertical well with constant injection rate. Radial flow results in isochronous time-of-flight contours that are circles.

3.3. Interval Source

The hydraulic fractures in our model were modeled by interval-sources, obtained by superposing an infinite number of point sources (Figure 4). Assuming an interval source centred at the origin, with total length L and strength m(t), the complex potential is [7]:

$$\mathsf{\Omega }(z,t)=\frac{m(t)}{2\pi L}\left\{(z+0.5L)\log (z+0.5L)-(z-0.5L)\log (z-0.5L)\right\},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$$

(6)

Rotation of the interval-source by angle β and shifting its centre from the origin to coordinate z_c (Figure 5) yields:

$$\begin{array}{ll}\mathsf{\Omega }(z,t)& =\frac{m(t)}{2\pi L}\{(e^{-i\beta }(z-z_c)+0.5L)\log (e^{-i\beta }(z-z_c)+0.5L)\\ & -(e^{-i\beta }(z-z_c)-0.5L)\log (e^{-i\beta }(z-z_c)-0.5L)\},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1})\end{array}$$

(7)

The complex potential for N interval-sources, each with their own angle β_k, centre z_ck, length L_k and strength m_k, reads [8]:

$$\begin{array}{ll}\mathsf{\Omega }(z,t)& ={\displaystyle \sum _{k=1}^N\frac{m_k(t)}{2\pi L_k}\{(e^{-i{\beta }_k}(z-z_{ck})+0.5L_k)\log (e^{-i{\beta }_k}(z-z_{ck})+0.5L_k)}\\ & -(e^{-i{\beta }_k}(z-z_{ck})-0.5L_k)\log (e^{-i{\beta }_k}(z-z_{ck})-0.5L_k)\},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1})\end{array}$$

(8)

Differentiation of the above expression with respect to z, results in the velocity field expression for N interval sources with arbitrary orientations [8]:

$$V(z,t)={\displaystyle \sum _{k=1}^N\frac{m_k(t)}{2\pi L_k}{e}^{-i{\beta }_k}\left\{\log (e^{-i{\beta }_k}(z-z_{ck})+0.5L_k)-\log (e^{-i{\beta }_k}(z-z_{ck})-0.5L_k)\right\}}, (\mathrm{m}·{\mathrm{s}}^{-1})$$

(9)

Streamlines and time-of-flight contours for the area drained by a vertical fracture are obtained with our streamline tracing algorithm and visualized in Figure 6. The isochronous time-of-flight contours are confocal ellipses.

3.4. Single Point Dipole

The micro-cracks in our study were modeled by line dipoles. The point dipole, also known as point doublet, singularity doublet and singularity dipole, is a common occurrence in amongst other electromagnetism [88], aerodynamics [89], and fluid mechanics [90]. Deriving the point dipole velocity field entails combining the velocity fields of a point source and point sink (Figure 7) in a limiting process. Placing a point sink and point source of equal but opposing strengths on the real axis at a distance of 2ε from each other, meaning their locations are ε and −ε respectively, yields the velocity field:

$$\begin{array}{ll}V(z,t)& =\frac{q(t)}{2\pi hn}\times \frac{1}{z+\epsilon }+\frac{-q(t)}{2\pi hn}\times \frac{1}{z-\epsilon }\\ & =\frac{-1}{2\pi hn}\times \frac{q(t)\times 2\epsilon }{(z-\epsilon )(z+\epsilon )},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1})\end{array}$$

(10)

The point dipole is then obtained by letting 2ε approach zero while the product q(t)∙2ε remains constant (meaning q(t) increases inversely proportional to distance 2ε). By defining this constant as υ(t) = q(t)∙2ε, we find that the velocity field for the point dipole, located at the origin, is:

$$V(z,t)=\frac{-1}{2\pi hn}\times \frac{\upsilon (t)}{z^2},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1})$$

(11)

Rotating the point dipole of expression (11) counter-clockwise by angle β and placing it at location z_k is possible by employing the conformal mapping f(z) = e^−βi(z − z_k) and evaluating V(f(z))·f’(z). This results in the velocity field formula:

$$V(f(z,t))·f^{\prime }(z,t)=\frac{-1}{2\pi hn}\times \frac{\upsilon (t)·e^{i\beta }}{{(z-z_k)}^2},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1}).$$

(12)

Note that υ is the strength of the point dipole and, as a result of the limiting process, its unit is (m⁴·s⁻¹). The complex potential of a point dipole is hence given by:

$$\mathsf{\Omega }(z,t)=\frac{1}{2\pi hn}\times \frac{\upsilon (t)·e^{i\beta }}{z-z_k},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1}).$$

(13)

3.5. Line-Dipole

The analytical point dipole in expression (12) has its source side pointing parallel to the real axis, to the right, if β = 0 (Figure 8). We derive a line-dipole, starting with β = 0 and its center at real axis coordinate x_c. The resulting velocity field is:

$$V(z,t)=\frac{-1}{2\pi hn}\times \frac{\upsilon (t)}{{(z-x_c)}^2},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1}).$$

(14)

Next combine j of the point dipoles, with a uniform strength υ(t):

$$V(z,t)={\displaystyle \sum _{k=0}^{j-1}\frac{1}{j}\times \frac{-1}{2\pi hn}\times \frac{\upsilon (t)}{{(z-x_k)}^2}},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1})$$

(15)

To transform this expression into a Riemann integral, the interval [−0.5L, 0.5L] is partitioned into j intervals, by defining the j + 1 points:

$${\widehat{x}}_k=-0.5L+\frac{L}{j}k,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le k\le j$$

(16)

Therefore, each point dipole is located at:

$$x_k=\frac{{\widehat{x}}_{k+1}-{\widehat{x}}_k}{2}=-0.5L+\frac{L}{j}k+\frac{L}{2j}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\le k\le j-1.$$

(17)

In order to turn expression (15) into a Riemann integral, the spacing between two point dipoles (Δx_k) needs to be defined:

$$\mathsf{\Delta }{x}_k=x_{k+1}-x_k=\frac{L}{j},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\le k\le j-2.$$

(18)

Multiplying expression (15) with L/L and splitting the sum then gives:

$$V(z,t)={\displaystyle \sum _{k=0}^{j-2}\frac{-1}{2\pi hnL}\times \frac{\upsilon (t)}{{(z-x_k)}^2}}\mathsf{\Delta }{x}_k+\frac{1}{j}\times \frac{-1}{2\pi hn}\times \frac{\upsilon (t)}{{(z-x_{j-1})}^2},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1})$$

(19)

By letting j increase without bound, the last term in expression (19) goes to zero and a Riemann integral is obtained because Δx_k goes to zero. The resulting integral and velocity field is:

$$V(z,t)={\displaystyle \underset{-0.5L}{\overset{0.5L}{\int }}\frac{-1}{2\pi hnL}\times \frac{\upsilon (t)}{{(z-x_k)}^2}d{x}_k}=\frac{-1}{2\pi hn}\times \frac{\upsilon (t)}{(z-0.5L)(z+0.5L)},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1})$$

(20)

Rotating this expression and shifting the center is again achieved with the conformal mapping f(z) = e^−βi(z − z_c) and by evaluating V(f(z))·f′(z). The velocity field for the line-dipole is therefore:

$$V(z,t)=\frac{-1}{2\pi hn}\times \frac{\upsilon (t)·e^{-i\beta }}{(e^{-i\beta }(z-z_c)-0.5L)(e^{-i\beta }(z-z_c)+0.5L)},\hspace{1em}(\mathrm{m}·{\mathrm{s}}^{-1}).$$

(21)

Expression (21) is similar to the one in Weijermars and Van Harmelen [79], though more explicitly formulated. The complex potential of a line-dipole is hence given by:

$$\mathsf{\Omega }(z,t)=\frac{-\upsilon (t)}{2\pi hnL}\left\{\log (e^{-i\beta }(\mathrm{z}-z_c)-0.5L)-\log (e^{-i\beta }(\mathrm{z}-z_c)+0.5L)\right\},\hspace{1em}({\mathrm{m}}^2·{\mathrm{s}}^{-1}).$$

(22)

3.6. History Matching and Flux Allocation

A general expression used for history matching of the prototype well rate is given by Duong [91] (using field units in bbls/day):

$$q_{well}(t)=q_{1}t(a,m)+q_{\infty },\hspace{1em}(\mathrm{stb}/\mathrm{day})$$

(23)

$$\mathrm{with} t(a,m)=t^{-m}{e}^{\frac{a}{(1-m)}(t^{1-m}-1)}$$

(24)

Cumulative oil production is:

$$N_p(t)=\frac{q_{well}}{a}{t}^m,\hspace{1em}(\mathrm{stb})$$

(25)

The next step is to allocate total well production rates to the individual fracs in our drainage model:

$$q_k(t)=C(1+WOR)q_{well}(t)(\frac{h_k{L}_k}{{\displaystyle {\sum }_k}{h}_k{L}_k}),\hspace{1em}({\mathrm{ft}}^3·{\mathrm{day}}^{-1})$$

(26)

C = 5.61458333 converts field input units of $q_{well}(t)$ in stb/day to field output units of $q_k(t)$ in ft³/day, and WOR = S_water/S_oil, the respective saturations of water and oil.

In our study, the total type curve output, q_well(t), was simply prorated to the stages studied (13 out of 33) and divided proportional to the surface area of each frac stage (as determined by the product of frac heights, h_k, and lengths L_k) relative to the total surface area of the 13 fracs, and correcting for the water to oil ratio, WOR:

$$q_k(t)=C(1+WOR)\times q_{welll}(t)\times \frac{13}{33}\times \frac{h_k\times L_k}{{\displaystyle {\sum }_k{h}_k\times L_k}},\hspace{1em}({\mathrm{ft}}^3·{\mathrm{day}}^{-1})$$

(27)

The prorated flux per frac over time was used in the drainage visualization of our study.

3.7. Oil in Place and Recovery Factor

The produced oil can be compared to the original oil in place:

$$\mathrm{OOIP}=\frac{7758Ahn(1-S_w)}{B},\hspace{1em}(stb)$$

(28)

With field units inputs used as follows: 7758 = barrels (bbls) of oil per (acre·ft); A = area in acres (TBD in acres); h = height in feet; n = porosity; S_w = water saturation; B = volume factor.

The oil recovery factor is:

F(t) = N_p(t)/OOIP

(29)

with F = Recovery factor, and N_p(t) the cumulative production (stb) till a certain time t.

3.8. Pressure Field

For a hydrostatic reservoir, the initial reservoir pressure P₀ is determined by the hydrostatic pressure gradient $\mathsf{\Delta }{P}_{Hydrostatic}/dL$ times reservoir depth (d):

$$P_{Hydrostatic}=(\mathsf{\Delta }{P}_{Hydrostatic}/dL)\times d,\hspace{1em}(\mathrm{Pa},\mathrm{psi})$$

(30)

The hydrostatic pressure gradient typically is $\mathsf{\Delta }{P}_{Hydrostatic}=0.45 \mathrm{psi}/\mathrm{ft}$. The pressure change due to production or injection is obtained from:

$$\mathsf{\Delta }P(t)=-real(\mathsf{\Omega }(t))\times \frac{\mu }{k},\hspace{1em}(\mathrm{Pa}, \mathrm{psi})$$

(31)

μ is the viscosity (psi·day or Pa·s, or Poise), k is the permeability (ft² or m² or Darcy).

The reservoir pressure consequent to the well flow will be:

$$P(t)=P_0+\mathsf{\Delta }P(t),\hspace{1em}(\mathrm{Pa},\mathrm{psi})$$

(32)

Each pressure field plot is generated by evaluating the pressure differential Equation (31) for a large number of complex-valued coordinates z, in conjunction with Matlab’s command ‘contourf’.

3.9. Streamline Tracing

The streamline tracing method uses a first order Eulerian displacement scheme. The initial particle position is z₀ (in complex coordinates) at time t₀ (unit in (s); see Weijermars et al. [7] for a non-dimensional approach). The velocity field V(z,t) needs to be solved to obtain the velocity vector components (v_x and v_y) for all positions of tracer elements at times t₁, denoted by z₁(t₁):

$$z_1(t_1)=z_0(t_0)+v(z_0(t_0))\times \mathsf{\Delta }t,\hspace{1em}(\mathrm{ft}).$$

(33)

The term v(z₀(t₀)) in expression (33) is the velocity of the traced particle at time t₀ and location z₀, for which the velocity field V(z,t) is used. The size of the time step is Δt. The locations of tracer elements at time t_j are:

$$z_j(t_j)=z_{j-1}(t_{j-1})+v(z_{j-1}(t_{j-1}))\times \mathsf{\Delta }t,\hspace{1em}(ft).$$

(34)

In Matlab we evaluate expression (34) by using the corresponding velocity field V(z,t) as follows:

$$z_j(t_j)=z_{j-1}(t_{j-1})+\left[real\left\{V(z_{j-1}(t_{j-1}))\right\}+1i\times imag\left\{-V(z_{j-1}(t_{j-1}))\right\}\right]\times \mathsf{\Delta }t,\hspace{1em}(ft).$$

(35)

In order to maintain smooth streamlines and time-of-flight contours (TOFCs), the time step Δt needs to be smaller when the strength of point sources/sinks increases. The velocity field plots in this study map the spatial variation in the magnitude of the fluid phase velocity, which is obtained by taking the absolute value |V(z,t)| at a given time t. Each velocity magnitude contour plot uses a large number of complex-valued coordinates z, in conjunction with Matlab’s command ‘contourf’.

4. Results

Our flow visualization method can show both the matrix drainage by the main fractures and any deeper reach into the matrix when micro-cracks occur as offshoots from the main fractures. First we show the drainage area when no micro-cracks develop. The locations of the main fractures for each stage of the multi-staged fractures in our prototype well were inferred from the micro-seismic image recorded by the operator during the frack job (Figure 9a). The images of Figure 9b–d show the reservoir response to drainage by a hydraulically fractured, horizontal well after 1 month of production. A summary of input data is given in

Table 1. Key parameters used in this study.

Hydraulic Fractures Parameters (Equation (27))					Duong Model: Total Well Strength Parameters (Equations (23–25))
Nr.	height hk (ft)	length Lk (ft)	azimuth	Center (ft)	$q_1$ = 9773 (stb/month)
1	48.7	898	122	190 + 960 × 1i	$q_{\infty }$ = 0 (stb/month)
2	53.2	1088	118	105 + 850 × 1i	$a$ = 1.133 (1/month)
3	63.6	1461	116	135 + 660 × 1i	$m$ = 1.2765
4	63.5	1134	119	175 + 445 × 1i
5	67.5	1100	127	70 + 305 × 1i	OOIP parameters (Equations (28) and (29)).
6	58.2	1223	123	60 + 145 × 1i	A = 320 acres
7	53.6	1295	118	10 − 5 × 1i	h = 57 ft
8	56.7	1328	124	−60 − 135 × 1i	n = 0.05
9	50.9	1248	116	−20 − 350 × 1i	Sw = 0.47
10	49.1	1171	116	−325 − 380 × 1i	B = 1.05
11	50.7	993	112	−390 − 560 × 1i
12	81.3	466	105	−170 − 850 × 1i	Reservoir parameters (for pressure scaling only, cf. Equation (31))
13	49	477	105	−170 − 1030 × 1i	Permeability $k$ = 1 μD
WOR = 1.1659					Viscosity $\mu$ = 1 cP

. The table includes a description of hydraulic fracture geometry and dimensions, the well productivity based on Duong parameters, original oil in place (OOIP) parameters, and other reservoir parameters.

The velocity plot (Figure 9b) reveals that the oil migrates fastest near the tips of the fractures, with an interior region between the fractures draining much slower. Classical analytical flow models [92] asserted that no fluid flow would occur from matrix into the fracture tips. Only flow of matrix normal to the fracture is assumed, which is clearly not what we see in our CAM simulations where flow near the fracture is fastest and not following sub-parallel or linear streamlines, but more complex particle paths. The simplification of linear flow would lead to overestimation of diffusion to the fracture, particularly when such fractures are short.

4.1. Velocity Field Versus Pressure Plots

The velocity field near hydraulic fractures is rarely visualized in commercial reservoir simulators [93,94]. Pressure depletion plots (Figure 9c) are commonly used as a proxy for drainage of the reservoir. The reservoir pressure when production starts is replaced by an evolving pressure gradient. The steepest pressure gradients occur at the onset of production, with local variations of higher and lower gradients. Figure 9c has the steepest pressure gradients between the red and blue zones (with tighter contour spacing), which broadly coincide with the region of higher flow velocities peripheral to the frac tips (Figure 9b). In commercial simulators, pressure plots resemble the analytical solution of Figure 9c but have limited resolution due to finite grid size and computational time limitations. Our model based on closed-form formulae provides infinite resolution and allows for the computation of a detailed pressure contour map for the central region between the fracs (Figure 9d). Due to relatively shallow pressure gradients, flow in the central zone is slowest. Contoured for smaller changes, pressure gradients still exist in the central region (Figure 9d), but pressure differences are more subtle than elsewhere in the fracked region. The shallow pressure gradient is the reason why the flow is nearly stagnant in the central region between the frac clusters (Figure 9b).

4.2. Implications for Production Efficiency

A prototype well that had 50 month historic production data was history matched with a decline curve using a Duong model [81,91,95]. The well type curve for Wolfcamp shale in the Midland Basin, West Texas, is given in Figure 10a. The production forecast was used for production allocation to individual fractures proportionally to each frac stage (see Section 3.6). The width of the drained matrix region around the fractures is visualized using the allocated well flux (Figure 10b). Even after 40 years of production, large portions of the matrix between the fractures remain un-drained. The drainage area and pressure plots were generated using a flow reversal principle: the produced fluid was injected back into the reservoir at the same rate as produced, which allows the construction of time-of-flight contours for the drained rock volume (Figure 10b). All pressure plots show positive pressure, but are equivalent to pressure depletion plots by taking the negative of the pressure scale.

4.3. Recovery Factors

Figure 10b shows the total drainage area for 40 years of production and the narrow region drained corresponds to the cumulative production given by the type curve (Figure 10a). OOIP of the prototype well was estimated for a well spacing of 320 acre/well resulting in a recovery factor of only 6% after 40 years, with 4% recovery already achieved after 5 years. Figure 10c plots the cumulative production and the corresponding recovery factor appreciation over a 40 year field life. These results are based on the history matched and drained rock volume (Figure 10b) obtained by our simulation method. Large undrained (dead) zones occur between the principal fracture zones. However, when micro-cracks are present to drain the deeper matrix, the dead zones may be narrower, which was further investigated.

4.4. Micro-Cracks

The CAM model can be adapted to account for drainage by micro-cracks, were these to occur, that reach deeper into the matrix. Such fractures are currently below the resolution of fracture diagnostics. Unknown matrix micro-cracks acting as hydraulic conduits may drain the matrix further away from the main fractures. Assuming micro-cracks may develop in hydraulically fractured shale reservoirs, the degree of micro-crack contribution to shale drainage will vary with the density and hydraulic conductivity of such cracks.

The flow support that hydraulic fractures receive from the microcracks was simulated using line dipoles (see Section 3.5). The strength of dipoles was systematically varied relative to the hydraulic fractures (modeled by interval sources) to show the effect on the drained region in our model. Examples of micro-crack assisted flow in the central region of the multistage fractured well of Figure 9 are shown in Figure 11. What the dipole fracs drain around the micro-crack tips as a consequence of their conductivity being higher than that of the matrix is transported toward the main hydraulic fractures. The depth of the drained area and length of the dipoles representing the micro-cracks is determined by the assigned strength of the dipole relative to that of the interval source used to model the hydraulic fracs. The combined dipole strength of micro-cracks that are connected to a particular hydraulic fracture is variable: equal to the hydraulic fracture strength (Figure 11a), five times the hydraulic fracture strength (Figure 11b), and 10 times the hydraulic fracture strength (Figure 11c).

The pressure plots (Figure 11a–c, row 1) show pressure contours begin to bungle near the micro-cracks when the hydraulic conductivity is larger. Such relatively steep pressure gradients are absent in Figure 9d, which is free of micro-cracks. The flow velocities are clearly enhanced around the micro-cracks (Figure 11a–c, row 2), due to the steep pressure gradients induced near the cracks (Figure 11a–c, row 1). As a result, the micro-cracks would extend the drainage region deeper into the matrix (Figure 11) as compared to situations where micro-cracks are absent (Figure 10b). The velocity peaks are confined to the immediate vicinity of the micro-cracks (Figure 11b,c, row 2), which concur with the locations were pressure gradient are steepest and pressure contours bungle (Figure 11b,c, row 1). The consequent drained region for each case is visualized in Figure 12. The higher resolution images show locally enhanced pressure gradients (Figure 11, row 1) and flow velocity peaks (Figure 11, row 2) occur near the micro-cracks, which facilitate the flow from the deeper matrix regions toward the main fractures.

5. Discussion

Our study suggests that matrix drainage by hydraulic fractures can be accelerated by micro-crack systems. Micro-crack tips drain deeper into the matrix as a consequence of their conductivity being higher than that of the matrix and thus may act as hydraulic fairways. The results show that micro-cracks will improve the drainage depth, but the existence of such micro-cracks is beyond the resolution of current fracture diagnostics. Although our method can visualize flow around the main fractures and any micro-cracks at high resolution, uncertainty remains about the precise location and geometric complexity of the micro-cracks. The relative strengths of dipoles, in our model, may be systematically varied relative to the interval source to show the effect of the micro-crack drainage with different hydraulic conductivity. Darcy’s law diffusion based flow descriptions was used to characterize the drainage of the SRV. This assumption implies nano-scale effects such as osmosis and capillary effects do not invalidate the macroscopic flow description. Future work should account for such nano-scale effects. In any case, micro-cracks would enhance and deepen the matrix area drained by the main hydraulic fractures.

Micro-cracks incorporated in the simulation method are synthetic and hypothetical. Such micro-cracks may contribute to the productivity of the well. However, any net gain in productivity is not accounted for in our model, because the net flux is history matched against the actual production, which fixes the 40 year production curve. What the model shows is that the drained rock volume, in the absence of micro-cracks, remains confined to the immediate vicinity of the hydraulic main fractures. Consequently, we are not underestimating well productivity or recovery factors. The main issue is that the region that will be drained may shift to the interior of the matrix (away from the main fractures) due the presence of micro-cracks.

Several simplifying assumptions are acknowledged. Multiphase flow is not included in our current model, for Wolfcamp B wells, reservoir pressures do not stay above the bubble-point for the lifetime of a well. However, we think the effect on well productivity in the oil window is immaterial due to the low gas content of the wells and late life well productivity contributes little to the cumulative production of the well. Although multiphase flow is not included in our current model, bubble-point effects can certainly be incorporated (subject of new work).

Our conclusion that the highest gradients will occur around the fracture tips, may be nuanced when proppant density varies and premature screen out occurs at fracture tips during the fracture treatment. When the flux around the tip is limited by the connectivity from the tip of the fracture to the wellbore, the flux will shift inward, which is accounted for in a forthcoming paper from our research group [96].

6. Conclusions

We conclude that even when micro-cracks develop minor splays into the matrix, the region between the larger master fractures will remain largely un-drained. The flow velocity and pressure gradients for both the main hydraulic fractures and micro-cracks will be largest near the fracture tips. The smallest velocities and flow stagnation occur in the zones between competing fractures, where the pressure gradients are shallower and become zero in pressure culmination points. The consequent oil entrapment between fractures would not occur in conventional, high-permeability reservoirs, because even slowly moving reservoir fluids would still reach the well on the time scale of the field life (5–50 years). However, unconventional reservoir rocks such as shale with permeabilities in the nano-Darcy range have ultra-low flow rates, even near the fracture tips, in the order of 10⁻⁷ m·s⁻¹, the rate at which finger nails grow.

The cluster spacing in the Wolfcamp well studied here is already tight (60 ft). Closer spacing of the fractures would accelerate early production but still creates flow stagnation points between the fractures. However, the un-drained oil occurring between the closely spaced fractures (in low permeability reservoirs with ultra-low flow rates) may be recovered by the refracking of shale wells with perforations placed midway between the prior cluster spacing and fracture initiation points. The second generation of hydraulic fractures will tap into the earlier, stagnant flow regions. The refracs will improve the recovery factors and this may boost well economics accordingly.