Shale Reservoir Drainage Visualized for a Wolfcamp Well (Midland Basin, West Texas, USA)
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.
Ruud Weijermars and Arnaud van Harmelen
Harold Vance Department of Petroleum Engineering, Texas A&M University, College Station, TX 77843-3116, USA
Received: 7 May 2018 / Accepted: 25 June 2018 / Published: 26 June 2018
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.
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 . 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 , 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 . 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 . 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 .
Figure 1. Plan view of reservoir showing streamlines (blue) and time-of-flight contours (red) toward a central vertical well. (a): Single well with a single hydraulic fracture. (b): Single well with multiple natural fractures in reservoir. After Van Harmelen and Weijermars .
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 , 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.
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 . 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.
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 . 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 ) will result in steeper pressure gradients.
Discrete Fracture Network Models and Limitations
The benefits of high-resolution modeling of the detailed fracture patterns in shale by numerical methods include : (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)  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.  describe additional issues of a DFN approach, which includes the major issue of upscaling  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  and Aimene and Ouenes  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  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 .
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 .
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 , parametric sensitivity analysis [51,52], genetic algorithms for derivative-free optimization [53,54], and particle swarms approach [55,56].
Early analytical models of flow through rock fractures are mainly limited to cubic law formulations for fluid flow transport through just a single fracture , 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 . Recent studies account for solute transport [65,66], and chemical interaction [67,68], all confined to a single fracture model.
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 . Our present approach applies calculus from complex analysis to provide closed-form solutions for single phase flow involving an unlimited number of interacting fractures .