AN ANALYTICAL INVESTIGATION OF BOUNDARIES IN NATURALLY FRACTURED UNCONVENTIONAL
RESERVOIRS
by
c
Copyright by Judson T. Greenwood, 2015 All Rights Reserved
A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial fulfillment of the requirements for the degree of Master of Science (Petroleum Engineering). Golden, Colorado Date Signed: Judson T. Greenwood Signed: Dr. Erdal Ozkan Thesis Advisor Golden, Colorado Date Signed: Dr. Erdal Ozkan Interim Department Head and Professor Department of Petroleum Engineering
ABSTRACT
This research presents a heuristic approach to develop an analytical model to study the effects of a stimulated zone in a fractured unconventional reservoir and the inherent boundaries that are observed. To simulate a stimulated reservoir volume (SRV) around the fractured horizontal well surrounded by a virgin outer reservoir, three separate solutions are generated and superimposed: Solution 1 – a multiply fractured, horizontal-well in an infinite-acting, homogeneous reservoir with the properties of the outer zone; Solution 2 – a multiply fractured, horizontal-well in a bounded, homogeneous (un-fractured) reservoir with the properties of the outer reservoir; and Solution 3 – a multiply fractured, horizontal-well in a bounded, naturally fractured reservoir with the properties of the stimulated zone. The solution for the composite reservoir consisting of a stimulated (naturally fractured) reservoir surrounded by an infinite acting, un-fractured (virgin) reservoir is obtained by subtracting Solution 2 from Solution 1 and then adding Solution 3. The same approach is also applied to develop a solution for the case where there is an additional transition zone between the SRV and the outer (virgin) reservoir. This method creates an approximate solution for the composite-reservoir system. Although the model is derived analytically, computations require numerical methods and the model is therefore referred to as semi-analytical.
The model is verified against literature models and an industry numerical simulator to find its limitations. This verification shows that the accuracy of the model is dependent on the size of the stimulated zone. For a large stimulated zone, because the flux profiles along the boundaries of the fractured (Solution 3) and un-fractured (Solution 2) reservoirs are not equal, the model over-predicts the drawdown pressure. However in real-world ex-amples of multiply fractured horizontal wells, the stimulated zone is much smaller and the model closely matches the drawdown pressures calculated in the numerical simulator. There-fore, the heuristic approach used in this work leads to an ad-hoc solution for the common
configurations of fractured horizontal wells in shale reservoirs.
Several synthetic examples are considered to show that the solution developed in this work can be used to identify the flow regimes after the effect of the stimulated reservoir boundary (that is, the fracture tip effects) are felt. This is an advantage over the commonly used trilinear model when the diffusivities of the stimulated and virgin reservoirs are comparable. Although not explored in this research, the ultimate utility of the proposed approach is in modeling multiple fractured-horizontal-wells to study the interference among SRVs. The fracture enhancement and extent influence the productivity of the well more than any other parameter and should be of utmost importance to characterize. And last, this model can be used with other tools to identify optimal full field development.
TABLE OF CONTENTS
ABSTRACT . . . iii
LIST OF FIGURES . . . viii
LIST OF TABLES . . . xi ACKNOWLEDGMENTS . . . xii DEDICATION . . . xiii CHAPTER 1 INTRODUCTION . . . 1 1.1 Motivation . . . 1 1.2 Objectives . . . 2 1.3 Method of Study . . . 3
1.4 Contributions of the Study . . . 5
1.5 Organization of the Thesis . . . 6
CHAPTER 2 LITERATURE REVIEW . . . 7
2.1 Shale Properties . . . 7
2.2 Dual Porosity . . . 10
2.3 Finite-Conductivity Vertical Fracture Model . . . 11
2.4 Multiply-Fractured Horizontal Well Model . . . 13
2.5 Trilinear Model . . . 15
CHAPTER 3 MATHEMATICAL MODEL . . . 18
3.1 Details of the Model Development . . . 19
3.1.2 Development and Calculation of the Composite-Reservoir Solution . . . 21
3.2 Model Verification . . . 27
3.2.1 Uniform-Flux (Single-Segment) Fracture Solution in a Bounded Reservoir . . . 27
3.2.2 Finite-Conductivity (Multiple Segments), Single-Fracture Solution . . . 29
3.2.3 Multiply-Fractured Horizontal Well Solution . . . 33
3.2.4 Composite-Reservoir Solution; Verification with Numerical Model . . . 35
3.3 Computational Considerations . . . 40
3.3.1 Number of Fracture Segments . . . 40
3.3.2 Infinite Sums and Late Time Calculation . . . 42
3.3.3 Finite Conductivity . . . 43
CHAPTER 4 ANALYSIS OF RESULTS . . . 45
4.1 Case 1: Naturally Fractured Inner Zone with Not Naturally Fractured Outer Zone . . . 51
4.2 Case 2: Inner Zone is Densely Fractured and Outer Zone is Sparsely Fractured . . . 53
4.3 Case 3: Not Stimulated Reservoir Volume in a Naturally Fractured Reservoir . 55 4.4 Case 4: Intermediate Stimulated Zone . . . 56
4.5 Case 5: Spacing Effects . . . 59
CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS . . . 61
5.1 Conclusions . . . 61
5.2 Recommendations for Future Work . . . 63
NOMENCLATURE AND ABBREVIATIONS . . . 64
LIST OF FIGURES
Figure 1.1 Example Case 1 - An example case diagram describing the layout of the inner and outer zone and their respective boundaries. . . 4 Figure 1.2 Example Case 2 - Another example case diagram building off of
Figure 1.1 with the addition of the transition zone and transition
boundary. . . 5 Figure 2.1 Shale basins in the United States with callouts for unconventional
formations . . . 9 Figure 2.2 Shale basins worldwide . . . 9 Figure 2.3 Warren and Root pictorial representation of an actual natural fractured
reservoir and a cube modeled dual porosity reservoir (1963). . . 10 Figure 2.4 Geometry of Cinco-Ley et al. vertical fracture finite conductivity model
(1978). . . 12 Figure 2.5 Cinco-Ley et al. figure showing that dimensionless conductivities,
labeled here as (kfbf)D, above 30 display near-infinite conductivity
behavior (1978). . . 13 Figure 2.6 Raghavan et al.’s nomenclature for the four flow regimes in a
multiply-fractured horizontal well (1997). . . 14 Figure 2.7 Raghaven et al.’s expanded Cinco-Ley et al. figure showing diminishing
returns for adding additional fractures (1997). . . 15 Figure 2.8 Brown and Ozkan’s trinlinear model geometry (2009). . . 16 Figure 3.1 Single fracture segment model verification against Ozkan (1988) Table
2.6.1. . . 29 Figure 3.2 Finite conductivity, single fracture, multiple segment model verification
against Cinco-Ley et al. (1978) Table 1. . . 33 Figure 3.3 Pictorial representation for the model verification of Raghavan et al.
Figure 3.4 Finite conductivity, five multiply-fractured horizontal well, multiple
segment model verification against Raghavan et al. (1997) Table 2. . . 35 Figure 3.5 Pictorial representation of the CMG (2013) numerical model verification. . 37 Figure 3.6 Numerical model verification for 0.1 and 0.5 ηD cases in the small SRV
case. . . 38 Figure 3.7 Numerical model verification for 0.1 and 0.5 ηD cases in the large SRV
case. . . 39 Figure 3.8 Effects of the number of fracture segments for an infinite conductivity
case. . . 41 Figure 3.9 Effects of the number of fracture segments for a finite conductivity case. . 41 Figure 3.10 Effect of finite conductivity in early time showing that infinite
conductivity behavior requires very high values of dimensionless
conductivity. . . 44 Figure 4.1 Case 1 - Natrurally fractured inner zone and non-naturally fractured
outer zone. . . 46 Figure 4.2 Case 2- Reservoir system in which the inner region is densely fractured
and the outer region is sparsely fractured. . . 47 Figure 4.3 Case 3 - No stimulated reservoir volume in a naturally fractured
reservoir. . . 48 Figure 4.4 Case 4- Intermediate sparsley fractured zone outside the densely
fractured zone in a non-naturally fractured reservoir. . . 49 Figure 4.5 Case 5- Well spacing examination using the mulitply hydraulic
fractured composite model. . . 50 Figure 4.6 Behavior of the composite model showing the individual terms from
Equation 3.1 and the divergence/convergence of each of those solutions. . 52 Figure 4.7 Replication of the trilinear model using no-flow boundaries placed at
the heel and toe of the horizontal well and comparision to the
composite model without the boundary. . . 54 Figure 4.8 The well response of a densely fractured SRV in a sparsely fractured
Figure 4.9 Well responses in a sparsely fractured reservior and one without the
fractures showing the importance that natural fractures can contribute. . 56 Figure 4.10 Early time plot showing the differences between the fractures systems
in the bounded and intermediate zones. . . 57 Figure 4.11 Late time plot showing the transistion to matrix flow. . . 58 Figure 4.12 Modeling different spacing units and their associated no-flow boundaries. . 60
LIST OF TABLES
Table 3.1 Single fracture segment model verification against Ozkan (1988) Table
2.6.1. . . 28 Table 3.2 Finite conductivity, single fracture, multiple segment model verification
against Cinco-Ley et al. (1978) Table 1. . . 30 Table 3.3 Finite conductivity, single fracture, multiple segment model verification
against Cinco-Ley et al. (1978) Table 1. (continued) . . . 31 Table 3.4 Finite conductivity, single fracture, multiple segment model verification
against Cinco-Ley et al. (1978) Table 1. (continued) . . . 32 Table 3.5 Input parameters for model verification against Raghavan et al. (1997)
Table 2. . . 34 Table 3.6 Finite conductivity, five multiply-fractured horizontal well, multiple
segment model verification against Raghavan et al. (1997) Table 2. . . 36 Table 3.7 CMG (2013) numerical model verification constant input parameters. . . . 38 Table 3.8 CMG (2013) Numerical model verification variable input parameters. . . . 38 Table 3.9 Case inputs to show different behavior with number of fracture segments. . 40 Table 4.1 Fracture and matrix properties for analysis for Chapter 4 cases. . . 51 Table 4.2 General properties for analysis for Chapter 4 cases. . . 51 Table A.1 FORTRAN 2003 files. . . 68
ACKNOWLEDGMENTS
I would like to thank several people for this thesis in no particular order. Thank you Dr. Erdal Ozkan for your patience with me over too many years for my degree. I am proud to have been able to work with you and have you as an advisor. My family and friends have always supported me incredibly over the years and are responsible for all the success I have enjoyed. Jeff, I would not have been able to write the program for this thesis without your help. You spent several days talking me through FORTRAN and helping to solve problems with my ineptitude of writing code. I will never forget your tremendous help. Aunt Heidi, thank you for pushing me to finish my degree when I was dismayed. Alissa, thank you for reminding me how much I wanted my masters and for making sure I finish on time. And thanks to all the people at Noble Energy for their professional and financial support.
I dedicate this work to my late grandfather, Donald Greenwood, Sr. Without him, I would not be an engineer, would not have a curiosity for how things work, and would not be the
CHAPTER 1 INTRODUCTION
The work in the following thesis was performed for the Marathon Center of Excellence for Reservoir Studies in the Petroleum Engineering Department of the Colorado School of Mines to fulfill the requirements for the degree of Master of Science in Petroleum Engineering. This research examines the effect of boundaries in unconventional reservoir systems.
1.1 Motivation
Hydrocarbon production continues to dominate the world supply of energy. Coal has been used for many centuries; oil production became popular worldwide in the late 19th century; and nowadays, natural gas is an important addition to hydrocarbon production. In the United States, the Energy Information Administration (US EIA, 2014), shows that 82% of energy consumption comes from hydrocarbons. Of that 82%, 18% comes from coal, 37% from oil, and 27% from natural gas.
Natural gas is becoming a very important factor for electricity generation. The EIA projects that coal electricity generation will decrease while natural gas electricity generation will increase (US EIA, 2014). The most important reason for this is environmental stress and supply. Coal is considered a much more dirty fuel compared to natural gas, and therefore many coal-fired power plants are either being shut down, or converted to natural gas-fired while natural gas reserves are simultaneously increasing.
Natural gas production in the US has increased mainly due to “unconventional resources.” The Society of Petroleum Engineers (SPE) defines unconventional resources as “hydrocar-bons from unconventional and more difficult to produce resources such as shale gas, shale oil, tight gas, and tight oil, coal seam gas/coalbed methane and hydrates” (SPE, 2014). These unconventional reservoirs account for 50% of US production. Of this group, shale gas has seen the biggest growth of all natural gas by doubling production in two years and doubling
technically recoverable reserves in one year. Shale gas is forecasted to become the dominant supplier of natural gas within a few years (US EIA, 2014).
While natural gas is becoming ever so prevalent in the US, the low price environment for natural gas and high prices for oil have shifted many companies to look for more oil-rich formations. In the last three years, most of the rigs in the US have been drilling for oil prospects and gas with associated condensate liquids. The US reduced their net oil imports to 40% in 2012 and expect further decline to 25% by 2016. Most of this incremental production is due to unconventional oil formations such as the Bakken, Eagleford, Niobrara, Wolfcamp and others which will rise to over 50% of total oil production by 2021. (US EIA, 2014)
This sharp increase in production in both oil and natural gas is due to the technology of multiply-fractured horizontal wells. This style of well design and completion has been in existence for over two decades. However, the technology has greatly improved over the past several years, therefore unlocking the key to shale gas and shale oil production. While petroleum engineers have started to master the drilling and completion processes, there are considerable challenges in the area of reservoir engineering. The motivation of this thesis is to examine the effects of the composite-reservoir structure of the drainage areas of multiply-fractured horizontal wells in tight, unconventional reservoirs. This includes the effects of both apparent boundaries between the stimulated inner reservoir region and the outer region as well as boundaries in the outer reservoir system caused by well-to-well interference. This examination includes investigation of under what reservoir and stimulation properties do boundary effects become significant.
1.2 Objectives
The goal of this research is to develop a practical solution to examine the flow behavior of multiply-fractured horizontal wells in unconventional reservoirs and, in particular, the flow behavior as boundaries are observed. As explained in in Chapter 3, the approach used to
first objective of the research is to examine the validity of the model. Each part of the model is tested against its respective literature examples. The final model is tested against a simple numerical model also. The second objective of this research is to build synthetic cases, which can exemplify the production in unconventional reservoirs and help us study the boundary effects.
1.3 Method of Study
This research uses analytical methods to model production from multiply-fractured hor-izontal wells in tight, unconventional reservoirs and study the reservoir flow behavior. The derivation of the solution is analytical, however, the computation of the solutions requires numerical methods therefore the solution is referred to as semi-analytical. Furthermore, the approach to develop the solution is heuristic and ad hoc. The results are verified through comparisons with the approximate cases and numerical simulations. It must also be empha-sized that this thesis does not consider stress-dependent properties; such considerations can be incorporated in the solution fairly easily, but the resulting non-linear expressions require development of iterative computational procedures, which are beyond the objectives of this research.
To study the reservoir behavior, several cases were built and will be explained in detail in Chapter 4. Figure 1.1 shows a sketch of the problem considered in this thesis. As seen in the figure, a horizontal well with transverse hydraulic fractures is placed within an inner reservoir zone. This inner zone has different reservoir properties than the outer zone. This inner zone is described as the stimulated reservoir volume (SRV); this thesis will assume that a stimulated reservoir is “stimulated” because of natural fractures that exist due to the hydraulic fracturing energy. This idea will be explained more in detail in Chapter 2.
This thesis uses dual porosity idealization to describe a naturally fractured reservoir, although other methods could be used. Dual porosity is an idealization of uniformly dis-tributed natural fractures on a continuum. If this assumption is not valid, another method, such as a discrete fracture network model, could be used. For the purposes of this thesis,
Figure 1.1: Example Case 1 - An example case diagram describing the layout of the inner and outer zone and their respective boundaries.
dual porosity was chosen because it is inherently simple for analytical models and has been used in the literature and reservoir engineering for decades and is well understood in the engineering community. This is beneficial because this research is intended to describe the behavior simply in well-understood methods so that it can be used and extrapolated to more complicated problems in future work.
Figure 1.2 shows another example case in which a transition zone is introduced. This transition zone will be stimulated also, but to a lesser extent than the inner zone. To clarify, the outer zone will be explored to be naturally fractured as well, however, the inner zone will be distinguished from the transition zone by a higher density of natural fractures. The distances to the boundaries shown in these figures in the x and y direction are arbitrary and several cases are explored with various distances. In the examples in Figure 1.1 and Figure 1.2, an arbitrary number of 10 fractures is shown. Neither this number nor any number of fractures in the rest of this thesis is intended to show an optimum well design but
rather is intended to explore the behavior of real world examples with correct physics.
Figure 1.2: Example Case 2 - Another example case diagram building off of Figure 1.1 with the addition of the transition zone and transition boundary.
Finally, it must be noted that the solutions are presented in the Laplace transform do-main. This approach is necessitated by the use of dual-porosity formulations and also pre-ferred by the elimination of the need for the discretization of the solutions in time. The results are obtained in the Laplace transform domain and inverted to the time domain by the Stehfest (1970a; 1970b) numerical inversion algorithm.
1.4 Contributions of the Study
The contribution of this research is primarily to provide a relatively simple but reasonably accurate solution to study flow behavior in naturally fractured, unconventional reservoirs. The primary focus of this revolves around accounting for the boundaries of composite re-gions, and how they affect the flow behavior and pressure distribution of multiply-fractured horizontal wells. Better understanding of the boundaries associated with these wells pro-vides a multitude of beneficial outcomes. First and foremost, this knowledge will help predict
long term flow performance, aiding in economic evaluation of projects and better estimated ultimate recovery (EUR) prediction. Second, understanding boundaries can aid in drilling programs and the desired well spacing to provide the most economic production acceleration and recovery factor of a field.
1.5 Organization of the Thesis
This thesis is organized into five chapters. The current chapter, Chapter 1, is the in-troduction chapter to the thesis. This chapter presents the purpose and motivation of this research, the approach and objectives intended, and the contribution to the field. Chapter 2 covers the background knowledge in the literature that lead up to this research. This includes the models from the literature which were adapted to build the composite-reservoir model for this thesis. Chapter 3 describes the development and the validation of the math-ematical model. Computational considerations encountered during the development of the model are also included in this chapter. Chapter 4 analyzes several cases pertinent to the reservoirs focused in the research. The last chapter, Chapter 5, presents the conclusions of the research and recommends directions for future work. And finally, Appendix A includes the FORTRAN code files created for this thesis
CHAPTER 2 LITERATURE REVIEW
This chapter covers some of the important background associated with unconventional reservoirs and current modeling approaches that leads up to this thesis. The first section covers properties of shales, and the subsequent sections cover the current state of analytical modeling for naturally fractured reservoirs and hydraulically fractured wells.
2.1 Shale Properties
An unconventional shale reservoir is defined as a reservoir system in which the formation itself is the source, seal, and reservoir. Most geologists would say that shale is composed of extremely fine-grained particles with high clay content, however many formations are named shales and do not have this attribute. Many shale reservoirs are fairly low in clay content, but were named so because of their dark color and negligible permeability compared to conventional reservoirs. The Bakken formation is a good example; it is primarily a low permeability dolomite, but it is named shale.
The most important factors for shale production are permeability, porosity, reservoir pressure, brittleness, kerogen content, and maturity. Permeability of shales are typically less than micro-Darcy range, and are typically measured in nano-Darcy. Porosity of shales can actually be quite high – as much as 30%; however, the majority of this is occupied by bound water, which leaves the effective porosity to be usually less than 10%. Shale matrix is composed of very fine grain particles and can be clay-rich, calcareous, or siliceous.
Because fractures, both natural and induced, are important for shale production, brittle-ness is a very important factor for shale development. A rock with a high Young’s modulus and a low Poisson’s ratio will be the most brittle and therefore easier to fracture. The miner-alogy of shale will define how brittle it is, and it is inversely proportional to the clay content. Productive shales usually contain a high kerogen content and therefore a high total organic
carbon (TOC) percentage. The maturity of the kerogen is also very important and is based on temperature, pressure, and time. Maturity is measured by the vitrinite reflectance (RO) with lower values indicating oil, medium values indicating gas, and high values indicating over maturity.
Some of the first gas wells in the mid-19th century were produced from shallow shale reservoirs. However, conventional gas resources in high permeability reservoirs became the primary target for gas production for nearly a century. Shale formations were considered to be seal and/or source rocks for conventional reservoirs until the early 1970’s with the OPEC oil embargo when the US pushed for energy independence. At that point in time, shale reservoirs were researched extensively and subsidized to find more reserves and production (Curtis, 2010).
Technology innovation was required to unlock the hydrocarbon potential of shales. Multiply-fractured horizontal wells have proven to be the best well design because it creates enormous surface area exposure to the reservoir allowing the drawdown in the reservoir to produce at economic rates. The first technically successful and expansive shale gas formation was the Barnett in Texas which has become the standard by which all other shales are measured. Following the Barnett was the Woodford in Oklahoma, Fayetteville in Arkansas, Haynesville in Louisiana and Texas, Marcellus in Pennsylvania and New York, as well as the shale oil play Bakken in North Dakota and Montana. Because of the success of these plays, newer shales such as the Eagleford in Texas, Niobrara in Colorado and Wyoming, Wolfcamp in Texas, Duverney in Alberta, Utica in Ohio as well as many others across the world are now being explored and developed. Figure 2.1 and Figure 2.2 show the shale plays in the United States and across the world, respectively. Many of these basins have existing shale produc-tion, and others have undeveloped shale resources that are exploitable and is being explored by many companies.
Figure 2.1: Shale basins in the United States with callouts for unconventional formations (Curtis, 2010).
2.2 Dual Porosity
Most models currently applied to shale reservoirs had been developed for different pur-poses long before the current completions in shale reservoirs. Modeling naturally fractured reservoirs started with Warren and Root (1963), who assumed pseudo-steady fluid transfer from matrix to fractures, and Kazemi (1969) who considered transient fluid transfer between matrix and fractures. Figure 2.3 from Warren and Root shows how the model is configured with cubic matrix blocks.
Figure 2.3: Warren and Root pictorial representation of an actual natural fractured reservoir and a cube modeled dual porosity reservoir (1963).
These two papers are the basis for most approaches to naturally fractured reservoirs. They define two parameters: storage capacity and flow capacity. In Laplace space, these models are the solutions of the diffusion equation which can be obtained from its counterpart for homogeneous reservoirs by replacing the Laplace transform parameter, s, with s · f (s), where f (s) is the transfer function which uses storativity, ω, and transmissivity, λ. Warren and Root (1963) defined the pseudo steady state transfer function as:
f (s) = ω(1 − ω) · s + λ
(1 − ω) · s + λ (2.1)
f (s) = 1 + s ˆ λˆω 3stanh s 3ˆωs ˆ λ (2.2)
The two models are different in their definitions of ω and λ. Brown (2009) examined the differences of the two models and discovered, for the purposes of unconventional reservoirs, the Kazemi transient model is much better at predicting reservoir pressures and produc-tion. In this case, the model assumes transient conditions within the matrix block. The matrix block sizes associated with unconventional reservoirs may take years for the pressure transients to reach the center of the matrix block, and therefore the pseudo-steady-state assumption of the Warren and Root model are not valid. Brown (2009) succinctly described the correct use of ˆω and ˆλ from Serra et al. (1983) for an alternating system of matrix and fracture slabs as:
ˆ ω = (φct)m (φct)f (2.3) ˆ λ = 12 l 2 h2 m ! kmhm kfhf ! (2.4) Other orientations of matrix blocks such as spheres, sticks, or parallelepipeds can be modeled with dual porosity through shape factors, but this discussion is outside of the scope of this thesis.
2.3 Finite-Conductivity Vertical Fracture Model
Interest in modeling vertically fractured wells started with the advent of fracture stimula-tion. Gringarten et al. (1975) and Cinco-Ley et al. (1978) presented models for infinite and finite-conductivity fractures, respectively. Figure 2.4 shows the configuration of the wellbore and fracture within the reservoir.
Cinco-Ley et al. defined a new parameter, dimensionless fracture conductivity. This parameter is described by:
Cf D =
kf · wf
km· xf
(a) 3-Dimensional representation of fracture model.
(b) Plan view representation with discretized fracture segments.
Figure 2.4: Geometry of Cinco-Ley et al. vertical fracture finite conductivity model (1978).
Where kf, wf, km, and xf are fracture permeability, fracture width, matrix permeability,
and fracture half length, respectively. Infinite conductivity is defined when Cf D ≥ 300
however a fracture exhibits near infinite conductivity when Cf D > 30. Figure 2.5 from
Cinco-Ley et al. (1978) shows this concept where above 30 Cf D, and certaintly above 300
Cf D, the curve is flat.
Cinco-Ley and Meng (1988) took the solution for finite conductivity fractures and put it in Laplace space. The benefit of this is that transfer functions, such as the dual porosity transfer function, can be applied easily. The model used is described in Cinco-Ley and Meng as ¯ pwd(s) − 1 4 ˆ xDf −xDf ¯ qf d(x0, s)K0(xD− x0) q sf (s)dx0+ π Cf D ˆ xD 0 ˆ x0 0 ¯ qf d(x”, s)dx”dx0 = πxD Cf D (2.6) The integral in this equation can be discretized as
Figure 2.5: Cinco-Ley et al. figure showing that dimensionless conductivities, labeled here as (kfbf)D, above 30 display near-infinite conductivity behavior (1978).
ˆ xDf −xDf ¯ qf d(x0, s)K0(xD − x0) q sf (s)dx0 = n X i=1 ¯ qf di(s) ˆ xDi+1 −xDi K0(xD − x0) q sf (s)dx0 (2.7)
2.4 Multiply-Fractured Horizontal Well Model
Raghavan et al. (1997) used the solution from Cinco-Ley and Meng (1988), and super-imposed multiple fractures next to each other to create a multiply-fractured horizontal well. Based on this model, Raghavan et al. described multiple flow regimes shown in Figure 2.6.
As described in their paper, during the late-time pseudo-radial flow, the fractured hori-zontal well may be thought as an equivalent hydraulic fracture for which the length is equal to the distance between the outermost fractures. Raghavan et al. expanded Figure 2.5 above from Cinco-Ley et al. (1978) to incorporate multiple fractures shown in Figure 2.7. This concept becomes very important because as more fractures are added, production enhance-ment is diminished due to interference between fractures. Examining the well orientation and
Figure 2.6: Raghavan et al.’s nomenclature for the four flow regimes in a multiply-fractured horizontal well (1997).
reservoir properties in shale reservoirs, the early radial flow regime is eliminated because the spacing between fractures is very close and therefore the fractures start to interfere with each other in early times. Also late-time radial flow can only occur many years, even decades, after initial production. Eliminating these two radial flow regimes leaves only three flow regimes – very early-time linear flow within the fracture, early-time linear flow in the reservoir which is perpendicular to the fracture, and late-time linear flow beyond the tips of the fractures, which is parallel to the fracture (referred to as “compound linear flow” by Raghavan et al.).
Figure 2.7: Raghaven et al.’s expanded Cinco-Ley et al. figure showing diminishing returns for adding additional fractures (1997).
2.5 Trilinear Model
For multiply-fractured horizontal wells in ultra-low permeability reservoirs, three linear flow regimes dominate the well’s lifetime production. This behavior of linear flow regimes is clearly described by Brown (2009). Figure 2.8 shows the description of three flow regimes that are used to create the trilinear model. In addition to the concept of three, mutually
perpendicular, linear flows, Brown also considered the stimulated reservoir volume (SRV) between hydraulic fractures. During the stimulation of the well, enough energy is put into the reservoir that the micro- and macro-natural fractures become activated and open. Because the stimulation opens up natural fractures, the inner zone between hydraulic fractures can be modeled as a dual porosity reservoir. The outer reservoir, which is not affected by the stimulation, is then modeled as a single porosity system.
Figure 2.8: Brown and Ozkan’s trinlinear model geometry (2009).
Using this trilinear model, Ozkan et al. (2009) was able to make several observations. Production is dominated to the inner reservoir in low matrix permeability reservoirs causing a near-no-flow boundary between the inner and outer reservoirs. Natural fracture density is more important than natural fracture permeability. Hydraulic fracture conductivity should be balanced with the flow capacity of the matrix and natural fractures. Decreasing hydraulic fracture spacing along the wellbore increases productivity albeit with diminishing returns with each increase of fractures. These observations lead to the conclusion that wellbore and
completion designs should be optimized to create complex hydraulic fractures as opposed to high conductivity planar fractures. In addition, spacing of the hydraulic stimulation should be close enough to stimulate the natural fracture network however not too close to create early time interference.
CHAPTER 3
MATHEMATICAL MODEL
This chapter presents the development and verification of the approximate mathematical model by superimposing independent solutions corresponding to flow in various reservoir sections as introduced in Chapter 1. The fundamental premise of the heuristic approach used here is that the solution for a composite reservoir, which consists of an SRV enclosing a multiply-fractured horizontal well and surrounded by a virgin reservoir shown in Fig-ure 1.1, can be derived by superimposing three solutions: Solution 1 – a multiply fractFig-ured, horizontal-well in an infinite-acting, homogeneous reservoir with the properties of the outer zone; Solution 2 – a multiply fractured, horizontal-well in a bounded, homogeneous (un-fractured) reservoir with the properties of the outer reservoir; and Solution 3 – a multiply fractured, horizontal-well in a bounded, naturally fractured reservoir with the properties of the stimulated zone. The solution for the composite reservoir consisting of a stimulated (nat-urally fractured) reservoir surrounded by an infinite acting, un-fractured (virgin) reservoir is, then, given by
Composite System Solution = Solution 1–Solution 2 + Solution 3 (3.1) The same approach is also applied to develop a solution for the case where there is an addi-tional transition zone between the SRV and the outer (virgin) reservoir shown in Figure 1.2 as follows
Composite System Solution =
Solution 1–Solution 2 + Solution 3–Solution 4 + Solution 5 (3.2) where the solutions are defined as follows: Solution 1 – a multiply fractured, horizontal-well in an infinite-acting, homogeneous reservoir with the properties of the outer zone; Solution 2
with the properties of the outer reservoir but the size of the transition region; Solution 3 – a multiply fractured, horizontal-well in a bounded, homogeneous (un-fractured) reservoir with the properties and the size of the transition region; Solution 4 – a multiply fractured, horizontal-well in a bounded, homogeneous (un-fractured) reservoir with the properties of the transition region but the size of the SRV; and Solution 5 – a multiply fractured, horizontal-well in a bounded, naturally fractured reservoir with the properties and the size of the SRV. We note that the solution is derived in the Laplace transform domain to avoid discretization of the solution in time for numerical evaluations and for the convenience of incorporating dual-porosity formulations in Laplace domain.
3.1 Details of the Model Development
The following sections describe the details of the development of the solution and the procedure to calculate it. As customary, the solution is derived in terms of dimensionless variables in this work. Therefore, first, the dimensionless variables are defined. Then the individual solutions to be superimposed are introduced and the computational procedure is discussed. These discussions include infinite and bounded reservoir solutions for hydraulically fractured wells, extensions to multiple, finite-conductivity fractures along a horizontal wells, and superimposing of the appropriate infinite and bounded reservoir solutions to obtain the composite-reservoir solution.
3.1.1 Dimensionless Variables
The standard definitions of dimensionless pressure and time used in petroleum engi-neering use the reservoir permeability. Because the solution procedure involves superposing individual solutions representing different flow zones of varying permeability, a common reference-permeability is used for the definition of the dimensionless variables. This intro-duces some scaling parameters into the dimensionless diffusion equations for each flow zone. For demonstration purposes, here we will use the application of dimensionless variables to the simple case of 1D diffusion equation,
∂2∆p ∂x2 = φctµ k ! ∂∆p ∂t (3.3)
with the definition of the following dimensionless variables pD = kh 141.2qoBoµ ∆p (3.4) xD = x xf (3.5) tD = η x2 f t (3.6) where η = 2.637 (10)−4 k φctµ0 (3.7) Inserting Equations 3.4, 3.5, and 3.6 into Equation 3.3 gives:
∂2pD ∂x2 D = ∂pD ∂tD (3.8) and in Laplace space is
∂2p¯D
∂x2 D
= spD(s) (3.9)
Now, take the simple version of Figure 1.1 where the only difference between the inner zone and the outer zone is two different permeabilties, k1 and k2. Let us select k = k1 as the
reference permeability to define dimensionless variables and note that
ηD = 2.637 (10)−4 k1 φctµ 2.637 (10)−4 k2 φctµ = k1 k2 (3.10) Then, we have for the inner and the outer zone, respectively,
∂2p¯D1 ∂x2 D = s ¯pD1 (3.11) and ∂2p¯D2 ∂x2 D = sηDp¯D2 (3.12)
3.1.2 Development and Calculation of the Composite-Reservoir Solution To calculate the solutions described in Equation 3.1 we must now start with the derivation described in Ozkan (1988). The pressure solution in Laplace time for a fully penetrating vertical fracture in an infinite reservoir is
¯ pD = 1 2s +xˆf D −xf D ¯ qf DK0 q (xD− xwD− α)2+ (yD − ywD)2 √ u dα (3.13)
As described by Cinco-Ley and Meng (1988), with the assumption that the flux, ¯qf D is
uniform over the surface of the fracture, this solution can be rearranged to yield
¯ pD = 1 2sqf D¯ +xˆf D −xf D K0 q (xD− xwD− α)2+ (yD − ywD)2 √ u dα (3.14)
Cinco-Ley and Meng have also introduced a discretized form of the solution by assuming that flux is uniform in each fracture segment but allowed to vary from one-segment to the next. This discretized solution can be used to obtain infinite- and finite-conductivity fracture results and is given by
¯ pD = 1 2s n X i=1 ¯ qf Di +xˆf Di+1 −xf Di K0 q (xDj− xwD− α)2+ (yDj− ywD)2 √ u dα (3.15)
Writing this equation for each segment j creates a system of n equations with n + 1 unknowns. The n + 1 equation is obtained from the condition that the sum of the fluxes equals the total flowrate of the well. In Laplace space this is described as
n X i=1 ¯ qf Di= 1 s (3.16)
For infinite-conductivity fractures, where the pressure is the same in all fracture segments and equal to the wellbore pressure, pwf, this formulation leads to the following system of
equations: [A] " ¯ pD ¯ qf Di # = [B] (3.17)
where the A matrix is n + 1 byn + 1 matrix filled in a way in which A1,1 to n= 1 (3.18) A1,n+1 = 0 (3.19) A2 to n+1,n+1= 1 (3.20) A2−n,1−n = − 1 2 +xˆf Di+1 −xf Di K0 q (xDj− xwD− α) 2 + (yDj− ywD) 2√ u dα (3.21)
and the B matrix is an n + 1 by 1 size matrix filled in a way that
B1 to n= 0 (3.22)
Bn+1=1/s (3.23)
Cinco-Ley and Meng’s model incorporated fracture segments only in the x direction in which the term (yDj− ywD)
2
is zero. Raghavan et al. (1997) first updated this model to incorporate multiple fractures in the y direction. This created the multiply-fracture horizontal well model. This is easier described in a symmetric case where m is the number of fractures and n is the number of fracture segments per fracture.
¯ pD = 1 2s m X k=1 n X i=1 ¯ qf Di,k +xf Di+1,kˆ −xf Di,k K0 q (xDj,k − xwD− α)2 + (yDj,k− ywD)2 √ u dα (3.24) However, the model does not need to be symmetric, each fracture can be modelled with different dimensions; it just needs to be described in each equation in the system of equations. This thesis will assume symmetric fractures unless otherwise noted. This equation is able to be solved in the same manner with the form of Equation 3.17 with the size of the A matrix now (n ∗ m) + 1 by (n ∗ m) + 1 and the size of the B matrix now (n ∗ m) + 1 by 1.
¯ pD Bounded = π ¯qf D xeDs " cosh (√u (yeD− |yD− ywD|)) + cosh ( √ u (yeD− |yD+ ywD|)) √ usinh (√uyeD) # +2xeDqf D¯ πs ∞ X k=1 1 ksin kπ 1 xeD cos kπxwD xeD cos kπ xd xeD cosh r u +kx22π2 eD (yeD− |yD − ywD|) + cosh r u + kx22π2 eD (yeD− |yD+ ywD|) r u +kx22π2 eD sinh r u +kx22π2 eD yeD (3.25)
As explained by Ozkan, the bounded reservoir solution is cumbersome and difficult to calculate with the hyperbolic functions and infinite sums. To improve the speed and accuracy of the numerical calculations at early times, Ozkan expressed the solution in the following form: ¯ pD Bounded= ¯pD1+ ¯pDb1+ ¯pDb2 (3.26) where ¯ pD1 = 2 ¯qf D s ∞ X k=1 1 ksin kπ 1 xeD cos kπxD xeD cos kπxwD xeD exp−√u + a |yD − ywD| √ u + a (3.27) ¯ pDb1= π ¯qf D xeDs √ u n exph−√u (yD+ ywD) i + exph−√u (2yeD− yD+ ywD) i +exph−√u |yD− ywD| i + exph−√u (2yeD − |yD − ywD|) io
( 1 + ∞ X m=1 h exp−2m√uyeD i ) (3.28) and ¯ pDb2= 2 ¯qf D s ∞ X k=1 1 k sinkπ 1 xeD coskπxwD xeD coskπ xd xeD √ u + a "( exph−√u + a (yD+ ywD) i + exph−√u + a (2yeD− yD+ ywD) i +exph−√u + a (2yeD− |yD− ywD|) i ) ( 1 + ∞ X m=1 h exp−2m√u + ayeD i ) +exph−√u |yD − ywD| i ( 1 + ∞ X m=1 h exp−2m√u + ayeD i )# (3.29) where a = k 2π2 x2 eD (3.30) Similarly, Ozkan suggested the following, numerically more efficient form of the bounded-system solution for late times:
¯
pD Bounded =pDinf¯ + ¯pDb1+ ¯pDb2+ ¯pDb3 (3.31)
where pDinf¯ is is the infinite reservoir solution given by equation 3.15 and pDb1¯ , and pDb2¯
are the solutions given by Equations 3.28 and 3.29, respectively. The termpDb3¯ in Equation 3.31 is given by ¯ pDb3 = 1 2s ˆ +xf D −xf D ¯ qf DK0 q (xD + xwD− α)2+ (yD − ywD)2 √ u dα
+1 2s ∞ X k=1 (ˆ +xf D −xf D ¯ qf D K0 q (xD− xwD− 2kxeD − α) 2 + (yD− ywD) 2√ u +K0 q (xD + xwD− 2kxeD− α)2+ (yD − ywD)2 √ u +K0 q (xD − xwD+ 2kxeD− α)2+ (yD − ywD)2 √ u +K0 q (xD+ xwD+ 2kxeD − α)2+ (yD− ywD)2 √ u dα −π exp (− √ u |yD − ywD|) xeDs √ u (3.32)
If evaluating on the fracture surface, (yD = ywD), then pDb3¯ can be simplified further to
¯ pDb3= 1 2s√u ˆ √ u(xD+xwD+xf D) 0 K0(z) dz − ˆ √ u(xD+xwD−xf D) 0 K0(z) dz ∞ X k=1 ˆ √ u(2kxeD−xD+xwD+xf D) 0 K0(z) dz − ˆ √ u(2kxeD−xD+xwD−xf D) 0 K0(z) dz ˆ √ u(2kxeD+xD−xwD+xf D) 0 K0(z) dz − ˆ √ u(2kxeD+xD−xwD−xf D) 0 K0(z) dz ˆ √ u(2kxeD−xD−xwD+xf D) 0 K0(z) dz − ˆ √ u(2kxeD−xD−xwD−xf D) 0 K0(z) dz ˆ √ u(2kxeD+xD+xwD+xf D) 0 K0(z) dz − ˆ √ u(2kxeD+xD+xwD−xf D) 0 K0(z) dz − π xeDs √ u (3.33)
The elegance in stating the model in the format described in Equation 3.26 and 3.31 is that solving of a system of fractures through Equation 3.17, A2−n,1−n may now be replaced
The next addition to the model is the consideration of finite fracture conductivity. Cinco-Ley and Meng (1988)derived the following solution for a finite-conductivity fracture in an infinite reservoir: ¯ pD −pDinf¯ + π Cf D xD ˆ xwD x0 ˆ 0 [ ¯qf Di(x”, s)] dx” dx0 = πxD Cf Ds (3.34) where thepDinf¯ term is the infinite system solution given by Equation 3.13. The pDinf¯ term
can be discretized as in Equation 3.15 and the double integral in the third term on the left side of Equation 3.34 can be discretized as follows:
π Cf D xD ˆ xwD x0 ˆ 0 [ ¯qf Di(x”, s)] dx” dx0 = π Cf D j−1 X i=1 ( ¯ qf Di " (∆x)2 2 + ∆x (xDj− i∆x) #) +(∆x) 2 8 qf Dj¯ (3.35)
The discretized equation can be evaluated at the center of each fracture segment and written in the form of Equation3.17. The evaluation of the solution follows the same lines as the evaluation of the infinite-conductivity solution discussed above. It must also be noted that replacing pDinf¯ in Equation 3.34 with pDbounded¯ given by Equations 3.26 and 3.31 leads
to the counterpart of the solution for a bounded, rectangular reservoir.
As noted earlier, the solution of the matrix problem is numerically inverted from Laplace space to time domain. For this thesis, Stehfest’s (1970a; 1970b) numerical inversion algorithm has been used. This algorithm expresses the numerical inversion, f (t), of ¯F at time t by:
f (t) ≈ ln(2) t N X i=1 ViF¯ i · ln(2) t ! (3.36) where Vi = (−1)( N/2+i) min(i,N/2) X k=(i+1)/2 " kN/2(2k)! (N/2− k)! k! (k − 1)! (i − k)! (2k − i)! # (3.37)
3.2 Model Verification
The process of building the code required several steps of model verification. The model was built by starting with the simplest of cases, a single-segment fracture in a bounded reservoir, and from there, every additional complexity added to the model was checked against the examples in the literature to confirm accuracy. Because no analytical/semi-analytical solution is available for the final composite-system solution in the literature, the comparison was made against a numerical model. In all of the verification plots in this section, the lines indicate the model output, and the symbols correspond to the literature or numerical simulation data used in the verification.
3.2.1 Uniform-Flux (Single-Segment) Fracture Solution in a Bounded Reservoir Ozkan (1988) provided the dimensionless pressures for an example case of uniform-flux (single-segment) fracture in a bounded, rectangular reservoir. This solution was useful for the verifications because the multi-segment (finite-conductivity) fracture code needed for the final solution uses the uniform flux solution for each fracture segment. Figure 3.1 and Table 3.1 show the error associated with this verification. In this figure dimensionless time is shown in terms of TDA which is defined as a function of the area of the bounded region.
This is described as TDA = TD AD (3.38) where AD = A x2 f (3.39) In this example the boundary edges, xeD and yeD, are equal to create a square boundary.
As can be seen in the figure and table, the four outputs are near identical for the literature and model output. Therefore, the error associated is acceptable to conclude that the code is running correctly.
Table 3.1: Single fracture segment model verification against Ozkan (1988) Table 2.6.1.
TDA
Model Output Ozkan Table 2.6.1 Output
pD pD pD pD pD pD pD pD
4 AD 16 AD 100 AD 400 AD 4 AD 16 AD 100 AD 400 AD
1.0E-4 3.54E-2 7.09E-2 1.77E-1 3.54E-1 3.54E-2 7.09E-2 1.77E-1 3.55E-1 1.5E-4 4.34E-2 8.68E-2 2.17E-1 4.34E-1 4.34E-2 8.68E-2 2.17E-1 4.34E-1 2.0E-4 5.01E-2 1.00E-1 2.51E-1 5.01E-1 5.01E-2 1.00E-1 2.51E-1 5.01E-1 3.0E-4 6.14E-2 1.23E-1 3.07E-1 6.11E-1 6.14E-2 1.23E-1 3.07E-1 6.11E-1 4.0E-4 7.09E-2 1.42E-1 3.54E-1 7.00E-1 7.09E-2 1.42E-1 3.55E-1 7.00E-1 5.0E-4 7.93E-2 1.59E-1 3.96E-1 7.76E-1 7.93E-2 1.59E-1 3.96E-1 7.76E-1 6.0E-4 8.68E-2 1.74E-1 4.34E-1 8.41E-1 8.68E-2 1.74E-1 4.34E-1 8.41E-1 8.0E-4 1.00E-1 2.01E-1 5.01E-1 9.51E-1 1.00E-1 2.01E-1 5.01E-1 9.52E-1 1.0E-3 1.12E-1 2.24E-1 5.59E-1 1.04E+0 1.12E-1 2.24E-1 5.59E-1 1.04E+0 1.5E-3 1.37E-1 2.75E-1 6.79E-1 1.21E+0 1.37E-1 2.75E-1 6.79E-1 1.21E+0 2.0E-3 1.59E-1 3.17E-1 7.76E-1 1.34E+0 1.59E-1 3.17E-1 7.76E-1 1.34E+0 3.0E-3 1.94E-1 3.88E-1 9.26E-1 1.53E+0 1.94E-1 3.88E-1 9.26E-1 1.53E+0 4.0E-3 2.24E-1 4.48E-1 1.04E+0 1.66E+0 2.24E-1 4.48E-1 1.04E+0 1.67E+0 5.0E-3 2.51E-1 5.01E-1 1.14E+0 1.77E+0 2.51E-1 5.01E-1 1.14E+0 1.77E+0 6.0E-3 2.75E-1 5.48E-1 1.21E+0 1.86E+0 2.75E-1 5.48E-1 1.21E+0 1.86E+0 8.0E-3 3.17E-1 6.30E-1 1.34E+0 2.00E+0 3.17E-1 6.30E-1 1.34E+0 2.00E+0 1.0E-2 3.54E-1 7.00E-1 1.44E+0 2.11E+0 3.55E-1 7.00E-1 1.44E+0 2.11E+0 1.5E-2 4.34E-1 8.41E-1 1.63E+0 2.31E+0 4.34E-1 8.41E-1 1.63E+0 2.31E+0 2.0E-2 5.01E-1 9.52E-1 1.77E+0 2.45E+0 5.01E-1 9.52E-1 1.77E+0 2.45E+0 3.0E-2 6.14E-1 1.12E+0 1.97E+0 2.65E+0 6.14E-1 1.12E+0 1.97E+0 2.65E+0 4.0E-2 7.09E-1 1.24E+0 2.11E+0 2.79E+0 7.09E-1 1.24E+0 2.11E+0 2.79E+0 5.0E-2 7.93E-1 1.35E+0 2.22E+0 2.91E+0 7.94E-1 1.35E+0 2.22E+0 2.91E+0 6.0E-2 8.70E-1 1.43E+0 2.31E+0 3.00E+0 8.71E-1 1.44E+0 2.31E+0 3.00E+0 8.0E-2 1.01E+0 1.59E+0 2.47E+0 3.16E+0 1.01E+0 1.59E+0 2.47E+0 3.16E+0 1.0E-1 1.15E+0 1.73E+0 2.61E+0 3.30E+0 1.15E+0 1.73E+0 2.61E+0 3.30E+0 1.5E-1 1.47E+0 2.05E+0 2.94E+0 3.62E+0 1.47E+0 2.05E+0 2.94E+0 3.63E+0 2.0E-1 1.78E+0 2.36E+0 3.25E+0 3.94E+0 1.78E+0 2.37E+0 3.25E+0 3.94E+0 3.0E-1 2.41E+0 2.99E+0 3.88E+0 4.57E+0 2.41E+0 2.99E+0 3.88E+0 4.57E+0 4.0E-1 3.04E+0 3.62E+0 4.51E+0 5.20E+0 3.04E+0 3.62E+0 4.51E+0 5.20E+0 5.0E-1 3.67E+0 4.25E+0 5.14E+0 5.82E+0 3.67E+0 4.25E+0 5.14E+0 5.83E+0 6.0E-1 4.29E+0 4.88E+0 5.76E+0 6.45E+0 4.29E+0 4.88E+0 5.77E+0 6.46E+0 8.0E-1 5.55E+0 6.13E+0 7.02E+0 7.71E+0 5.55E+0 6.14E+0 7.02E+0 7.71E+0 1.0E+0 6.81E+0 7.39E+0 8.28E+0 8.97E+0 6.81E+0 7.39E+0 8.28E+0 8.97E+0 1.5E+0 9.95E+0 1.05E+1 1.14E+1 1.21E+1 9.95E+0 1.05E+1 1.14E+1 1.21E+1 2.0E+0 1.31E+1 1.37E+1 1.46E+1 1.52E+1 1.31E+1 1.37E+1 1.46E+1 1.53E+1 3.0E+0 1.94E+1 2.00E+1 2.08E+1 2.15E+1 1.94E+1 2.00E+1 2.08E+1 2.15E+1 4.0E+0 2.57E+1 2.62E+1 2.71E+1 2.78E+1 2.57E+1 2.62E+1 2.71E+1 2.78E+1
Figure 3.1: Single fracture segment model verification against Ozkan (1988) Table 2.6.1.
3.2.2 Finite-Conductivity (Multiple Segments), Single-Fracture Solution
The next complexity to add to the model is the effect of finite fracture conductivity (multiple fracture segments) in a single fracture. To verify this model, Cinco-Ley et al.’s (1978) Table 1 example was used. Cinco-Ley et al. did not explicitly provide the number of fracture segments used in this example, although they stated that the solutions did not change appreciably when more than 20 fracture segments are used. However, our computa-tions indicated that the results in their Table 1 were obtained with more than 20 segments (probably 40 fracture segments). Because of this, the results from the model developed in this study for 40 fracture segments were compared to those from the study of Cinco-Ley et al. are shown in Figure 3.2 and Table 3.2, Table 3.3, and Table 3.4. See Section 3.3.1 for more detail on the number of fracture segments. The error is extremely low in all six cases compared in Table 3.2, Table 3.3, and Table 3.4.
Table 3.2: Finite conductivity, single fracture, multiple segment model verification against Cinco-Ley et al. (1978) Table 1.
TDA
Model Output Cinco-Ley et al. Table 1 Output
pD pD pD pD pD pD
0.2π Cf D π Cf D 2π Cf D 0.2π Cf D π Cf D 2π Cf D
1.E-3 5.41E-1 2.45E-1 1.74E-1 5.45E-1 2.44E-1 1.73E-1 2.E-3 6.38E-1 2.91E-1 2.07E-1 6.38E-1 2.88E-1 2.06E-1 3.E-3 7.02E-1 3.22E-1 2.30E-1 7.02E-1 3.19E-1 2.29E-1 4.E-3 7.51E-1 3.45E-1 2.49E-1 7.52E-1 3.43E-1 2.48E-1 5.E-3 7.92E-1 3.65E-1 2.64E-1 7.93E-1 3.63E-1 2.63E-1 6.E-3 8.26E-1 3.82E-1 2.78E-1 8.27E-1 3.81E-1 2.77E-1 7.E-3 8.56E-1 3.97E-1 2.90E-1 8.58E-1 3.96E-1 2.89E-1 8.E-3 8.83E-1 4.11E-1 3.01E-1 8.85E-1 4.10E-1 3.01E-1 9.E-3 9.07E-1 4.23E-1 3.11E-1 9.09E-1 4.22E-1 3.11E-1 1.E-2 9.29E-1 4.35E-1 3.21E-1 9.31E-1 4.34E-1 3.21E-1 2.E-2 1.09E+0 5.21E-1 3.95E-1 1.08E+0 5.18E-1 3.94E-1 3.E-2 1.19E+0 5.81E-1 4.49E-1 1.19E+0 5.79E-1 4.48E-1 4.E-2 1.26E+0 6.29E-1 4.92E-1 1.26E+0 6.27E-1 4.92E-1 5.E-2 1.33E+0 6.69E-1 5.30E-1 1.33E+0 6.68E-1 5.30E-1 6.E-2 1.38E+0 7.05E-1 5.63E-1 1.38E+0 7.04E-1 5.63E-1 7.E-2 1.43E+0 7.36E-1 5.92E-1 1.43E+0 7.36E-1 5.93E-1 8.E-2 1.47E+0 7.65E-1 6.20E-1 1.47E+0 7.65E-1 6.20E-1 9.E-2 1.50E+0 7.92E-1 6.45E-1 1.50E+0 7.92E-1 6.46E-1 1.E-1 1.54E+0 8.17E-1 6.68E-1 1.54E+0 8.17E-1 6.69E-1 2.E-1 1.77E+0 1.00E+0 8.45E-1 1.77E+0 1.00E+0 8.45E-1 3.E-1 1.92E+0 1.13E+0 9.68E-1 1.93E+0 1.13E+0 9.69E-1 4.E-1 2.04E+0 1.23E+0 1.06E+0 2.04E+0 1.23E+0 1.06E+0 5.E-1 2.13E+0 1.31E+0 1.14E+0 2.13E+0 1.31E+0 1.14E+0 6.E-1 2.21E+0 1.38E+0 1.21E+0 2.21E+0 1.38E+0 1.21E+0 7.E-1 2.27E+0 1.44E+0 1.27E+0 2.28E+0 1.44E+0 1.27E+0 8.E-1 2.33E+0 1.49E+0 1.32E+0 2.34E+0 1.49E+0 1.32E+0 9.E-1 2.38E+0 1.54E+0 1.37E+0 2.39E+0 1.54E+0 1.37E+0 1.E+0 2.43E+0 1.58E+0 1.41E+0 2.44E+0 1.59E+0 1.42E+0 2.E+0 2.75E+0 1.89E+0 1.71E+0 2.76E+0 1.89E+0 1.72E+0 3.E+0 2.95E+0 2.08E+0 1.90E+0 2.95E+0 2.08E+0 1.90E+0 4.E+0 3.09E+0 2.21E+0 2.03E+0 3.09E+0 2.21E+0 2.04E+0 5.E+0 3.19E+0 2.32E+0 2.14E+0 3.20E+0 2.32E+0 2.14E+0 6.E+0 3.28E+0 2.41E+0 2.23E+0 3.29E+0 2.41E+0 2.23E+0 7.E+0 3.36E+0 2.48E+0 2.30E+0 3.37E+0 2.48E+0 2.31E+0 8.E+0 3.43E+0 2.55E+0 2.37E+0 3.43E+0 2.55E+0 2.37E+0 9.E+0 3.48E+0 2.60E+0 2.42E+0 3.49E+0 2.61E+0 2.43E+0
Table 3.3: Finite conductivity, single fracture, multiple segment model verification against Cinco-Ley et al. (1978) Table 1. (continued)
TDA
Model Output Cinco-Ley et al. Table 1 Output
pD pD pD pD pD pD
0.2π Cf D π Cf D 2π Cf D 0.2π Cf D π Cf D 2π Cf D
2.E+1 3.88E+0 3.00E+0 2.82E+0 3.89E+0 3.00E+0 2.82E+0 3.E+1 4.08E+0 3.20E+0 3.02E+0 4.07E+0 3.20E+0 3.02E+0 4.E+1 4.22E+0 3.34E+0 3.16E+0 4.23E+0 3.34E+0 3.16E+0 5.E+1 4.34E+0 3.45E+0 3.27E+0 4.34E+0 3.45E+0 3.28E+0 6.E+1 4.43E+0 3.54E+0 3.36E+0 4.43E+0 3.55E+0 3.37E+0 7.E+1 4.50E+0 3.62E+0 3.44E+0 4.51E+0 3.62E+0 3.44E+0 8.E+1 4.57E+0 3.69E+0 3.51E+0 4.58E+0 3.69E+0 3.51E+0 9.E+1 4.63E+0 3.74E+0 3.56E+0 4.64E+0 3.75E+0 3.57E+0 1.E+2 4.68E+0 3.80E+0 3.62E+0 4.69E+0 3.80E+0 3.62E+0 2.E+2 5.03E+0 4.14E+0 3.96E+0 5.03E+0 4.15E+0 3.97E+0 3.E+2 5.23E+0 4.35E+0 4.16E+0 5.24E+0 4.35E+0 4.17E+0 4.E+2 5.37E+0 4.49E+0 4.31E+0 5.38E+0 4.49E+0 4.31E+0 5.E+2 5.49E+0 4.60E+0 4.42E+0 5.49E+0 4.60E+0 4.42E+0 6.E+2 5.58E+0 4.69E+0 4.51E+0 5.58E+0 4.70E+0 4.52E+0 7.E+2 5.65E+0 4.77E+0 4.59E+0 5.66E+0 4.77E+0 4.59E+0 8.E+2 5.72E+0 4.84E+0 4.66E+0 5.73E+0 4.84E+0 4.66E+0 9.E+2 5.78E+0 4.89E+0 4.71E+0 5.79E+0 4.90E+0 4.72E+0
TDA
Model Output Cinco-Ley et al. Table 1 Output
pD pD pD pD pD pD
10π Cf D 20π Cf D 100π Cf D 10π Cf D 20π Cf D 100π Cf D
1.E-3 8.66E-2 7.17E-2 5.88E-2 8.66E-2 7.18E-2 5.90E-2 2.E-3 1.10E-1 9.45E-2 8.14E-2 1.10E-1 9.46E-2 8.14E-2 3.E-3 1.28E-1 1.12E-1 9.85E-2 1.28E-1 1.12E-1 9.86E-2 4.E-3 1.42E-1 1.26E-1 1.13E-1 1.42E-1 1.27E-1 1.13E-1 5.E-3 1.55E-1 1.39E-1 1.25E-1 1.55E-1 1.39E-1 1.26E-1 6.E-3 1.67E-1 1.50E-1 1.37E-1 1.67E-1 1.51E-1 1.37E-1 7.E-3 1.77E-1 1.61E-1 1.47E-1 1.77E-1 1.61E-1 1.47E-1 8.E-3 1.87E-1 1.70E-1 1.56E-1 1.87E-1 1.71E-1 1.57E-1 9.E-3 1.96E-1 1.79E-1 1.65E-1 1.96E-1 1.80E-1 1.66E-1 1.E-2 2.04E-1 1.87E-1 1.74E-1 2.05E-1 1.88E-1 1.74E-1 2.E-2 2.71E-1 2.54E-1 2.40E-1 2.72E-1 2.55E-1 2.41E-1 3.E-2 3.21E-1 3.04E-1 2.89E-1 3.22E-1 3.05E-1 2.90E-1 4.E-2 3.62E-1 3.45E-1 3.30E-1 3.63E-1 3.46E-1 3.31E-1 5.E-2 3.98E-1 3.80E-1 3.65E-1 3.97E-1 3.81E-1 3.66E-1 6.E-2 4.29E-1 4.11E-1 3.96E-1 4.30E-1 4.12E-1 3.94E-1 7.E-2 4.57E-1 4.39E-1 4.24E-1 4.59E-1 4.41E-1 4.26E-1
Table 3.4: Finite conductivity, single fracture, multiple segment model verification against Cinco-Ley et al. (1978) Table 1. (continued)
TDA
Model Output Cinco-Ley et al. Table 1 Output
pD pD pD pD pD pD
10π Cf D 20π Cf D 100π Cf D 10π Cf D 20π Cf D 100π Cf D
8.E-2 4.83E-1 4.65E-1 4.50E-1 4.85E-1 4.69E-1 4.51E-1 9.E-2 5.07E-1 4.89E-1 4.74E-1 5.09E-1 4.91E-1 4.75E-1 1.E-1 5.30E-1 5.11E-1 4.96E-1 5.32E-1 5.13E-1 4.98E-1 2.E-1 7.00E-1 6.80E-1 6.64E-1 7.02E-1 6.82E-1 6.66E-1 3.E-1 8.19E-1 7.99E-1 7.82E-1 8.21E-1 8.01E-1 7.85E-1 4.E-1 9.12E-1 8.92E-1 8.75E-1 9.14E-1 8.94E-1 8.77E-1 5.E-1 9.89E-1 9.68E-1 9.52E-1 9.92E-1 9.71E-1 9.54E-1 6.E-1 1.06E+0 1.03E+0 1.02E+0 1.06E+0 1.04E+0 1.02E+0 7.E-1 1.11E+0 1.09E+0 1.08E+0 1.12E+0 1.10E+0 1.08E+0 8.E-1 1.16E+0 1.14E+0 1.13E+0 1.17E+0 1.15E+0 1.13E+0 9.E-1 1.21E+0 1.19E+0 1.17E+0 1.22E+0 1.19E+0 1.18E+0 1.E+0 1.25E+0 1.23E+0 1.22E+0 1.26E+0 1.24E+0 1.22E+0 2.E+0 1.55E+0 1.53E+0 1.51E+0 1.55E+0 1.53E+0 1.52E+0 3.E+0 1.74E+0 1.71E+0 1.70E+0 1.74E+0 1.72E+0 1.70E+0 4.E+0 1.87E+0 1.85E+0 1.83E+0 1.87E+0 1.85E+0 1.83E+0 5.E+0 1.98E+0 1.95E+0 1.93E+0 1.98E+0 1.96E+0 1.94E+0 6.E+0 2.06E+0 2.04E+0 2.02E+0 2.07E+0 2.05E+0 2.03E+0 7.E+0 2.14E+0 2.11E+0 2.10E+0 2.14E+0 2.12E+0 2.10E+0 8.E+0 2.20E+0 2.18E+0 2.16E+0 2.21E+0 2.18E+0 2.17E+0 9.E+0 2.26E+0 2.24E+0 2.22E+0 2.26E+0 2.24E+0 2.22E+0 1.E+1 2.31E+0 2.29E+0 2.27E+0 2.31E+0 2.29E+0 2.27E+0 2.E+1 2.65E+0 2.63E+0 2.61E+0 2.66E+0 2.63E+0 2.62E+0 3.E+1 2.85E+0 2.83E+0 2.81E+0 2.86E+0 2.83E+0 2.82E+0 4.E+1 2.99E+0 2.97E+0 2.95E+0 3.00E+0 2.98E+0 2.96E+0 5.E+1 3.10E+0 3.08E+0 3.06E+0 3.11E+0 3.09E+0 3.07E+0 6.E+1 3.20E+0 3.17E+0 3.15E+0 3.20E+0 3.18E+0 3.16E+0 7.E+1 3.27E+0 3.25E+0 3.23E+0 3.28E+0 3.26E+0 3.24E+0 8.E+1 3.34E+0 3.32E+0 3.30E+0 3.34E+0 3.32E+0 3.30E+0 9.E+1 3.40E+0 3.38E+0 3.36E+0 3.40E+0 3.38E+0 3.36E+0 1.E+2 3.45E+0 3.43E+0 3.41E+0 3.46E+0 3.43E+0 3.41E+0 2.E+2 3.80E+0 3.77E+0 3.76E+0 3.80E+0 3.78E+0 3.76E+0 3.E+2 4.00E+0 3.98E+0 3.96E+0 4.00E+0 3.96E+0 3.96E+0 4.E+2 4.14E+0 4.12E+0 4.10E+0 4.15E+0 4.13E+0 4.11E+0 5.E+2 4.25E+0 4.23E+0 4.21E+0 4.26E+0 4.24E+0 4.22E+0 6.E+2 4.35E+0 4.32E+0 4.30E+0 4.35E+0 4.33E+0 4.31E+0 7.E+2 4.42E+0 4.40E+0 4.38E+0 4.43E+0 4.40E+0 4.39E+0
Figure 3.2: Finite conductivity, single fracture, multiple segment model verification against Cinco-Ley et al. (1978) Table 1.
3.2.3 Multiply-Fractured Horizontal Well Solution
Raghavan et al. (1997) provided an example that could be used to verify the model for a multiply-fractured horizontal well. This case requires superimposing fracture segments both in the x- and y-directions. The case inputs are shown in Table 3.5 with a pictorial representation in Figure 3.3, and the results are shown in Figure 3.4 and Table 3.6. In their paper, Raghavan et al. (1997) defined dimensionless conductivity differently from this thesis. They calculated an effective fracture conductivity that incorporates all of the fracture segments together in the following equation.
˜ Cf D = Pn j=1 q kf jwf j nkPn j=1Lf j (3.40) However, for the objectives of this thesis, using the Cinco-Ley et al. (1978)definition of individual fracture conductivity is more appropriate. In this definition, Cf D > 30 indicates
an infinite-conductivity and Cf D ≤ 30 indicates a finite-conductivity fracture. See Section
Table 3.5: Input parameters for model verification against Raghavan et al. (1997) Table 2. X Boundary, xeD 100 Y Boundary, yeD 80 Number of Fractures 5 Horizontal Length, LwD 20 Combined Fracture 100 Dimensionless Conductivity,Cef D Single Fracture 500 Dimensionless Conductivity, Cf D
Figure 3.3: Pictorial representation for the model verification of Raghavan et al. (1997) Table 2.
As in the Cinco-Ley et al.’s example (1978), Raghavan et al.’s example (1997) does not describe the number of segments that are used in the model. They state elsewhere in the paper that 40-80 fracture elements are required for each fracture depending on the Cf D.
Larger conductivities require less elements than smaller conductivities. For the comparisons provided in this thesis, the model was run with 40 fracture segments per fracture. See Section 3.3.1 for more detail on the number of fracture segments.
Figure 3.4: Finite conductivity, five multiply-fractured horizontal well, multiple segment model verification against Raghavan et al. (1997) Table 2.
3.2.4 Composite-Reservoir Solution; Verification with Numerical Model
The main motivation of this research was the absence of a relatively simple and rea-sonably accurate analytical/semi-analytical solution for a multiply-fractured horizontal well in a composite-reservoir system. Therefore, the final solution developed in this thesis was verified against a simple numerical model built in the Computer Modeling Group (CMG) software package, IMEX (CMG, 2013). Two cases were examined: a small and a large SRV. It was the original hypothesis that this method would not be correct for a large SRV because the flux distribution along the SRV boundary would depend on the contrast between the petrophysical characteristics of the SRV and the virgin reservoir as well as the size of the SRV. Specifically, if the flow at the boundary of the SRV is dominated by linear convergence,
Table 3.6: Finite conductivity, five multiply-fractured horizontal well, multiple segment model verification against Raghavan et al. (1997) Table 2.
TD
Model Output Raghavan et al. Table 2 Output
pD dpD/dLnT pD dpD/dLnT
1.0E-2 3.66E-2 1.66E-2 3.59E-2 1.65E-2 2.0E-2 5.04E-2 2.26E-2 4.94E-2 2.26E-2 4.0E-2 6.92E-2 3.04E-2 6.76E-2 3.03E-2 7.0E-2 8.83E-2 3.80E-2 8.66E-2 3.80E-2 1.0E-1 1.03E-1 4.35E-2 1.01E-1 4.34E-2 4.0E-1 1.79E-1 6.74E-2 1.78E-1 6.74E-2 1.0E+0 2.48E-1 8.19E-2 2.46E-1 8.19E-2 4.0E+0 3.85E-1 1.26E-1 3.83E-1 1.26E-1 1.0E+1 5.28E-1 1.90E-1 5.26E-1 1.90E-1 4.0E+1 8.74E-1 3.12E-1 8.72E-1 3.12E-1 1.0E+2 1.20E+0 3.90E-1 1.19E+0 3.90E-1 4.0E+2 1.80E+0 4.81E-1 1.80E+0 4.81E-1 1.0E+3 2.32E+0 7.18E-1 2.32E+0 7.18E-1 4.0E+3 4.24E+0 2.51E+0 4.23E+0 2.51E+0 1.0E+4 8.08E+0 6.55E+0 8.08E+0 6.55E+0 4.0E+4 2.87E+1 2.83E+1 2.87E+1 2.83E+1 1.0E+5 7.21E+1 7.33E+1 7.21E+1 7.33E+1
then the flux distribution would be more uniform and the accuracy of the solution would improve. This is more likely to be the case for an SRV, which does not extend much beyond the fracture tips. In contrast, a larger SRV would exhibit radial flow and more non-uniform flux distribution along the SRV boundary, which would result in the deterioration of the solution’s accuracy.
For these verification models, several simplifications were made such as constant PVT and rock properties and the stimulated region was modelled with a simple permeability contrast. Pictorial representations of the model cases are shown in Figure 3.5 and input data for the models are shown in Table 3.7 and Table 3.8 which show the constant and varied parameters across all four verification runs. Each case was run with 2 different permeability contrast (ηD) values for comparison as shown in Figure 3.6 and Figure 3.7.
(a) Small SRV case. (b) Large SRV case.
Figure 3.5: Pictorial representation of the CMG (2013) numerical model verification.
In the small SRV case (Figure 3.6), the semi-analytical model matches the numerical model output very well in both pressure and pressure-derivative; however in the large SRV case (Figure 3.7), the results of the two models diverge when the boundary is reached.
Table 3.7: CMG (2013) numerical model verification constant input parameters. Constant parameters Half length, xf 100 ft SRV permeability, kSRV 1 µD Porosity, φ 5% Oil Saturation, SO 100% Oil Viscosity, µO 0.5 cp
Formation Volume Factor, BO 1.2ST BRB
Total Compressibility, ct 3.0E-6 psi1
Sandface Flow Rate, qBH 1 RBday
Initial Reservoir Pressure, pi 5000 psi
Table 3.8: CMG (2013) Numerical model verification variable input parameters. Varied Parameters
Parameter Small SRV Case Large SRV Case
SRV x Boundary, xSRV 200 ft 400 ft
SRV y Boundary, ySRV 100 ft 400 ft
Low Permeability Contrast, ηD1 0.5 0.5
High Permeability Contrast, ηD2 0.1 0.1
Figure 3.7: Numerical model verification for 0.1 and 0.5 ηD cases in the large SRV case.
Also, the difference between the semi-analytical and numerical results is larger for the large permeability contrast (ηD = 0.1) than that for the low permeability contrast.
It should be noted that the numerical simulation of the multiply-fractured horizontal well surrounded by an SRV was a difficult endeavor. The first attempts using a multiply-fractured well and a dual porosity SRV led to very erratic results. Hence, a single fracture and a homogeneous SRV with higher effective permeability than the virgin (outer) reservoir were used in the subsequent verification runs. Even with these simplifications, the numerical model was still fairly unstable. As shown in Figure 3.6 and Figure 3.7, for example, the pressure-derivatives of the numerical model for both small and large SRV display an erratic behavior, especially after the boundary is reached. However this is a relic of the time steps used in the numerical simulations; using different time steps changed this behavior in magnitude and timing, but it could not be completely eliminated. The fact that the pressure-derivative behavior is an essential element for flow regime identifications in pressure-transient analysis and numerical solutions are not intended to provide accuracy at that scale should further support our interest in a semi-analytical solution.
3.3 Computational Considerations
There are several computation considerations that were discovered during the develop-ment of the model.
3.3.1 Number of Fracture Segments
The number of fracture segments is very important for the accuracy of the results. A single fracture segment assumes uniform flux across the entire fracture plane. Discretizing the fracture in multiple segments allows for the flux to be variable across the fracture creating a more accurate solution for both infinite and finite conductivity fractures. For early time calculations in finite conductivity fractures, there needs to be a very large number of segments for the solution to be accurate.
Figure 3.8 and Figure 3.9 exemplify this behavior for the two cases shown in Table 3.9. In the figures, the infinite conductivity case is less effected by the number of segments than the finite conductivity case. In both cases however, more fracture segments cause a decrease in the error in pressure and its derivative . Keep in mind that these plots are all log-log plots, which compress the error. For instance, in the infinite-conductivity case, at tD = 0.1, the difference between the dimensionless pressures of the 2-segment and 100-segment cases is 10%, but this error cannot be deduced from Figure 3.8.
Table 3.9: Case inputs to show different behavior with number of fracture segments. Parameter Infinite Conductivity Finite Conductivity
Fractures 10 10
Horizontal Well Length,LwD 18 18
Dimensionless Conductivity, Cf D ∞ 10
For the purposes of accuracy, larger numbers of fracture segments is better; however, there is a significant computational tradeoff to create that accuracy. Due to the nature of the matrix calculation, the computation time is approximately equivalent to the squared number of fractures. So, doubling the number of fractures will quadruple the computational
Figure 3.8: Effects of the number of fracture segments for an infinite conductivity case.
time, and ten times the fractures will increase computational time one hundred fold. When analyzing a single fracture, this is near insignificant. For example, one hundred fracture segments might only be 1 minute or 10 minutes of CPU time depending on the computer and the complexity of the system. But for a horizontal well with 20 fractures and each fracture consisting of 100 fracture-segments, the computational time would be 400 times longer (over 6.6 hours to over 66 hours). In the course of this thesis research, 10 fracture segments per fracture seemed to be appropriate to compute the pressure and pressure-derivatives with approximately 2-4% error. On a log-log plot, 2% error is hard to identify and therefore sufficient accuracy for most applications. For quantitative analysis, 40 fractures will be required to reduce error under 1% in many cases. In an example of calculating cumulative production in a rate transient model, 1% error may not be sufficient in which 100 fractures should probably be used.
3.3.2 Infinite Sums and Late Time Calculation
In the calculation of the boundary dominated responses at late times, evaluation of the infinite sums in the solution become very important for the accuracy of the results. For example, in Equation 3.28, these sums require significant computation time for convergence especially at very late times. In the computational code used for this thesis, the infinite sums were evaluated by running calculations in the background and the summations were assumed to have converged when the contribution of the additional terms became insignificant to the total summation. This was achieved by checking the convergence at specific intervals, such as every 8th summation. The convergence criteria used in this work was 10−4.
We have experienced extreme difficulties in computing the very late time behavior for the cases where the outer reservoir is infinite. This difficulty arises from the fact that the stimulated and non-stimulated inner reservoir solutions, which nearly cancel each other at late times, are many magnitudes greater than the infinite-reservoir solution. For example, the bounded inner reservoir solutions can easily be 3 orders of magnitude larger or more
