Impact on Drained Rock Volume (DRV) of Storativity and Enhanced Permeability in Naturally Fractured Reservoirs: Upscaled Field Case from Hydraulic Fracturing Test Site (HFTS), Wolfcamp Formation, Midland Basin, West Texas

Hydraulic fracturing for economic production from unconventional reservoirs is subject to many subsurface uncertainties. One such uncertainty is the impact of natural fractures in the vicinity of hydraulic fractures in the reservoir on flow and thus the actual drained rock volume (DRV). We delineate three fundamental processes by which natural fractures can impact flow. Two of these mechanisms are due to the possibility of natural fracture networks to possess (i) enhanced permeability and (ii) enhanced storativity. A systematic approach was used to model the effects of these two mechanisms on flow patterns and drained regions in the reservoir. A third mechanism by which natural fractures may impact reservoir flow is by the reactivation of natural fractures that become extensions of the hydraulic fracture network. The DRV for all three mechanisms can be modeled in flow simulations based on Complex Analysis Methods (CAM), which offer infinite resolution down to a micro-fracture scale, and is thus complementary to numerical simulation methods. In addition to synthetic models, reservoir and natural fracture data from the Hydraulic Fracturing Test Site (Wolfcamp Formation, Midland Basin) were used to determine the real-world impact of natural fractures on drainage patterns in the reservoir. The spatial location and variability in the DRV was more influenced by the natural fracture enhanced permeability than enhanced storativity (related to enhanced porosity). A Carman–Kozeny correlation was used to relate porosity and permeability in the natural fractures. Our study introduces a groundbreaking upscaling procedure for flows with a high number of natural fractures, by combining object-based and flow-based upscaling methods. A key insight is that channeling of flow through natural fractures left undrained areas in the matrix between the fractures. The flow models presented in this study can be implemented to make quick and informed decisions regarding where any undrained volume occurs, which can then be targeted for refracturing. With the method outlined in our study, one can determine the impact and influence of natural fracture sets on the actual drained volume and where the drainage is focused. The DRV analysis of naturally fractured reservoirs will help to better determine the optimum hydraulic fracture design and well spacing to achieve the most efficient recovery rates.

Beyond natural field examples, mineback and laboratory experiments by Huang and Kim [18] confirmed that man-made hydraulic fractures do not always propagate linearly perpendicularly to the direction of minimum far-field stress and are not always planar in shape. This evidence provides support for the idea that hydraulic fractures are complex branching structures and models that do not honor this geometry may give incorrect results. The assumption of planar hydraulic fractures brings with it the creation of stagnation zones between individual fractures due to flow interference [19,20]. Such dead zones are detrimental to reservoir drainage and affect well productivity. If these hydraulic fractures are not planar, then the influence of stagnation zones changes [21]. It is now obvious how crucial it is to accurately model fracture geometry to ensure the optimization of wells under production.
The present study looks at the flow interaction in a reservoir between natural fracture sets and hydraulic fractures, with an emphasis on how these fracture networks influence the development of the drained rock volume (DRV). A series of methodical simulations allows us to understand how the natural fractures impact the DRV evolution. The effect on the DRV of complex fracture branching has been modeled elsewhere [21]. The present study uses closed-form analytical solutions based on complex analysis methods (CAM) to model flow in a 2D model of both the natural and hydraulic fractures. In this method, the hydraulic fractures are represented as line sinks and the natural fractures are modeled as an infinite number of line doublets using a new algorithm [22,23]. The interaction of the natural fractures and the hydraulic fractures is modeled in CAM to determine the flow response and pressure changes in the reservoir. Based on these responses, Eulerian particle tracking can then quantify the impact of natural fractures and hydraulic fractures on the DRV. Insights generated from the models can be used to optimize well production and recovery factors in unconventional reservoirs. This paper is organized as follows. We first highlight how natural fractures can impact flow in the subsurface via three major mechanisms (Section 2). These mechanisms are expounded upon and the ways in which they can affect the drainage area in unconventional reservoirs are discussed. A review of the previous ways natural fractures were modeled (in classical numerical reservoir simulations) together with a brief analysis of the drawbacks of such methods when compared to CAM is next presented. Following this is a discussion about the relation of natural fracture porosity and permeability and the surprising lack of data on these crucial properties. The next section (Section 3) looks at the methodology of using CAM to model natural fractures and hydraulic fractures in 2D. The following section (Section 4) describes the additional methods used to incorporate and model the natural fracture mechanisms that impact subsurface flow in the CAM workflow. Once the methods used are explained, we then move to the results from various flow simulations (Section 5). For the results, modeling begins with simple representative elementary volume (REV) models to show the impact of natural fracture with altered porosity and permeability. Subsequently, the models are Beyond natural field examples, mineback and laboratory experiments by Huang and Kim [18] confirmed that man-made hydraulic fractures do not always propagate linearly perpendicularly to the direction of minimum far-field stress and are not always planar in shape. This evidence provides support for the idea that hydraulic fractures are complex branching structures and models that do not honor this geometry may give incorrect results. The assumption of planar hydraulic fractures brings with it the creation of stagnation zones between individual fractures due to flow interference [19,20]. Such dead zones are detrimental to reservoir drainage and affect well productivity. If these hydraulic fractures are not planar, then the influence of stagnation zones changes [21]. It is now obvious how crucial it is to accurately model fracture geometry to ensure the optimization of wells under production.
The present study looks at the flow interaction in a reservoir between natural fracture sets and hydraulic fractures, with an emphasis on how these fracture networks influence the development of the drained rock volume (DRV). A series of methodical simulations allows us to understand how the natural fractures impact the DRV evolution. The effect on the DRV of complex fracture branching has been modeled elsewhere [21]. The present study uses closed-form analytical solutions based on complex analysis methods (CAM) to model flow in a 2D model of both the natural and hydraulic fractures. In this method, the hydraulic fractures are represented as line sinks and the natural fractures are modeled as an infinite number of line doublets using a new algorithm [22,23]. The interaction of the natural fractures and the hydraulic fractures is modeled in CAM to determine the flow response and pressure changes in the reservoir. Based on these responses, Eulerian particle tracking can then quantify the impact of natural fractures and hydraulic fractures on the DRV. Insights generated from the models can be used to optimize well production and recovery factors in unconventional reservoirs. This paper is organized as follows. We first highlight how natural fractures can impact flow in the subsurface via three major mechanisms (Section 2). These mechanisms are expounded upon and the ways in which they can affect the drainage area in unconventional reservoirs are discussed. A review of the previous ways natural fractures were modeled (in classical numerical reservoir simulations) together with a brief analysis of the drawbacks of such methods when compared to CAM is next presented. Following this is a discussion about the relation of natural fracture porosity and permeability and the surprising lack of data on these crucial properties. The next section (Section 3) looks at the methodology of using CAM to model natural fractures and hydraulic fractures in 2D. The following section (Section 4) describes the additional methods used to incorporate and model the natural fracture mechanisms that impact subsurface flow in the CAM workflow. Once the methods used are explained, we then move to the results from various flow simulations (Section 5). For the results, modeling begins with simple representative elementary volume (REV) models to show the impact of natural fracture with altered porosity and permeability. Subsequently, the models are extended to synthetic cases of

Natural Fracture and Hydraulic Fracture Models
The importance of natural fractures and their possible effects on fluid flow and well production in unconventional reservoirs has been noted by many authors, the work of whom is outlined (Sections 2.1 and 2.2). We propose three major mechanisms by which natural fracture properties may impact flow and describe each mechanism in detail (Section 2.1). We also look at previous attempts to model natural fractures in reservoir simulations and the drawbacks of the prior attempts (Section 2.2). The relationship between natural fracture porosity and permeability is reviewed (Section 2.3). These properties can affect reservoir drainage patterns, which is why field data were obtained from an unconventional shale reservoir in West Texas (Section 2.3).

Natural Fracture and Hydraulic Fracture Interaction Mechanisms
Numerous authors have stated that the presence of natural fractures increase production in hydraulically fractured wells in unconventional reservoirs [2,24]. Such a broad statement neglects the intricacies in natural fracture morphology, distribution and its ability to impact production. Gale et al. [25] stated that "fracture systems in shales are heterogeneous; they can enhance or detract from producibility, augment or reduce rock strength and have the propensity to interact with hydraulic fracture stimulation". Of importance is whether the natural fractures are sealing or not, which depends on the degree of cementation in the natural fractures. For natural fractures lacking cement, with no natural proppant (as used in hydraulic fractures), a significant reduction in permeability is possible but will not result in complete closure and thus, permeability in the natural fracture would still be above that for the intact host rock [26]. Another factor to consider is the connectivity of the natural fracture system. Cross-cutting and abutting fracture systems of different ages may not be hydraulically connected, depending on the degree of sealing. Here, it is possible that hydraulic fracturing can be beneficial for the reactivation of these natural fracture systems, which may lead to natural fracture networks becoming connected to the hydraulic fractures for the first time. In this study, we modelled natural fracture systems with an enhanced conductivity, i.e., cementation is not a hindrance to the flow potential within the system.
Although some ambiguity remains on the true nature of natural fractures influence on well production, research using static, object-based permeability suggests that natural fractures would enhance well productivity [2]. Three major mechanisms for the increase in productivity due to natural fractures have been put forward by Weijermars and Khanal [23]. These three production enhancement mechanisms related to natural fractures involve (1) equivalent permeability enhancement, (2) storage effects, due to enhanced porosity in natural fractures, (3) connection of hydraulic to natural fractures. Each mechanism is further discussed below.
(1) Equivalent permeability enhancement: The presence of a natural fracture system open to flow (uncemented) with higher permeability than the matrix, would increase the equivalent permeability of the overall reservoir. This enhanced equivalent permeability will result in a corresponding higher flow rate towards the hydraulically fractured well increasing the well productivity. (2) Storage effects due to natural fracture enhanced porosity: Natural fracture porosity may differ from the matrix either on initial formation of the fracture or due to later dissolution of precipitated minerals in the fracture space [25]. Due to size dependent sealing patterns, larger natural fractures are believed to have greater porosity [27] and as such, porosity in natural fractures is thought to be underestimated in most models. A greater porosity in the natural fractures than in the matrix may affect the extent of the drained area because porosity is a major control on time of flight for particles traveling along streamlines [28]. If the porous fractures are more fluid-filled than the surrounding matrix, storage effects will affect the well productivity. Uncemented fractures with enhanced porosity will allow for storage of hydrocarbons that, when tapped by the hydraulic fractures, will flow readily towards the well. (3) Connection of hydraulic fractures to natural fractures: Hydraulic fractures will propagate preferentially along planes of weakness in the reservoir such as those created by natural fracture systems. If a hydraulic fracture reactivates and connects to the natural fracture system, this connection leads to the natural fractures essentially becoming a direct extension of the hydraulic fracture pressure sink. The connection of both fracture systems correspondingly increases the total fracture surface area that is in contact with the reservoir matrix and improves the production rate of such wells.

Natural Fracture Modeling
Numerous attempts have been made to properly model fractured reservoirs that can accurately account for flow in such fractured porous media. The earliest attempt was made by Warren and Root [29] by using the dual-porosity model. Irregular natural fractures were modeled by using homogenous matrix blocks that are separated by orthogonal uniform natural fractures with fluid communication between the isotropic and matrix blocks governed by the inter-porosity flow coefficient (λ) and fracture storage capacity ratio (ω). Starting with this model, Kazemi et al. [30] introduced modifications that allowed for multiphase flow as well as the introduction of a new matrix shape factor. In addition to this work, numerous other authors have tried to adapt the Warren and Root [29] model to account for changes in matrix block geometry with new methods moving from double-porosity models to triple-porosity models [31,32]. Drawbacks of dual and multi porosity-based fracture models are that discrete fractures are not included and actual fracture density is not accounted for. Dual-porosity models also do not account for the flow paths followed when the fluid exchange occurs between the matrix and fractures, which can thus lead to inaccurate modeling of complex flow behaviors and can result in the wrong calculation of pressure gradients [33]. Another method to model naturally fractured reservoirs has been the use of Discrete Fracture Networks (DFN). For this model, fluid flow in the medium is represented through a system of connected natural fractures embedded within the rock matrix. This technique was first introduced by Long et al. [34] and has evolved over the years and seen increased use to model flow in conventional and unconventional naturally fractured hydrocarbon reservoirs [35,36]. The DFN method is typically used when (1) simulations done on a small scale where fracture dominance would otherwise result in an invalid upscaled continuum approximation, (2) in simulations on a larger scale where fracture dominance is small and the upscaled continuum model with only the largest fractures accounted for is valid [37]. Drawbacks of DFN modeling comes from the lack of data for the detailed inputs needed for the model such as fracture orientation, length, aperture and transmissibility along the (natural) fractures. The use of field analogs in surface outcrops may help fill these data gaps but there is no consensus on how accurately these measurements from outcrops match the subsurface. To combat this downside, current modeling attempts use a stochastic approach based on probability density functions to determine parameters of interest. This stochastic realization method can be used to create multiple realizations of the natural fracture patterns with fracture lengths following a power-law distribution [38]. The DFN method is also computationally intensive (and therefore expensive) as it requires very fine grids, which is particularly the case for multi-stage wells in unconventional reservoirs with numerous perforation zones per well [33].

Natural Fracture Porosity and Permeability
The effect of natural fractures on fluid flow is highly dependent on the reservoir type. Four major naturally fractured reservoir types have been identified by Nelson [39] based on the extent that fractures have altered the reservoir characteristics. Type 1 reservoirs have natural fractures that provide the bulk of the reservoir storage capacity and permeability, and typically have very high Energies 2019, 12, 3852 6 of 36 natural fracture density. In Type 2 reservoirs, permeability is essentially provided by the fractures while the matrix is responsible for the bulk of porosity. For Type 3, the reservoir matrix has high permeability and porosity but the permeability is further enhanced by the natural fracture system and can result in very high flow rates. Type 4 naturally fractured reservoirs have fractures that provide no additional porosity or permeability enhancement due to the fractures being filled with impermeable minerals. Natural fractures in Type 4 reservoirs are actually detrimental to fluid flow as they create significant reservoir anisotropy, which acts as barriers to flow [40].
One could say that Nelson's classification is mostly valid for conventional reservoirs and less applicable to shale reservoirs. Unconventional shale reservoirs have the majority of porosity contained within the rock matrix, while hydraulic fracturing is needed to create high enough permeability pathways for economical production. The majority of shale reservoirs also exhibit a high degree of natural fracturing. Due to the described attributes, unconventional shale reservoirs can be considered to range between Type 1 and Type 2 classification of naturally fractured reservoirs, with an example of Type 2 being the Spraberry reservoir in West Texas [40]. The extent to which the natural fracture systems in shale reservoirs affect hydrocarbon production due to enhanced storage and permeability is yet unclearly defined and remains nebulous [25]. There is consensus that hydraulic fracture propagation needs to take into account the impact of natural fractures on this propagation [41] but the impact that natural fractures have independently on production is not well constrained [25]. This is because core observations tend to show cemented natural fractures giving lower permeability and porosity measurements. However, field tests indicate much higher values for both permeability and porosity of natural fractures. Soeder [42] stated "typical natural fractures that enhance reservoir permeability to the point of commercial production are probably not obvious lithological features, such as near-wellbore calcite mineralized joints". Description of natural fractures in the Barnett shale show completely cemented fractures before hydraulic fracturing that subsequently became open and might demonstrate stress sensitivity [43]. The cited evidence shows that there is a strong possibility of natural fracture systems with enhanced porosity and permeability in shale reservoirs potentially high enough to impact fluid flow, which is crucial to accurately capture in any flow models.
Important characteristics of natural fractures include fracture length, aperture, orientation, density, spacing, porosity and permeability. Values for most of these parameters are difficult to obtain from the subsurface. Outcrops can give some indication of fracture length, density and spacing there is reason to believe that limited outcrop data do not give a proper representation of subsurface features that lie deeper within the earth [25]. What we do know is that many shales exhibit a wide range of fracture sizes and properties. The larger the natural fracture, the greater the porosity because of size-dependent sealing patterns [27] and it is believed that underestimation of natural fracture porosity may have occurred (due to this phenomenon) in some case studies. A value of 2% or less for the porosity of a natural fracture system is considered typical, however, field data from the Monterey shale Formation using samples from highly fractured parts, have shown values as high as 6% for natural fracture porosity [44]. Studies conducted by Weber and Bakker [45] as well as Lee et al. [46] give values of 2% to 7% for natural fracture porosities of the Marcellus shale [25].
The present study makes use of detailed core descriptions from the Hydraulic Fracturing Test Site (HFTS) for accurate natural fracture property and distribution data for our field case model (Sections 5.3 and 5.4). These descriptions come from six cores from a slanted well that sampled the rock volume around a hydraulically fractured well. These cores were located in the Upper and Lower Wolfcamp formation and this data (type based in origin, dip and dip direction of the fractures) was previously used to visualize fracture orientation, types of fracture and perforation clusters by Shrivastava et al. [12]. We make use of this data for a more realistic representation of the natural fracture system present in the subsurface in our flow models to determine the impact of this system on the DRV and its implication for well productivity. An essential corollary of our model is the introduction of a new upscaling method for natural fractures, which reduces the number of fractures to the critical ones, while maintaining the same equivalent permeability as the prototype. The upscaled model still contains discrete fractures to reveal their key impact on the flow. The revolutionary upscaling method makes use of a combination of object-based and flow-based upscaling techniques (Appendix C).

CAM Solution for Hydraulic Fractures and Natural Fractures
This work uses the Complex Analysis Method (CAM) to model fluid flow in the reservoir. CAM is an analytical method that was originally limited to the use of integral solutions to model streamlines for steady state flow [47][48][49]. This approach was expanded upon to include the Eulerian particle tracking of time-dependent flows with applications to natural lava flows and other gravity flows [50,51]. The basic methodology to describe gravity flows was then applied to flow in subsurface reservoirs, properly accounting for permeability, porosity and fluid volume factors, while benchmarking tests of the CAM using numerical streamline simulations proved highly successful [33,52]. We propose the use of CAM because of it being a gridless meshless method that allows for infinite resolution at the fracture scale and also has faster computational times as compared to gridded numerical simulators.
For the hydraulic fractures that are connected to the horizontal wellbore, we use the concept of line sources and sinks. The complex potential along an interval [a,b] is given by Potter [53] with time dependent strength m(t) as follows: Differentiation of this equation with respect to z (the complex coordinate that represents a point in the reservoir space) gives Equation (2), which represents the corresponding velocity field: In general terms for multiple line interval sources (k) with time dependent strength m k (t) (see Appendix A the instantaneous velocity field at time t can be calculated from: For the line interval source solution, we are also able to determine the corresponding pressure depletion due to fluid withdrawal from the reservoir. The real part of the complex potential can then be evaluated to quantify the pressure change ∆P(z, t) at location z at a given time t: The potential function represented by φ(z,t) has its pressure scaled based on fluid viscosity µ and permeability k of the reservoir. The actual pressure can be calculated from the initial pressure (P o ) plus the pressure change: A new algorithm for the modeling of natural fractures was proposed by Van Harmelen and Weijermars [22] that makes use of a complex potential function created by the superposing of an infinite amount of line doublets and is: [ft 2 /day] (6) Similarly to the solution for a line source (Equations (1) and (2)), differentiation of the specific complex potential equation of Equation (6) yields Equation (7), which gives the instantaneous velocity field in the natural fractures at time t.
Here, υ(t)(ft 4 /day) is the strength of the natural fracture, which scales the permeability contrast with the matrix (as further detailed in Section 4.1). The height, width and length of the natural fracture are denoted by h, W and L (ft) respectively, n is porosity, γ is the tilt angle of the natural fracture as shown in Figure 2. The variables z a1 , z a2 , z b1 , and z b2 give the corner points of the natural fracture domain.  (1) and (2)), differentiation of the specific complex potential equation of Equation (6) yields Equation (7), which gives the instantaneous velocity field in the natural fractures at time t. 21 Here, υ(t)(ft 4 /day) is the strength of the natural fracture, which scales the permeability contrast with the matrix (as further detailed in Section 4.1). The height, width and length of the natural fracture are denoted by h, W and L (ft) respectively, n is porosity, γ is the tilt angle of the natural fracture as shown in Figure 2. The variables za1, za2, zb1, and zb2 give the corner points of the natural fracture domain. Figure 2. Natural fracture model. L and W are the length and width; zc is the center; za1, za2, zb1, and zb2 are the corners; β is the wall angles, while γ is the rotation angle of the natural fracture. The blue arrows give the direction of the flow [22].
As for boundary and initial conditions, CAM can be used to model both steady-state flow as well as transient flow as shown in our models. The initial REV models used in Section 5.1 to demonstrate the fundamental impacts of natural fractures on flow, use constant rate boundary conditions (using a constant far field flow of 2.5 ft/year). For the hydraulic fracture line sink models, we were able to introduce transient flow by the use of a declining flow strength based on the declining rate of the forecasted well production that is allocated back into each hydraulic fracture segment. All of our models were coded using the Matlab programming language.

Modeling of Natural Fracture Interaction Mechanisms
The major controls on fluid flow propagation in porous media are the porosity and permeability of the domain. For a naturally fractured reservoir, one may consider two domains for flow, the unfractured rock matrix and the natural fractures present within the reservoir. This assumes that the natural fractures are uncemented and allow for flow. For streamline simulations, the flow paths (FP) and time of flight (TOF) of fluids being transported in porous media due to pressure sources/sinks are calculated by the equation of motion, which is intrinsically dependent on porosity and permeability in the reservoir. A study by Zuo and Weijermars [28] led to the creation of two fundamental rules for FP and TOF in porous media. The first rule shows that an increase in permeability decreased the time of flight, and conversely, an increase in porosity increased the time Figure 2. Natural fracture model. L and W are the length and width; z c is the center; z a1 , z a2 , z b1 , and z b2 are the corners; β is the wall angles, while γ is the rotation angle of the natural fracture. The blue arrows give the direction of the flow [22].
As for boundary and initial conditions, CAM can be used to model both steady-state flow as well as transient flow as shown in our models. The initial REV models used in Section 5.1 to demonstrate the fundamental impacts of natural fractures on flow, use constant rate boundary conditions (using a constant far field flow of 2.5 ft/year). For the hydraulic fracture line sink models, we were able to introduce transient flow by the use of a declining flow strength based on the declining rate of the forecasted well production that is allocated back into each hydraulic fracture segment. All of our models were coded using the Matlab programming language.

Modeling of Natural Fracture Interaction Mechanisms
The major controls on fluid flow propagation in porous media are the porosity and permeability of the domain. For a naturally fractured reservoir, one may consider two domains for flow, the unfractured rock matrix and the natural fractures present within the reservoir. This assumes that the natural fractures are uncemented and allow for flow. For streamline simulations, the flow paths (FP) and time of flight (TOF) of fluids being transported in porous media due to pressure sources/sinks are calculated by the equation of motion, which is intrinsically dependent on porosity and permeability in the reservoir. A study by Zuo and Weijermars [28] led to the creation of two fundamental rules for FP and TOF in porous media. The first rule shows that an increase in permeability decreased the time of flight, and conversely, an increase in porosity increased the time of flight. The second rule states that Armed with the above rules, we now proceed to explain the three principal mechanisms by which natural fractures may impact fluid flow in the reservoir. Natural fractures may result in localized discrete changes in both permeability (Section 4.1) and porosity or storativity (Section 4.2) in the reservoir domain, creating a direct impact on reservoir drainage patterns and drained areas. The third possibility is the reactivation and connection of natural fractures to the hydraulic fracture network (Section 4.3), which functions as an extension of the hydraulic fracture pressure sink.

Equivalent Permeability Enhancement
This mechanism is due to the difference of permeability within the natural fracture and the surrounding rock matrix. In unconventional reservoirs, the natural fracture permeability (k f ) is typically greater than that of the rock permeability (k m ). Weijermars and Khanal [23] show via explicit derivations how the permeability ratio (R k ) directly impacts the strength of flow in natural fractures as follows: The fracture hydraulic conductivity (C f ) is determined by the product of its fracture aperture (w f ) and its permeability (k f ): From this conductivity, we are able to define and scale the strength of the natural fracture segment (υ f ) in terms of corresponding permeability contrast with the matrix as follows: The length dimensions for the natural fractures (h f -natural fracture height, L f -natural fracture length, w f -fracture aperture) are directly specified in the CAM models and matrix flow velocity (v m ) can be measured near the fracture in the simulation. By fixing the constituent parameters at time t the equation for R k thus becomes Thus, from the above equation, we can set the permeability ratio using an assigned strength in the natural fractures in our CAM model.
The most important aspect of the permeability enhancement mechanism is that natural fractures do not act as fluid sinks (see also detailed examples in Weijermars and Khanal, [23]). Mechanism 1 assumes the natural fractures are not connected to the hydraulic fractures (unlike mechanism 3). Instead the natural fractures act as zones of flow acceleration and preferentially drain matrix fluid further away from the well at the end of the highly conductive natural fractures rather than from the nearby lower permeability matrix. Change in permeability in the natural fractures impacts both streamline patterns as well as time of flight. This mechanism was thoroughly modeled and investigated by Weijermars and Khanal [23], using a variety of natural fracture parameters and readers are referred to this seminal work for further detail. Although prior studies (e.g., Aguilera, [2]) that use static object-based permeability scaling also give results that natural fractures can enhance well productivity, the method employed by Weijermars and Khanal is based on dynamic flow-based upscaling and is believed to be more accurate. Flow-based upscaling of permeability explicitly shows how for a fractured medium the equivalent permeability increases greatly when compared to similar porous media that are non-fractured. It is this overall increase in equivalent permeability (due to the enhanced permeability of the natural fractures) that leads to a higher flow rate towards the well and thus higher recovery during the economic life of such wells. In this study, we extend this work to investigate the implications of equivalent permeability enhancement due to natural fractures on DRV in conjunction with porosity changes in natural fractures. A new upscaling method for discrete fractures is given in Appendix C, allowing for simultaneous changes to fracture permeability and porosity, which has not been investigated previously.

Natural Fracture Storativity Effect
In addition to effecting localized permeability changes, natural fractures have the ability to alter porosity. Shale reservoirs tend to exhibit a wide range of fracture sizes. Due to the industries limited data of natural fracture porosity, the effects of this on flow alteration in the subsurface has not been previously studied in any detail. We present a set of high-resolution simulations with altered porosity in the natural fractures to quantify how this parameter affects drainage in the subsurface (see the results in Section 5).
As with the change in permeability, we are now able to define a porosity ratio (R n ) for the porosity change inside of a natural fracture (n f ) compared to the matrix porosity (n m ) surrounding it, given by the following equation: For the CAM analytical solution, natural fracture alignment can be defined in relation to the hydraulic fracture. Equation (7) assumes that the porosity across both the fracture zone and matrix remains the same. If we remove this assumption, based on the evidence presented on porosity differences in natural fractures when compared to the reservoir matrix, Equation (7) can be locally modified to take into account the altered natural fracture porosity as follows: This equation will now account for both the altered porosity and permeability within the natural fracture domain. As we manually define the boundaries of the natural fractures, the tracked particles that are displaced based on the time dependent strength of the flow in the reservoir will have velocity increased or decreased based on the porosity and permeability once the fluid particles enter the natural fracture domain. The trajectories of these particles are set by the permeability in the reservoir matrix and natural fractures [28]. Based on Rule 2 for flight paths and time of flight contours in porous media [28]. The time of flight will be slower in natural fractures with a higher porosity than the matrix (the streamline patterns will not be affected). Thus, for a hydraulically fractured well, the presence or absence of natural fractures with different porosities (that may be in situ porosity or increased porosity due to natural fracture reactivation) will affect how far the matrix is drained (i.e., the shape and location of the DRV will be affected), which knowledge is relevant for fracture treatment design and well spacing decisions.

Natural Fractures as Extension to Hydraulic Fracture Network
The third mechanism that may cause natural fractures to increase well productivity occurs when natural fractures become extensions of the hydraulic fracture system. This can lead to the creation of complex fracture networks, defined as non-planar, branching fracture geometries that are caused by either strong stress shadow effects or by the interactions with natural fractures [54]. Wu and Olsen [54] further state that the efforts to study interaction between natural fractures and hydraulic fractures have taken various forms of theoretical, experimental and numerical work. From this research, they propose three possibilities due to the intersection of natural fractures and man-made hydraulic fractures. The first possibility is that the created hydraulic fracture propagates along its original directions and crosses the natural fracture with no change in orientation. A second possibility is that the hydraulic fracture could be arrested by the natural fracture and then continue to propagate along the natural fracture to finally exit at the tip of the natural fracture. Deflection of the hydraulic fracture into the natural fracture, followed by re-initiating out of the natural fracture at a point of weakness is given as the third possibility [55]. No matter the propagation due to the interaction, the overall effect is that the natural fractures that intersect with the hydraulic fractures become extensions of the pressure sink imposed on the reservoir due to the connection of the fracture network to the wellbore. One way to model these interactions is via the use of fractal theory to replicate the branching fracture geometry that can then be modeled using CAM. The changes and implications on DRV due to such complex branching fracture networks was modeled in detail in Nandlal and Weijermars [21], and as such, is not the subject of focus in the current work. The reader is referred to the previous paper with the major conclusion being that the presence of a complex hydraulic fracture network will increase the initial production rate from the reservoir.

Results
Using the CAM approach, we investigated systematically the effects of porosity and permeability alterations within natural fractures on fluid flow using a range of model designs. The changes in these two crucial parameters were studied to determine the effect on the drainage area in the reservoir. Obviously, a proper understanding of the DRV development in naturally fractured reservoirs has implications for production from both conventional and unconventional oil and gas reservoirs.
We adjust the fracture strength and porosity ratio to determine the impact on drained areas in the reservoir. Modeling starts with a simple planar fracture with varying porosity ratios as well as different natural fracture configurations (Section 5.1). The effects of natural fracture storativity and enhanced permeability on DRV are demonstrated and proved (Section 5.2). We investigated the flow patterns near hydraulic fractures (modeled as line sinks using CAM), and how the presence of natural fractures and their corresponding porosity and permeability may change the DRV. It should be noted that the CAM models used in these flow simulations assume hydraulic fractures of infinite conductivity. The creation of hydraulic fracture via geo-mechanical simulations and corresponding coupling to flow simulation is outside the scope of this work.
The results in Sections 5.1 and 5.2 are for synthetic models, intended to systematically demonstrate the effects of natural fractures via the natural fracture interaction mechanisms explained in Section 4. The idealistic representative elementary volumes (REV) and simple fracture models assume porosity changes are independent of any permeability changes. In reality, this may not be true and there are many established correlations that relate increases in porosity with corresponding increases in permeability. We make use of field data for natural fractures to determine the DRV in an actual reservoir (Sections 5.3 and 5.4). Field data obtained from cores in the Hydraulic Fracturing Test Site (HFTS) as well as porosity-permeability correlations are used to determine the impact of natural fractures in the case study. By incorporating real data in our models, we can more accurately determine the impact of natural fractures on the DRV in the field. This is relevant to next proposed methods for optimization of recovery in both highly fractured unconventional and conventional reservoirs.

Representative Elementary Volume (REV) Models
To properly understand the effects on fluid flow, we started with the modeling of a simple representative elementary volumes (REV) that use a constant far field flow. A representative elementary volume (REV) is defined as a volume over which a measurement can be made that is representative of the whole. Using the REV allows for the understanding of the physics behind any changes in drainage patterns (before moving on to more complex situations). The first model provided is a base case which we use to compare all subsequent models. In this model (Figure 3) we show a reservoir space in 2D with five natural fractures represented by discrete elements that have the same porosity and permeability as the reservoir space. Using Eulerian particle tracking, we determined the flow path based on a constant far field flow. Flight paths (FP) are displayed in blue ( Figure 3, left image) with the corresponding time-of-flight contours (TOFC) shown in red ( Figure 3, right image). The base model represents a flow time of 30 years with each TOFC representing the fluid displacement after 3 years with reservoir porosity of 5%. Referring back to the two fundamental rules for FP and TOF [28], we observed that with no change in porosity and permeability in the natural fractures, the FP and TOFC remain constant. changes in drainage patterns (before moving on to more complex situations). The first model provided is a base case which we use to compare all subsequent models. In this model (Figure 3) we show a reservoir space in 2D with five natural fractures represented by discrete elements that have the same porosity and permeability as the reservoir space. Using Eulerian particle tracking, we determined the flow path based on a constant far field flow. Flight paths (FP) are displayed in blue ( Figure 3, left image) with the corresponding time-of-flight contours (TOFC) shown in red ( Figure 3, right image). The base model represents a flow time of 30 years with each TOFC representing the fluid displacement after 3 years with reservoir porosity of 5%. Referring back to the two fundamental rules for FP and TOF [28], we observed that with no change in porosity and permeability in the natural fractures, the FP and TOFC remain constant. Porosity effects: The next REV model ( Figure 4) highlights the effect of systematically increasing the porosity in the natural fractures. As stated previously, fracture system porosities of 2% or less are considered typical [44], but values as high as 7% for natural fracture porosity in shale formations have been reported [25,46]. With numbers still based on very limited datasets, it is possible that porosity changes in natural fractures can be higher than the values reported thus far. Therefore, we model porosity changes by up to 15% to observe the impact on flow. The initial models decouple the correlation between increased porosity and permeability and such that there is no permeability change in the natural fracture relative to the matrix. When we use the term porosity, we mean connected porosity. Figure 4 shows the effects of increasing NF porosity on the FP and TOFC in the reservoir space. The reservoir porosity is kept at 5%, while NF porosity changes incrementally from being equal to the reservoir space to a high of 15%. The results clearly show that the change in porosity within the natural fracture has no effect on the streamline flow paths but does affect the time-of-flight contours. In NF 1, the porosity is the same as the reservoir and as such, there is no slowdown in the TOFC. From NF 2 to NF 5, we progressively increased the porosity to 6%, 8%, 10% and finally, 15%. The model shows that for each successive porosity increase in the natural fractures, the FP stays constant but the TOF increases. As we are using a constant run time for all models, the increase in TOF results in flow not reaching as far into the reservoir space for the natural fractured with higher porosity. With no porosity change, flow reaches out to approximately 75 ft in the reservoir space. With a porosity change from 5% to 15% in the natural fractures, flow is retarded and reaches only approximately 44 ft out into the reservoir space. Thus, a 10% increase in natural fracture porosity results in a 40% reduction in lateral flow extent. This result can have great implications for accurately determining the DRV in the subsurface when the reservoir rock has a high density of natural fractures with variable porosity. Porosity effects: The next REV model ( Figure 4) highlights the effect of systematically increasing the porosity in the natural fractures. As stated previously, fracture system porosities of 2% or less are considered typical [44], but values as high as 7% for natural fracture porosity in shale formations have been reported [25,46]. With numbers still based on very limited datasets, it is possible that porosity changes in natural fractures can be higher than the values reported thus far. Therefore, we model porosity changes by up to 15% to observe the impact on flow. The initial models decouple the correlation between increased porosity and permeability and such that there is no permeability change in the natural fracture relative to the matrix. When we use the term porosity, we mean connected porosity.  Permeability effects: The next REV model shows the impact of change in natural fracture permeability on the FP and TOFC after simulation for 30 years. For this model, the porosity in the natural fractures are kept constant with the reservoir to allow for a detailed investigation of the flow effects due to only the permeability change in the fractures. Using CAM, we model a higher permeability in the natural fractures by assigning (scaling with) a particular fracture strength (Equation (10)). As discussed in Section 4.1, an increase in strength can be related back to the natural fracture permeability using the permeability ratio Rk. The REV model ( Figure 5) uses a far field flow of 2.5 ft/year (after being scaled by the reservoir space porosity of 5%). The strengths for NF 1 to NF 5 are increased, respectively, from 0.1 ft 4 /year to 40, 160, 500, and 1000 ft 4 /year.
The results in Figure 5 show that keeping porosity constant in the natural fractures while increasing the natural fracture strength (and thus NF permeability) leads to a change in both the FP and TOF. This is in line with what is expected from the first fundamental rule for FP and TOFC [28]:  Figure 4 shows the effects of increasing NF porosity on the FP and TOFC in the reservoir space. The reservoir porosity is kept at 5%, while NF porosity changes incrementally from being equal to the reservoir space to a high of 15%. The results clearly show that the change in porosity within the natural fracture has no effect on the streamline flow paths but does affect the time-of-flight contours. In NF 1, the porosity is the same as the reservoir and as such, there is no slowdown in the TOFC. From NF 2 to NF 5, we progressively increased the porosity to 6%, 8%, 10% and finally, 15%. The model shows that for each successive porosity increase in the natural fractures, the FP stays constant but the Energies 2019, 12, 3852 13 of 36 TOF increases. As we are using a constant run time for all models, the increase in TOF results in flow not reaching as far into the reservoir space for the natural fractured with higher porosity. With no porosity change, flow reaches out to approximately 75 ft in the reservoir space. With a porosity change from 5% to 15% in the natural fractures, flow is retarded and reaches only approximately 44 ft out into the reservoir space. Thus, a 10% increase in natural fracture porosity results in a 40% reduction in lateral flow extent. This result can have great implications for accurately determining the DRV in the subsurface when the reservoir rock has a high density of natural fractures with variable porosity.
Permeability effects: The next REV model shows the impact of change in natural fracture permeability on the FP and TOFC after simulation for 30 years. For this model, the porosity in the natural fractures are kept constant with the reservoir to allow for a detailed investigation of the flow effects due to only the permeability change in the fractures. Using CAM, we model a higher permeability in the natural fractures by assigning (scaling with) a particular fracture strength (Equation (10)). As discussed in Section 4.1, an increase in strength can be related back to the natural fracture permeability using the permeability ratio R k . The REV model ( Figure 5) uses a far field flow of 2.5 ft/year (after being scaled by the reservoir space porosity of 5%). The strengths for NF 1 to NF 5 are increased, respectively, from 0.1 ft 4 /year to 40, 160, 500, and 1000 ft 4 /year. The results in Figure 5 show that keeping porosity constant in the natural fractures while increasing the natural fracture strength (and thus NF permeability) leads to a change in both the FP and TOF. This is in line with what is expected from the first fundamental rule for FP and TOFC [28]: permeability changes affect the FP and thus, the path of the streamlines is altered. Fluid is seen funneled into the higher permeability natural fractures while the TOF correspondingly decreases. Using the constant run-time of 30 years, this decrease in TOF results in fluid flow reaching further out into the reservoir space. As more of the fluid flow is funneled into the NF due to increasing strength, less of the fluid is transported in the inter-fracture domain (space between the natural fractures). In the space between NF 2 and 3 (though the FP are altered due to the increased NF permeability), fluid still flows in the inter-fracture space as shown by the streamlines. However, in the space between NF 4 and 5 (which are assigned much greater strengths) almost all the fluid flow is funneled into the natural fractures, with most of the inter-fracture space receiving no fluid. Permeability effects: The next REV model shows the impact of change in natural fracture permeability on the FP and TOFC after simulation for 30 years. For this model, the porosity in the natural fractures are kept constant with the reservoir to allow for a detailed investigation of the flow effects due to only the permeability change in the fractures. Using CAM, we model a higher permeability in the natural fractures by assigning (scaling with) a particular fracture strength (Equation (10)). As discussed in Section 4.1, an increase in strength can be related back to the natural fracture permeability using the permeability ratio Rk. The REV model ( Figure 5) uses a far field flow of 2.5 ft/year (after being scaled by the reservoir space porosity of 5%). The strengths for NF 1 to NF 5 are increased, respectively, from 0.1 ft 4 /year to 40, 160, 500, and 1000 ft 4 /year.
The results in Figure 5 show that keeping porosity constant in the natural fractures while increasing the natural fracture strength (and thus NF permeability) leads to a change in both the FP and TOF. This is in line with what is expected from the first fundamental rule for FP and TOFC [28]: permeability changes affect the FP and thus, the path of the streamlines is altered. Fluid is seen funneled into the higher permeability natural fractures while the TOF correspondingly decreases. Using the constant run-time of 30 years, this decrease in TOF results in fluid flow reaching further out into the reservoir space. As more of the fluid flow is funneled into the NF due to increasing strength, less of the fluid is transported in the inter-fracture domain (space between the natural fractures). In the space between NF 2 and 3 (though the FP are altered due to the increased NF permeability), fluid still flows in the inter-fracture space as shown by the streamlines. However, in the space between NF 4 and 5 (which are assigned much greater strengths) almost all the fluid flow is funneled into the natural fractures, with most of the inter-fracture space receiving no fluid.  The relationships between the natural fractures input parameters used in Figure 5 and the approximate equivalent natural fracture permeability (based on Equation (11)) are given in Table 1. Fracture input properties used in all subsequent flow models with enhanced natural fracture permeability in this study are included in Table 1. A matrix permeability of 100 nD is assumed. R k is calculated from Equation (11) with natural fracture permeability then back-calculated from Equation (8)  Once again, we artificially separate the effects of porosity and permeability to investigate each parameter individually. Figure 6a shows the result for completely open natural fractures set within a reservoir space of 5% porosity. The FP is unchanged but the TOF in the fractures increases dramatically. The fluid drawn from the open fracture does not require long travel paths (due to 100% fluid fill) and drawing the same amount of fluid from the inter-fracture matrix regions requires much longer travel paths in those regions outside the NF. Figure 6b shows the effect of natural fractures with very high permeability as compared to the reservoir space. The natural fractures in this model have a strength of 1000 ft 4 /year (while porosity is kept the same as that of the reservoir matrix) and fluid flow is simulated for a run-time of 30 years. The marked effect of the change in permeability is seen in the alteration of the FP as well as the decrease in TOF. With such a high fracture strength (high R k ) almost all flow is funneled through the natural fractures, with no fluid being transported via the inter-fracture domain. A similar conclusion was reported in earlier work by our research group [23].
The relationships between the natural fractures input parameters used in Figure 5 and the approximate equivalent natural fracture permeability (based on Equation (11)) are given in Table 1. Fracture input properties used in all subsequent flow models with enhanced natural fracture permeability in this study are included in Table 1.  (11) with natural fracture permeability then back-calculated from Equation (8) using the assumed matrix permeability.
Open fracture: A final scenario investigated with the REV model was the effect of a natural fracture with 100% porosity. Theoretically, this can be thought of as an open fracture in the subsurface. Once again, we artificially separate the effects of porosity and permeability to investigate each parameter individually. Figure 6a shows the result for completely open natural fractures set within a reservoir space of 5% porosity. The FP is unchanged but the TOF in the fractures increases dramatically. The fluid drawn from the open fracture does not require long travel paths (due to 100% fluid fill) and drawing the same amount of fluid from the inter-fracture matrix regions requires much longer travel paths in those regions outside the NF. Figure 6b shows the effect of natural fractures with very high permeability as compared to the reservoir space. The natural fractures in this model have a strength of 1000 ft 4 /year (while porosity is kept the same as that of the reservoir matrix) and fluid flow is simulated for a run-time of 30 years. The marked effect of the change in permeability is seen in the alteration of the FP as well as the decrease in TOF. With such a high fracture strength (high Rk) almost all flow is funneled through the natural fractures, with no fluid being transported via the inter-fracture domain. A similar conclusion was reported in earlier work by our research group [23]. All previous REV models have considered the varying effects of porosity and permeability independently of each other. Figure 7 investigates the effect of simultaneous changes of natural fracture porosity on flow, while the permeability contrast with the matrix exists (Rk > 1). In this model, we systematically change the porosity within the NF from initially being equal to that of the reservoir space of 5% ( Figure 7a) to a high of 30% (Figure 7d), all the while keeping a constant enhanced permeability in the natural fractures.
The results from the models in Figure 7 show the competing effects between porosity and permeability as defined in the fundamental rules for FP and TOF by Zuo and Weijermars [28]. Figure  7a shows the alteration in FP and decrease in TOF (fast travel times via the fractures) due to the enhanced natural fracture permeability. The successive models (Figure 7b-d) with gradually increasing porosity in the natural fracture conversely increase the TOF and thus reduce the lateral distance reached by the fluid flow in the given run-time. Although the porosity change negates the effect of the enhanced permeability in terms of lateral distance reached, the alteration of the FP by the permeability still occurs. This proves that permeability is responsible for the particle paths while both the permeability and porosity inversely affect the TOF (as stated in Zuo and Weijermars, [28]). All previous REV models have considered the varying effects of porosity and permeability independently of each other. Figure 7 investigates the effect of simultaneous changes of natural fracture porosity on flow, while the permeability contrast with the matrix exists (R k > 1). In this model, we systematically change the porosity within the NF from initially being equal to that of the reservoir space of 5% ( Figure 7a) to a high of 30% (Figure 7d), all the while keeping a constant enhanced permeability in the natural fractures. All previous REV models have considered the varying effects of porosity and permeability independently of each other. Figure 7 investigates the effect of simultaneous changes of natural fracture porosity on flow, while the permeability contrast with the matrix exists (Rk > 1). In this model, we systematically change the porosity within the NF from initially being equal to that of the reservoir space of 5% (Figure 7a) to a high of 30% (Figure 7d), all the while keeping a constant enhanced permeability in the natural fractures.
The results from the models in Figure 7 show the competing effects between porosity and permeability as defined in the fundamental rules for FP and TOF by Zuo and Weijermars [28]. Figure  7a shows the alteration in FP and decrease in TOF (fast travel times via the fractures) due to the enhanced natural fracture permeability. The successive models (Figure 7b-d) with gradually increasing porosity in the natural fracture conversely increase the TOF and thus reduce the lateral distance reached by the fluid flow in the given run-time. Although the porosity change negates the effect of the enhanced permeability in terms of lateral distance reached, the alteration of the FP by the permeability still occurs. This proves that permeability is responsible for the particle paths while both the permeability and porosity inversely affect the TOF (as stated in Zuo and Weijermars, [28]). The results from the models in Figure 7 show the competing effects between porosity and permeability as defined in the fundamental rules for FP and TOF by Zuo and Weijermars [28]. Figure 7a shows the alteration in FP and decrease in TOF (fast travel times via the fractures) due to the enhanced natural fracture permeability. The successive models (Figure 7b-d) with gradually increasing porosity in the natural fracture conversely increase the TOF and thus reduce the lateral distance reached by the fluid flow in the given run-time. Although the porosity change negates the effect of the enhanced permeability in terms of lateral distance reached, the alteration of the FP by the permeability still occurs. This proves that permeability is responsible for the particle paths while both the permeability and porosity inversely affect the TOF (as stated in Zuo and Weijermars, [28]).

Synthetic Hydraulic Fracture Models
Using the CAM model, hydraulic fractures can be modeled as either line sinks or as line sources, which is used in this study by applying the principle of flow reversal. Line sinks can show fluid withdrawal contours being forward modeled by line sources (a simple sign reversal in our equations). The effects of fluid flow of enhanced permeability and porosity in natural fractures of an otherwise homogenous reservoir space was modeled in the previous section using a constant far field flow. Models are now presented to demonstrate how natural fracture will alter fluid flow around a single hydraulic fracture. Time-dependent production data from a well completed in the Wolfcamp Formation is used in these models and prorated for fluid allocation produced by a single hydraulic fracture stage ( Figure A1). The relatively wide zones (10 ft) of altered permeability and porosity used in these models represent the effect of upscaling numerous smaller individual natural fractures (a detailed upscaling procedure is given in Appendix C). The effect of such altered zones can be clearly demonstrated visually. Each naturally fractured zone has dimensions of 10 ft width by 20 ft in length and the zones are angled at values of 45 • and 135 • from the hydraulic fracture.
The first model looks at the effect of a synthetic, single hydraulic fracture surrounded by six natural fracture zones having a higher porosity than the reservoir matrix ( Figure 8). For this model, the natural fracture zones do not attribute any additional permeability change, only porosity enhancement. Figure 8a,b has a progressively increasing porosity in the natural fracture from left to right, starting with a NF porosity of 10% in 8a and increases to 15% and 20% in Figure 8b,c. The models show that as porosity increases in the natural fracture zone, there is a decrease in the distance drained. In other words, as porosity increases, the time-of-flight also increases. The major observation from these models is that the presence of naturally fractured zones with increased porosity (and assumed fluid storage in those fractures) will decrease the distance drained away from the hydraulic fracture. The next property investigated is the effect of increased permeability (by changing the strength of the natural fractures as compared to rest of the reservoir matrix) (Figure 9). The porosity in the NF zones is kept the same as for the reservoir matrix. Therefore, we can focus solely on the permeability effect. From left to right, the strength in the natural fracture zones is progressively increased from 1000 ft 4 /day in Figure 9a to 5000 ft 4 /day and 10,000 ft 4 /day respectively in Figure 9b,c. Once again, as demonstrated in the REV models of Section 5.1, the streamlines converge into the high permeability zones and lead to larger drainage regions in the direction of the higher permeability zones. One additional point of note is that the direction of the angle of these zones in conjunction with the streamline direction influences how much effect there is on the drainage. If the naturally fractured zones are angled in the same direction as the streamlines, the effect is more pronounced than if they occur at a larger angle to the principal flow direction induced by the hydraulic fracture.  The previous models investigated the effect of altered porosity and permeability in naturally fractured zones around a hydraulic fracture independently. Figure 10 looks at the competing effects of altered porosity and permeability together. Figure 10a shows the case with an enhanced permeability in the natural fractured zones, while the porosity is kept the same as the reservoir porosity of 5%. The results show the convergence of the streamlines into these zones, resulting in a lateral extension of the DRV beyond. As we progress from left to right, Figure 10a-c shows the effect of increasing porosity in the NF zones while also having an enhanced permeability. Figure 10b has the same enhanced permeability as in Figure 10a The previous models investigated the effect of altered porosity and permeability in naturally fractured zones around a hydraulic fracture independently. Figure 10 looks at the competing effects of altered porosity and permeability together. Figure 10a shows the case with an enhanced permeability in the natural fractured zones, while the porosity is kept the same as the reservoir porosity of 5%. The results show the convergence of the streamlines into these zones, resulting in a lateral extension of the DRV beyond. As we progress from left to right, Figure 10a-c shows the effect of increasing porosity in the NF zones while also having an enhanced permeability. Figure 10b has the same enhanced permeability as in Figure 10a, but now the porosity in the natural fractured zones is increased from 5% (same as the reservoir matrix) to 10%. This model shows that although the streamlines once again converge into the zones of higher permeability, the lateral extent of the DRV is now slightly reduced due to the increased porosity. The enhanced DRV from Figure 10a has now been reduced in Figure 10b to an extent smaller (due to the porosity effect) than if there were no natural fractures. If the natural fracture porosity is increased further to 20%, the extent of the drained area shrinks much further (Figure 10c). The large changes in lateral extent and the spatial location of the DRV due to natural fractures may have significant implications for fracture and well spacing for optimum drainage. The limiting factor for improving models is the lack of fracture diagnostics for field cases (in particular the fracture permeability and porosity values). In the next section, detailed field data abstracted from the Hydraulic Fracture Test Site will be used to constrain fluid withdrawal patterns near the hydraulic fractures that drain the Wolfcamp reservoir space.
is now slightly reduced due to the increased porosity. The enhanced DRV from Figure 10a has now been reduced in Figure 10b to an extent smaller (due to the porosity effect) than if there were no natural fractures. If the natural fracture porosity is increased further to 20%, the extent of the drained area shrinks much further (Figure 10c). The large changes in lateral extent and the spatial location of the DRV due to natural fractures may have significant implications for fracture and well spacing for optimum drainage. The limiting factor for improving models is the lack of fracture diagnostics for field cases (in particular the fracture permeability and porosity values). In the next section, detailed field data abstracted from the Hydraulic Fracture Test Site will be used to constrain fluid withdrawal patterns near the hydraulic fractures that drain the Wolfcamp reservoir space.

Field Models Using Data from the Hydraulic Fracture Test Site (HFTS)
Data from the Hydraulic Fracture Test Site (HFTS; Midland Basin, West Texas) were used because the natural fracture network present in the subsurface has been characterized in prior studies for this real field case [12,56]. Six cores obtained from the Wolfcamp Formation within the Stimulated Rock Volume (SRV) near to a hydraulically fractured well were studied in detail [12]. One of the aims of the core description was to understand the primary origins of fractures in terms of hydraulic, natural and reactivated natural fractures. The density of the individual types of fractures along the core depths and the dominant orientations of the fractures obtained by Shrivastava et al. [12] were used in our study for a field-based simulation of the impact of natural fractures on the DRV development.
For the present study, certain mean values for natural fracture lengths and aperture were assumed in our models because natural fracture length and aperture values from the HFTS core samples were poorly constrained [12]. In their approximation, the latter authors used a power-law relation to generate a range for natural fracture lengths and the fracture apertures were estimated using a geomechanical fracture propagation simulator. In the present study, we constrain the fracture length to 30 ft (Table 1), corresponding to the maximum value used by Shrivastava et al. [12].

Field Models Using Data from the Hydraulic Fracture Test Site (HFTS)
Data from the Hydraulic Fracture Test Site (HFTS; Midland Basin, West Texas) were used because the natural fracture network present in the subsurface has been characterized in prior studies for this real field case [12,56]. Six cores obtained from the Wolfcamp Formation within the Stimulated Rock Volume (SRV) near to a hydraulically fractured well were studied in detail [12]. One of the aims of the core description was to understand the primary origins of fractures in terms of hydraulic, natural and reactivated natural fractures. The density of the individual types of fractures along the core depths and the dominant orientations of the fractures obtained by Shrivastava et al. [12] were used in our study for a field-based simulation of the impact of natural fractures on the DRV development.
For the present study, certain mean values for natural fracture lengths and aperture were assumed in our models because natural fracture length and aperture values from the HFTS core samples were poorly constrained [12]. In their approximation, the latter authors used a power-law relation to generate a range for natural fracture lengths and the fracture apertures were estimated using a geomechanical fracture propagation simulator. In the present study, we constrain the fracture length to 30 ft (Table 1), corresponding to the maximum value used by Shrivastava et al. [12]. Additionally, the DRV model requires inputs, for every natural fracture, of permeability and porosity. However, almost no data is present in the literature for relating in situ natural fracture porosity with permeability in the subsurface, which is why a Carman-Kozeny (CK) relation was used in our study (Appendix B).
An example of the impact of the Carman-Kozeny porosity-permeability correlation in the natural fractures, but for a still unscaled model, is given in Figure 11. The effect of the enhanced permeability in the natural fractures (Figure 11b) as compared to a single hydraulic fracture without any natural fractures nearby (Figure 11a) is to channel fluid flow faster through these high-speed zones. The effect of the enhanced permeability for this synthetic case completely outweighs any impact of the increased porosity in the natural fracture, which actually increases the time of flight (TOF) and leads to narrowly spaced TOF contours. porosity with permeability in the subsurface, which is why a Carman-Kozeny (CK) relation was used in our study (Appendix B).
An example of the impact of the Carman-Kozeny porosity-permeability correlation in the natural fractures, but for a still unscaled model, is given in Figure 11. The effect of the enhanced permeability in the natural fractures (Figure 11b) as compared to a single hydraulic fracture without any natural fractures nearby (Figure 11a) is to channel fluid flow faster through these high-speed zones. The effect of the enhanced permeability for this synthetic case completely outweighs any impact of the increased porosity in the natural fracture, which actually increases the time of flight (TOF) and leads to narrowly spaced TOF contours. The analysis of the HFTS natural fracture field data suggests that a dense network of natural fractures occurs around the hydraulic fractures [12]. The natural fracture density model based on HFTS field data generated by a discrete fracture network contained over 40,500 individual natural fractures distributed over a domain of 300 m by 300 m [56]. For tractable run times with our smaller model, the number of natural fractures can be reduced by upscaling. A similar approach was used by Kumar et al. [56], whereby the permeability tensor for the entire stimulated rock volume was determined from flowback for input in a discrete fracture network model.
The upscaling method used in the present study sought to reduce the overall number of fractures to be modeled by upscaling the natural fracture widths and fracture permeabilities (strengths) for a dense natural fracture network. Original natural fracture apertures in the subsurface were assumed to be 5 mm (0.2 inches), which follows from core observations that kinematic apertures were estimated to have been more than 1 mm wide [57]. A combination of object-based and flow-based upscaling was developed for this study, with an in-depth discussion of this topic given in Appendix C. The proposed upscaling method was applied to produce field models for DRV around a single hydraulic fracture with a representative, upscaled natural fracture distribution of the HFTS. Using the data input ranges ( Table 2) for natural fractures in conjunction with the Carman-Kozeny correlation, the final model was simulated to determine the real-life impact of natural fractures on the DRV.  The analysis of the HFTS natural fracture field data suggests that a dense network of natural fractures occurs around the hydraulic fractures [12]. The natural fracture density model based on HFTS field data generated by a discrete fracture network contained over 40,500 individual natural fractures distributed over a domain of 300 m by 300 m [56]. For tractable run times with our smaller model, the number of natural fractures can be reduced by upscaling. A similar approach was used by Kumar et al. [56], whereby the permeability tensor for the entire stimulated rock volume was determined from flowback for input in a discrete fracture network model.
The upscaling method used in the present study sought to reduce the overall number of fractures to be modeled by upscaling the natural fracture widths and fracture permeabilities (strengths) for a dense natural fracture network. Original natural fracture apertures in the subsurface were assumed to be 5 mm (0.2 inches), which follows from core observations that kinematic apertures were estimated to have been more than 1 mm wide [57]. A combination of object-based and flow-based upscaling was developed for this study, with an in-depth discussion of this topic given in Appendix C. The proposed upscaling method was applied to produce field models for DRV around a single hydraulic fracture with a representative, upscaled natural fracture distribution of the HFTS. Using the data input ranges ( Table 2) for natural fractures in conjunction with the Carman-Kozeny correlation, the final model was simulated to determine the real-life impact of natural fractures on the DRV. From the upscaling of the original natural fracture density, the outcome is a model with 12 natural fractures around the single hydraulic fracture. These 12 natural fractures are stochastically placed around the hydraulic fracture using the relevant field data all other parameters needed. Once again, the CK correlation was used to relate natural fracture permeability to porosity. Simulation of this model in CAM gives the representative DRV when affected by natural fractures (Figure 12a). Figure 12a shows that fluid is preferentially channeled through the natural fractures for the HFTS field case models. The DRV in the upscaled HFTS model is highly convolute (Figure 12a) with numerous undrained matrix zones occurring between the upscaled natural fractures created from field data. Any storativity effects of the enhanced porosity in the natural fractures remain obscured by the enhanced flow due to the enhanced permeability of the natural fractures. For comparison, the pressure plot after 1 month of production was generated using CAM (Figure 12b). Pressure was calculated in CAM by extracting the potential function from the complex potential and normalizing by the ratio of reservoir permeability and fluid viscosity [58]. For the plot presented, the pressure scale was normalized by the maximum pressure present in the reservoir at 1 month production. The lowest pressures occur near the hydraulic fractures. We utilized the process of flow reversal, which means that the highest pressures occur at the hydraulic fractures (which can be simply corrected by flipping the scale in Figure 12b). Anomalous high pressures at the tips of the natural fractures are due to singularities and associated branch cut effects occurring when high permeability contrasts (R k ) are used [59]. The progressive distortion of the pressure field near a hydraulic fracture due to the presence of natural fractures is further discussed in Section 6.2 (see also Figure 15).
The overall pressure field is greatly altered by the presence of natural fractures due to their impact on the flow pattern. The results presented here confirm that the calculated DRV do not conform 1:1 to the pressure field, making the use of pressure plots very poor proxies for reservoir drained areas.  [12]. d Values obtained from upscaling (Appendix C).
From the upscaling of the original natural fracture density, the outcome is a model with 12 natural fractures around the single hydraulic fracture. These 12 natural fractures are stochastically placed around the hydraulic fracture using the relevant field data all other parameters needed. Once again, the CK correlation was used to relate natural fracture permeability to porosity. Simulation of this model in CAM gives the representative DRV when affected by natural fractures (Figure 12a). Figure 12a shows that fluid is preferentially channeled through the natural fractures for the HFTS field case models. The DRV in the upscaled HFTS model is highly convolute (Figure 12a) with numerous undrained matrix zones occurring between the upscaled natural fractures created from field data. Any storativity effects of the enhanced porosity in the natural fractures remain obscured by the enhanced flow due to the enhanced permeability of the natural fractures. For comparison, the pressure plot after 1 month of production was generated using CAM (Figure 12b). Pressure was calculated in CAM by extracting the potential function from the complex potential and normalizing by the ratio of reservoir permeability and fluid viscosity [58]. For the plot presented, the pressure scale was normalized by the maximum pressure present in the reservoir at 1 month production. The lowest pressures occur near the hydraulic fractures. We utilized the process of flow reversal, which means that the highest pressures occur at the hydraulic fractures (which can be simply corrected by flipping the scale in Figure 12b). Anomalous high pressures at the tips of the natural fractures are due to singularities and associated branch cut effects occurring when high permeability contrasts (Rk) are used [59]. The progressive distortion of the pressure field near a hydraulic fracture due to the presence of natural fractures is further discussed in Section 6.2 (see also Figure 15).
The overall pressure field is greatly altered by the presence of natural fractures due to their impact on the flow pattern. The results presented here confirm that the calculated DRV do not conform 1:1 to the pressure field, making the use of pressure plots very poor proxies for reservoir drained areas.

HFTS Full Well Model and Implications
The previous section analyzed the impact that natural fractures modeled from field data have on the DRV around an individual hydraulic fracture. This concept is now expanded upon to determine the impact of natural fractures on DRV across multiple fracture stages representative of an entire hydraulically fractured well. The Wolfcamp production used effectively in these models had 22 stages with each stage spanning 300 ft with a total of 131 individual fracture clusters along the entire lateral. Our modeled DRV around a single hydraulic fracture is assumed to be representative of the collated drainage for all the fracture clusters per stage. Each stage has six fracture initiation points (clusters) with 50 ft spacing. The results thus show the total drainage of these six clusters when upscaled to one single hydraulic fracture.
The first model investigates the drainage based on the given 50 ft cluster spacing (corresponding to the stage spacing of 300 ft) with the assumption of a homogenous reservoir with no natural fractures (Figure 13a). Based on this stage spacing and from the DRV calculated, the multi-stage plot shows large undrained areas in between the existing DRV's after 30 years forecasted production. Results indicate that a maximum distance of 50 ft is drained perpendicularly away from the hydraulic fractures, which represents the drainage of all six fracture clusters. The plots (Figure 13a,b) show this stage spacing was sub-optimal due to the large undrained areas that can be targeted for refracs. For comparison, we model the same number of stages but now including the impact of reservoir heterogeneity using the HFTS field data on natural fractures (Figure 13b). When compared to the case with no natural fractures, the maximum area drained perpendicular to the hydraulic fracture increases from 50 ft to approximately 80 ft. Figure 13b shows that even though there is a shift in the spatial location of the DRV due to the natural fractures, this increase in lateral drainage is not enough to efficiently drain in between the fractures at this stage spacing. Assuming a modified initial fracture cluster spacing of 25 ft, down from 50ft (which corresponds to a stage spacing of 150 ft instead of the field value of 300 ft), the DRV's were modeled using CAM to investigate cases of a homogenous reservoir and heterogeneous reservoir with natural fractures (Figure 14a,b). The first case for a homogenous reservoir (Figure 14a) suggests that the reduction of the cluster spacing based on the upscaled DRV for a single stage, allows for more efficient drainage along the length of the lateral. This decrease in spacing to a more optimal value would lead to enhanced well productivity. Our method visualizes the exact DRV and the new spacing does not create adverse flow interference. In fact, the model shows that the spacing can be further optimized to slightly less than 150 ft per stage due to there still being undrained areas between the hydraulic fractures. The introduction of natural fracture heterogeneity reveals a different finding when the stage spacing is decreased to 150 ft. Natural fractures with enhanced permeability when properly oriented to the hydraulic fracture extend the lateral drained areas as shown in our models (Figure 14b). Although the natural fractures extend the drained areas, at the new stage spacing of 150 ft, there is now nearly an overlapping of the DRVs from each stage (shown by dashed black ellipses in Figure 14b). The proximity of these DRVs implies that reduction of the stage spacing to less than 150 ft will lead to flow interference that will reduce the overall recovery from the well. The conclusion from this is that when natural fractures are present, fracture stage treatment with a spacing of less than 150 ft will now be sub-optimal. These results show the importance of accounting for and properly modeling natural fractures, particularly in flow simulations for unconventional reservoirs.

Discussion
Proper modeling and forecasting of production from unconventional reservoirs needs to take into account important reservoir heterogeneity such as the presence and the impact of natural fractures. Numerous authors have noted the possible impact that natural fractures can have on production and well performance [2,24], but very few have sought to succinctly delineate and differentiate the ways in which this is possible. The present study puts forward three major mechanisms by which natural fractures can impact well productivity. Natural fractures present in the subsurface can affect well productivity via (1) enhanced permeability, (2) enhanced storativity, and (3) reactivation of natural fractures as extensions to the created hydraulic fracture network. By the use of a simple analytical streamline simulator based on complex analysis methods (CAM), we visualized the drainage patterns around hydraulic fractures by Eulerian particle tracking. The effects of natural fractures, in particular the enhanced permeability and storativity were investigated systematically and the results show that the drainage patterns (DRV) can be greatly altered by the presence of these reservoir heterogeneities.

Storativity Impact of Natural Fractures
Natural fractures present in the subsurface show a range of measured porosity from 2% to 7% [25] but these measured data sets are very limited in sample size and it is believed that porosity ranges may include even higher values. The altered mineralogy in these natural fractures can lead to a porosity and permeability that is vastly different to that of the unfractured reservoir matrix. With regards to natural fractures present in the Permian Basin, Forand et al. [24] stated that "despite natural fractures having a calcite fill, the permeability contrast between the fracture and matrix is likely high enough that the healed fractures may be preferential hydrocarbon pathways. Combining this dominant character with the orientation of natural fractures to maximum horizontal principal stress has the potential to affect the efficiency of hydraulic fractures and the size of the total connected and stimulated rock volume." The change in permeability will also result in an increased porosity, which we see as a cause of enhanced storativity for reservoir fluids.
Enhanced storativity can contribute to better well performance as these naturally fractured regions will have a larger hydrocarbon fluid supply that may last longer [23]. The impact of enhanced storativity in natural fractures on the drainage area around a well is for the first time visualized in our results. Starting with a simple REV model (Figure 4), the effect of increased porosity is seen to slow the time-of-flight (TOF) in the natural fracture as compared to the matrix. Once again, this proves that porosity changes do not affect streamline patterns but only the time-of-flight [28]. When applied to naturally fractured zones around a hydraulic fracture (Figure 8), the increase in the TOF results in a slower expansion of the DRV in the natural fracture zones compared to the rest of the matrix with a lower porosity. This leads to a decrease in the lateral distance drained away from the hydraulic fracture and can thus impact the optimum fracture cluster spacing distance. For a highly naturally fractured reservoir with higher storativity, the well spacing could be decreased compared to a reservoir with no natural fractures, as the drained area laterally would be smaller. This ability to increase the number of wells without introducing interference effects (by draining the same area with multiple hydraulic fractures) will lead to higher recoveries per acreage.

Enhanced Permeability vs. Enhanced Storativity
For natural fractures with higher permeability, fluid moves preferentially through these high-velocity conduits. REV models for natural fractures with various permeabilities (Figure 5), modeled by individually specified natural fracture strengths in our CAM simulation, show that as fluid moves via the natural fractures, some of the matrix areas between the natural fractures are bypassed or left undrained. When applied to flow around a single hydraulic fracture ( Figure 9) the preference for flow through the higher permeability zones creates enhanced lateral drainage in the areas where the drainage plumes near the tips of the natural fractures reach deeper into the lateral reservoir space. Our results show that altered permeability impacts both the streamline patterns (convergence into natural fractures) and TOF. For a greater permeability, the TOF reduces in the natural fractures as compared to the TOF in the matrix. Thus, natural fractures with enhanced permeability can lead to greater lateral drainage with the caveat that there is the possibility of bypassed areas between the natural fractures that can still contain hydrocarbons.
The synthetic models all assumed variations in the porosity being possible independent of permeability changes. In reality, this is not the case as an increase in the effective porosity commonly correlates to an increase in permeability. Nonetheless, the synthetic examples clearly highlight that increased porosity leads to an increase of the TOF (i.e., flow is slowed down in the higher porosity region), whereas increased permeability reduces the TOF (i.e., flow if quickened). The latter also alters the flow paths in the reservoir. This leads to a competing effect of higher porosity reducing the lateral DRV, with greater permeability increasing the lateral DRV assuming otherwise similar production (as used in our models).
The key questions now become: (1) which parameter (permeability vs. porosity) has the more dominant impact on the drainage pattern? and (2) how can one correlate any increases in porosity with permeability, and vice versa? Data for natural fracture porosity values are very limited and any natural fracture permeability values are for typically reactivated fractures that connect directly to the hydraulic fracture. Due to this paucity of data, this paper made use of the commonly used Carman-Kozeny (CK) correlation for determining permeability based on a given natural fracture porosity. The results show ( Figure 11) that using this correlation with a limited number of natural fractures, the permeability effect far outweighs the storativity of the enhanced porosity.
The HFTS case (Figure 12), using field data for natural fracture representation (based on natural fracture upscaling), shows that once the CK correlation is used, the impact of the natural fracture enhanced permeability (lateral extension of DRV and undrained matrix between natural fractures), vastly outweighs the storativity effect of said natural fractures. The DRV and pressure field distortion for the HFTS (Figure 12a,b) provide a specific example of what is a generic effect. For example, Figure 15a-d show the pressure field around a single hydraulic fracture without any natural fractures present (Figure 15a) and the stepwise distortion of the associated pressure field due to the presence of one, two and six natural fractures (Figure 15b-d). It should be noted that our models have the highest pressures at the hydraulic fracture due to the flow reversal modeling used (whereby fluid is placed back into the reservoir via the hydraulic fractures at the same rate as that produced).
dominant impact on the drainage pattern? and (2) how can one correlate any increases in porosity with permeability, and vice versa? Data for natural fracture porosity values are very limited and any natural fracture permeability values are for typically reactivated fractures that connect directly to the hydraulic fracture. Due to this paucity of data, this paper made use of the commonly used Carman-Kozeny (CK) correlation for determining permeability based on a given natural fracture porosity. The results show (Figure 11) that using this correlation with a limited number of natural fractures, the permeability effect far outweighs the storativity of the enhanced porosity.
The HFTS case (Figure 12), using field data for natural fracture representation (based on natural fracture upscaling), shows that once the CK correlation is used, the impact of the natural fracture enhanced permeability (lateral extension of DRV and undrained matrix between natural fractures), vastly outweighs the storativity effect of said natural fractures. The DRV and pressure field distortion for the HFTS (Figure 12a,b) provide a specific example of what is a generic effect. For example, Figure  15a-d show the pressure field around a single hydraulic fracture without any natural fractures present (Figure 15a) and the stepwise distortion of the associated pressure field due to the presence of one, two and six natural fractures (Figure 15b-d). It should be noted that our models have the highest pressures at the hydraulic fracture due to the flow reversal modeling used (whereby fluid is placed back into the reservoir via the hydraulic fractures at the same rate as that produced).

Model Strengths and Limitations
The CAM models presented here are grid-less and meshless, unlike the more often used numerical methods in industry. Due to their being grid-less, CAM is much less computationally intensive than finite-volume/difference numerical methods with the added advantage of high resolution at the scale of the hydraulic and much smaller natural fractures. Other strengths of the CAM model to accurately determine the impact of natural fractures on drained rock volumes comes in the form of this analytical method having closed form solutions as well transparency in all steps of the methodology [23]. The present study is limited to flow in 2D as well as only modeling single phase fluid flow. As the natural fractures are modeled as individual discrete elements, the model would become cumbersome to use and computationally expensive if large scale, stochastically generated natural fracture networks are taken as inputs. This is the rationale behind the use of upscaling methods to represent natural fractures used in the field scale models. In reality, the geometry of both the natural fractures (in terms of inclination angle in 3D) and the hydraulic fractures (as fractal networks instead of simple bi-planar features) are much more complex than that represented here. In spite of these simplifying and reductionist model assumptions (as all other models also have), the CAM tool developed in this paper to include the impact of natural fractures can be used as a quick and simple method to screen optimum hydraulic fracture spacing and to support and direct well spacing decisions in naturally fractured reservoirs. What the 2D studies provide are very valuable systematic insight that will benefit the improvement of 3D model studies as well. Accounting for 3D dimensionality may make for more realistic models, but when coupled with flow, may also disguise some of the systematic effects visualized in our 2D models of flow in hydraulically and naturally fractured reservoirs.

Practical Implications
The impacts of natural fractures on production in unconventional wells are still debated. However, the interaction of the in-situ stress, hydraulic fractures and natural fractures could be leveraged to optimize well path planning and completions designs [24]. In this study, we distinguished three major mechanisms via which natural fractures may impact flow and, implicitly, acreage productivity. Flow models based on CAM show that enhanced natural fracture permeability and porosity can alter the DRV shapes and spatial location greatly. This can have implications for the spacing of both hydraulic fractures and wells once the nature of the natural fracture network in the subsurface has been accurately characterized. For formations with highly permeable natural fractures, well spacing should be slightly increased to avoid interference as the DRVs would otherwise overlap.
However, this assumes that the spacing is based on DRV modeling. If based on pressure interference models only, our previous work [33,60] argues that such pressure interference occurs for much larger well spacing and fracture spacing. However, such pressure interference should not be used as the sole criterion for well and fracture spacing decisions because of the over one order of magnitude time-lag between the pressure front and the tracer front propagation in ultra-low permeability reservoirs [61].
The models presented emphasize how the spatial orientation, location and lateral extent of the DRV are vastly impacted by the presence of natural fractures. Fluid flows preferentially through the highly conductive natural fractures, altering the shape of the DRV around hydraulic fractures. Any undrained matrix zones that have been bypassed due to flow channeling into the natural fractures with high flow rates can then be preferentially targeted for refracturing. For rock formations where the stress regimes preferentially allow for reactivation of natural fractures to form an extension of the hydraulic fracture, cluster spacing can be decreased to allow for the creation of the largest, most complex fracture network that gives the greatest access to the hydrocarbons trapped in the low permeability reservoir rock [21].

Conclusions
Natural fractures present in the subsurface are a major form of heterogeneity in both conventional and unconventional hydrocarbon reservoirs. Highly conductive natural fractures may provide preferential pathways for fluid withdrawal to the production wells, which is why natural fractures are highly crucial for well design decisions (especially in unconventional reservoirs). The major conclusions from our analysis on the impact of natural fractures on subsurface flow are (1) Natural fractures can affect reservoir flow through three major mechanisms: (i) by enhancing permeability, (ii) by altering the porosity in the fractures, leading to increased storativity, and (iii) by becoming extensions of the hydraulic fracture network due to reactivation. (2) Enhanced permeability in natural fractures creates high velocity flow zones which preferentially channel fluid flow through them. At high enough permeabilities (or natural fracture strengths as used in our models), this preferential pathway to flow leads to bypassed regions in the matrix blocks between the natural fractures, which are left undrained. These undrained matrix regions can then be targeted by refracturing to improve recovery factors from hydraulically fractured horizontal wells. (3) Altered porosity or enhanced storativity (due to natural fractures with a higher porosity than the reservoir matrix as investigated in synthetic models) leads to a decrease in the lateral extent of the DRV. The impact of both natural fracture storativity and permeability greatly affect the shape and extent of the DRV around the hydraulic fractures. (4) The Carman-Kozeny (CK) relation was used to determine the relative impacts of the correlated porosity and permeability in natural fractures on the DRV development. Results based on the CK correlation show that the enhanced flow due to permeability far outweighs any storativity effects (even if natural fractures were to have a higher porosity than the reservoir matrix). (5) Use of a hybrid object-based and flow-based method for upscaling allows for the modeling of a high-density natural fracture network. Upscaling is needed to reduce the number of natural fractures modeled while keeping the equivalent permeability the same. (6) Field data on in-situ natural fracture characteristics such as porosity and permeability is sparse and lacking in the literature. Industry needs to ensure collection of such data for use in reservoir models to accurately determine subsurface flow and drainage volumes. (7) Proper analysis of natural fracture data and the predominant mechanism by which it will affect flow will lead to accurate DRV calculations in the subsurface. From these determined DRV (based on a well type curve), fracture cluster spacing and well spacing could possibly be optimized.

Appendix A. Flux Modeling and Production Allocation for Hydraulic Line Sink Models
Flux allocation was proportional to the relative surface areas of each hydraulic fracture. The flux allocation algorithm used is as follows: Z is a conversion factor of 5.61 to convert from barrels to ft 3 as q well is given in bbls/day and q k in ft 3 /day. The WOR though seemingly very high (Table A1) is usual for these Wolfcamp completed wells [9]. S is the prorated factor to scale the total well production, for one hydraulic fracture stage from this well (which had a total of 22 stages): After calculation of the flux using the given algorithm, we next determine the appropriate strength based on the time-dependent flow to use in the velocity and pressure potential equations. This strength is scaled by reservoir properties such as the formation volume factor (B), porosity (n), residual oil saturation (R o ) [60] and hydraulic fracture height (H k ) and is determined as follows: Production from a well was history matched using the Arp's hyperbolic decline curve method and then forecasted for the 30 year production life of the well ( Figure A1). This total well production was then allocated back into a hydraulics fracture representative of a single stage of the well. The production well used has 22 stages with a total of 131 individual fracture clusters. As such our DRV models are representative for the upscaled production from each stage.  Production from a well was history matched using the Arp's hyperbolic decline curve method and then forecasted for the 30 year production life of the well ( Figure A1). This total well production was then allocated back into a hydraulics fracture representative of a single stage of the well. The production well used has 22 stages with a total of 131 individual fracture clusters. As such our DRV models are representative for the upscaled production from each stage. Figure A1. Arps hyperbolic decline curve model used to history match field data to give: (1) Production rate (STB/day, left scale, red curve) and (2) cumulative production (STB, right scale, green curve) for well after forecasted time of 30 years (10,958 days).

Appendix B. Carman-Kozeny Relation for Estimating Natural Fracture Permeability from Porosity
For the field models looking at use of the natural fracture data and its impact on DRV, the Carman-Kozeny correlation was used to determine an effective porosity-permeability relationship. The generic Carman-Kozeny correlation is given by [62]: Figure A1. Arps hyperbolic decline curve model used to history match field data to give: (1) Production rate (STB/day, left scale, red curve) and (2) cumulative production (STB, right scale, green curve) for well after forecasted time of 30 years (10,958 days).

Appendix B. Carman-Kozeny Relation for Estimating Natural Fracture Permeability from Porosity
For the field models looking at use of the natural fracture data and its impact on DRV, the Carman-Kozeny correlation was used to determine an effective porosity-permeability relationship. The generic Carman-Kozeny correlation is given by [62]: This well-known correlation seeks to link the permeability of a porous medium (in our case natural fractures with a predetermined porosity) to the porosity along with other rock properties. β represents the shape factor of the rock and is a constant characteristic for a particular type of granular material, S is known as the specific surface area and is the ratio of the total interstitial surface area to the bulk volume [62]. T is the hydraulic tortuosity defined as by the equation: where λ represents the mean length of fluid particle paths and the variable − L gives the straight-line distance through the medium in the direction of macroscopic flow. We adopted a T value of 1.41 [62] and a β of 3 for the pore shape coefficient for thin cracks [63]. The specific surface area by volume (S) is calculated from the specific surface area by weight and the average density using data from Wolfcamp formation samples. Specific surface areas are given by Tinni et al. [64] for various particle sizes in the Wolfcamp formation with an average specific surface area of 9.36 m 2 /g. Using this value in conjunction with the average Wolfcamp formation density of 2.73 g/cm 3 [65], S is calculated at 2.55 × 10 7 m −1 . Using these values with a given natural fracture porosity, natural fracture permeability is then calculated and converted to the equivalent strength using Equation (11) for use in the CAM models. An example of the correlation is given in Table A2 with the first row values used for Figure 11b.   Our present formulation for upscaling the permeability in fractured porous media is a hybrid between the object-based and flow-based upscaling methods. The object-based upscaling (Equations (A6)-(A12)) is first used to reduce the total number of natural fractures used in the model (essentially decreasing the natural fracture density). Next, the flow-based method (Equations (A13) and (A14)) is used with the upscaled fracture density to ensure the equivalent permeability for the REV of concern remains identical to the prototype.

Object-based upscaling step:
To demonstrate the validity of the proposed procedure we consider two similar REV's ( Figure A3). Our present formulation for upscaling the permeability in fractured porous media is a hybrid between the object-based and flow-based upscaling methods. The object-based upscaling (Equations (A6)-(A12)) is first used to reduce the total number of natural fractures used in the model (essentially decreasing the natural fracture density). Next, the flow-based method (Equations (A13) and (A14)) is used with the upscaled fracture density to ensure the equivalent permeability for the REV of concern remains identical to the prototype.

Object-based upscaling step:
To demonstrate the validity of the proposed procedure we consider two similar REV's ( Figure  A3). For REV 2: Assuming k f and w REV are constant and equating the equations for REV 1 and REV 2 we arrive at; The number of fractures in REV 1 can be determined from the natural fracture density and REV width and length (L REV ); N 1 = NF density1 × w REV × L REV Substituting for N 1 ; Based on a user defined value for a new natural fracture width we can upscale from N 1 fractures to a lower value of N 2 natural fractures which is more practical for use in discrete natural fracture models, including CAM used in our study

Validation of object-based upscaling step:
The proposed object-based upscaling (reduction) of the number of natural fractures in a given reservoir area was validated using the flow-based upscaling method. For the models with N 1 and N 2 fractures, the velocities are calculated in and outside of the natural fractures and the permeability tensor ellipses are generated. To properly account for the reduction of the number of fracture and equivalent upscaling, the assigned natural fracture permeabilities of the original prototype (υ f 1 w f 1 ) and upscaled models (υ f 2 w f 2 ) needed to maintain the same equivalent permeability are given by; where υ f 1 (t) is the original strength prior to upscaling, and υ f 2 (t) is the new strength (which are proxies to the permeability in our models) to be used after upscaling the number of natural fractures with the corresponding fracture width change. This procedure is demonstrated via the upscaling of natural fractures at an angle of 45 • to the far field flow starting with Figures A4 and A5, up to the final upscaled REV in Figure A6. The proposed object-based upscaling (reduction) of the number of natural fractures in a given reservoir area was validated using the flow-based upscaling method. For the models with N1 and N2 fractures, the velocities are calculated in and outside of the natural fractures and the permeability tensor ellipses are generated. To properly account for the reduction of the number of fracture and equivalent upscaling, the assigned natural fracture permeabilities of the original prototype ( 11 ff w  ) and upscaled models ( 22 ff w  ) needed to maintain the same equivalent permeability are given by;      The above results show that with a reduction in the number of natural fractures by object-based upscaling within a defined REV, using the appropriate upscaling for fracture width and permeability in the natural fractures, the equivalent permeability remains constant. By using this method, we can upscale a realistic fracture density to a manageable number of natural fractures for use in the CAM models for DRV calculations. This upscaling methodology was applied in the next section to field data from the Hydraulic Fracturing Test Site (Midland Basin, West Texas, with completions in the Wolfcamp Formation).

Application of object-based and flow-based upscaling to HFTS field model:
This section makes use of the proposed combination of object-based and flow-based upscaling to reduce the natural fracture density used by Shrivastava et al. [12] in their model to match the data from the HFTS. Selecting a REV located around a hydraulic fracture of 125 ft in length by 45 ft in height above the hydraulic fracture corresponds a true density of 210 natural fractures with an assumed width of 0.2 inches. The 210 fractures are reduced in the proposed upscaling procedure, making use of Equation (A17), and adopting an upscaled natural fracture width of 6 inches (based on object-based upscaling), results in 6 natural fractures of length 30 ft. These 6 natural fractures have fracture centers and angles (kept in range of HFTS data) that are stochastically generated within the specified REV both below and above the hydraulic fracture. This results in a total of 12 upscaled natural fractures that are used in the final HFTS field model ( Figure A7). The CK correlation was used with a final upscaled strength of 155 ft 4 /year, which gives a corresponding porosity of 7.32% within the natural fractures. The above results show that with a reduction in the number of natural fractures by object-based upscaling within a defined REV, using the appropriate upscaling for fracture width and permeability in the natural fractures, the equivalent permeability remains constant. By using this method, we can upscale a realistic fracture density to a manageable number of natural fractures for use in the CAM models for DRV calculations. This upscaling methodology was applied in the next section to field data from the Hydraulic Fracturing Test Site (Midland Basin, West Texas, with completions in the Wolfcamp Formation).

Application of object-based and flow-based upscaling to HFTS field model:
This section makes use of the proposed combination of object-based and flow-based upscaling to reduce the natural fracture density used by Shrivastava et al. [12] in their model to match the data from the HFTS. Selecting a REV located around a hydraulic fracture of 125 ft in length by 45 ft in height above the hydraulic fracture corresponds a true density of 210 natural fractures with an assumed width of 0.2 inches. The 210 fractures are reduced in the proposed upscaling procedure, making use of Equation (A17), and adopting an upscaled natural fracture width of 6 inches (based on object-based upscaling), results in 6 natural fractures of length 30 ft. These 6 natural fractures have fracture centers and angles (kept in range of HFTS data) that are stochastically generated within the specified REV both Energies 2019, 12, 3852 33 of 36 below and above the hydraulic fracture. This results in a total of 12 upscaled natural fractures that are used in the final HFTS field model ( Figure A7). The CK correlation was used with a final upscaled strength of 155 ft 4 /year, which gives a corresponding porosity of 7.32% within the natural fractures.