Integration Of Shale-gas-production Data And Microseismic For Fracture And Reservoir Properties With The Fast Marching Method.pdf

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

Stage:

Page: 347

Total Pages: 13

Integration of Shale-Gas-Production Data and Microseismic for Fracture and Reservoir Properties With the Fast Marching Method Jiang Xie, SPE, Changdong Yang, SPE, Neha Gupta, SPE, Michael J. King, SPE, and Akhil Datta-Gupta, SPE, Texas A&M University

Summary We present a novel approach to calculate drainage volume and well performance in shale gas reservoirs by use of the fast marching method (FMM) combined with a geometric pressure approximation. Our approach can fully account for complex fracture-network geometries associated with multistage hydraulic fractures and their impact on the well pressure and rates. The major advantages of our proposed approach are its simplicity, intuitive appeal, and computational efficiency. For example, we can compute and visualize the time evolution of the well-drainage volume for multimillion-cell geologic models in seconds without resorting to reservoir simulation. A geometric approximation of the drainage volume is then used to compute the well rates and the reservoir pressure. The speed and versatility of our proposed approach make it ideally suited for parameter estimation by means of the inverse modeling of shale-gas performance data. We use experimental design to perform the sensitivity analysis to identify the “heavy hitters” and a genetic algorithm (GA) to calibrate the relevant fracture and matrix parameters in shale-gas reservoirs by history matching of production data. In addition to the production data, microseismic information is used to help us constrain the fracture extent and orientation and to estimate the stimulated reservoir volume (SRV). The proposed approach is applied to a fractured shale-gas well. The results clearly show reduced ranges in the estimated fracture parameters and SRV, leading to improved forecasting and reserves estimation. Introduction Shale-gas production has become an important share of US energy supply, driven, to a large extent, by the advances in horizontal-well-completion technology and multistage hydraulic fracturing. To reliably estimate shale-gas reserves and ultimate recoveries, it is important to predict shale-gas well performance accounting for the relevant reservoir and fracture parameters. Currently, decline-curve analysis (Fetkovich 1980; Valko and Lee 2010) and pressure/rate-transient analysis (Ilk et al. 2010; Song and Ehlig-Economides 2011; Clarkson et al. 2012) are two types of widely used analytical methods for production forecasting in shale-gas well development. The methods in decline-curve analysis are largely curve-fitting techniques and are used to forecast production by means of extrapolation and to obtain the estimated ultimate recovery. In pressure/rate-transient analysis, fracture properties are first estimated from identified flow regimes, and then well production is predicted with the estimated properties. In addition to the commonly used analytical methods, numerical simulations were also used (Cipolla et al. 2009; Fan et al. 2010). The advantage of numerical simulation is that it can account for complex fracture geometry and heterogeneity. However, numerical simulation can be very time-consuming, particularly the high levels C 2015 Society of Petroleum Engineers Copyright V

This paper (SPE 161357) was accepted for presentation at the SPE Eastern Regional Meeting, Lexington, Kentucky, USA, 3–5 October 2012, and revised for publication. Original manuscript received for review 25 March 2013. Revised manuscript received for review 17 February 2014. Paper peer approved 22 April 2014.

of grid refinement are used to accurately model flow in the vicinity of the hydraulically fractured wells. Also, numerical models need to be calibrated with available production/pressure data to improve confidence in the prediction. The model calibration typically involves matching production/pressure history to identify fracture/ matrix parameters, either manually or through inverse modeling (Mayerhofer et al. 2006; Cipolla et al. 2009; Ghods and Zhang 2010; Yin et al. 2011). Manual history matching can be time-consuming and typically results in a single history-matched model with unknown reliability. Compared with manual or deterministic history-matching methods, stochastic search approaches such as GA and evolutionary algorithms (Schulze-Riegert et al. 2002; Yin et al. 2011) naturally provide multiple updated models that are consistent with well-performance history. Multistage hydraulic fracturing is essential to the economic recovery of shale-gas resources (Cramer 2008; Holditch 2010; King 2010). Hydraulic fracturing in shale-gas reservoirs not only creates high-conductivity transverse (primary) fractures, but also reopens pre-existing natural (secondary) fractures in the vicinity of the hydraulic fractures. This ultimately generates a complex fracture network or a stimulated region surrounding each stage of hydraulic fracturing (Fisher et al. 2005; Maxwell et al. 2002). The reservoir volume associated with the stimulated region is defined as the SRV (Mayerhofer et al. 2010). Because of the ultralow matrix permeability, much of the fluid flow in shale-gas reservoirs occurs in the interconnected fracture networks. Proper modeling of the orientation, distribution, and connectivity of the complex fracture network is critical to shale-gas-reservoir simulation and forecasting (Cipolla et al. 2011b). The common practice in the industry is to model the fracture network as evenly spaced planar fractures with a stimulated region around the fractures (Cipolla et al. 2011a). Planar-fracture models may be too simple to represent the complex-fracture network. To account for the complexfracture network, additional information about geometry and the complexity of the fracture network is required. Microseismic information, if available, can help constrain the fracture-network geometry, in particular the fracture extent and orientation (Mayerhofer et al. 2006; Cipolla et al. 2011b, 2012). During a hydraulic-fracturing job, a shear slippage occurs in the reservoir rock as the pore pressure and stress increase. This small rock failure results in the release of energy in the form of seismic waves and is known as a microseismic event (Albright and Pearson 1982; Rutledge and Philips 2003; Warpinski et al. 2004). Because of the physical nature of the microseismic event, its pattern can be used to interpret fracture geometry and azimuth and to estimate the stimulated region during reservoir characterization (Mayerhofer et al. 2010; Cipolla et al. 2012). A critical issue is how to incorporate microseismic events into the modeling of complex fracture networks. Cipolla et al. (2011b) developed an unconventional fracture model (UFM), which uses a discrete-fracture-network (DFN) model to generate natural fractures and to calibrate hydraulic fractures constrained to the treatment data and microseismic events. Also, Xu et al. (2009) and Cipolla et al. (2011b) developed a wiremesh model, which essentially is an ellipse with orthogonal fractures, calibrated to microseismic events.

April 2015 SPE Journal

347 ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

The objectives of this paper are threefold. First, we present a novel and efficient approach to calculate drainage volume and well performance in shale-gas reservoirs with an FMM combined with a geometric pressure approximation. The proposed method can be orders of magnitude faster than the conventional numerical reservoir simulator. Second, we integrate microseismic events into a history-matching workflow to calibrate fracture/matrix parameters for improved production forecast. Third, we estimate the SRV regions with the calibrated models. The integration of microseismic events in addition to rate and pressure data provides additional constraints during history matching and is thus expected to reduce parameter ranges and lead to better estimation of SRV. Methodology Our proposed approach introduces a novel forward model to predict well production given fracture and reservoir parameters. It is applied to an inverse problem to calibrate those parameters and improve production forecasts. The proposed forward model is based on the fast drainage-volume calculation (Datta-Gupta et al. 2011) and pseudosteady-state approximation that link production data to drainage volume (Xie et al. 2012), now extended to variable well rates and to highly-compressible-gas systems. Drainage-Volume Calculations With FMM. Our concept of drainage volume relies on the definition of the radius of investigation given by Lee (1982). That author defines the radius of investigation as the propagation distance of a “peak” pressure disturbance (we will call pressure “front”) for an impulse source or sink. For simplified flow geometries and homogeneous reservoir conditions, one can calculate the radius of investigation analytically. For example, we could write analytical solutions of the radius of investigation for different flow patterns as follows: r¼

pffiffiffiffiffiffi cat; . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð1Þ

where r and t are propagation distance and time of the pressure front, respectively, and hydraulic diffusivity a is defined as a ¼ k=ð/lct Þ. Moreover, c denotes different constants for different flow patterns (Kim et al. 2009). For example, for linear flow, c ¼ 2; for radial flow, c ¼ 4; and for spherical flow, c ¼ 6. However, such analytic solutions are severely limited for heterogeneous and fractured reservoirs, particularly for unconventional reservoirs with multistage hydraulic fractures. To generalize the concept to heterogeneous reservoirs, we introduce a variable called the diffusive time of flight (TOF), s. Because the pressure-front propagation has the scaling behavior of square root of time, we define the diffusive TOF as s¼

pffiffiffiffi ct: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð2Þ

Combining Eqs. 1 and 2, we obtain pffiffiffiffiffiffiffiffiffi r ¼ að~ x Þs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð3aÞ or ds 1 ¼ pffiffiffiffiffiffiffiffiffi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ð3bÞ dr að~ xÞ Consider the 1D problem with reservoir heterogeneity: Hydraulic diffusivity is now a function of location. We can integrate Eq. 3b along a trajectory to calculate diffusive TOF: sðrÞ ¼

ðr

1 pffiffiffiffiffiffiffiffiffiffi dr 0 : aðr 0 Þ 0

For 2D or 3D problems with reservoir heterogeneity, Eq. 3a can be generalized as pffiffiffiffiffiffiffiffiffi að~ x Þjrsð~ x Þj ¼ 1: . . . . . . . . . . . . . . . . . . . . . . . . ð4aÞ

Stage:

Page: 348

Total Pages: 13

For a heterogeneous reservoir, the diffusivity term a (shown in Eq. 4b as a function of spatial coordinates) takes into account the reservoir heterogeneity mainly from reservoir permeability and porosity: að~ xÞ ¼

kð~ xÞ : . . . . . . . . . . . . . . . . . . . . . . . . . . ð4bÞ /ð~ x Þlct

Eq. 4a describes the propagation of the pressure front and is called the Eikonal equation. Vasco et al. (2000) and Kulkarni et al. (2000) derived the Eikonal equation and introduced the concept of the diffusive TOF with asymptotic ray theory from geometric optics and seismology. The solution satisfies a variational principle in which the multidimensional trajectory is the one that minimizes the diffusive TOF integral. These trajectories were approximated by either ray paths (Vasco et al. 2000) or by convective streamlines (Kulkarni et al. 2000). In this paper, we do not compute the diffusive TOF by integrating along the pressure trajectories. Instead, the Eikonal equation is solved very efficiently with a class of front-tracking methods called the FMM (Sethian 1999a). The solution of the Eikonal equation provides us with the diffusive TOF, without the intermediary definition of a trajectory. We illustrate the method on a rectangular orthogonal mesh, where a finite-difference approximation is used to calculate the s gradient. We start with the discretization of the Eikonal equation on a 2D rectangular grid as shown in Eq. 5: y þy 2 2 þx maxðDx ij s; Dij s; 0Þ þ maxðDij s; Dij s; 0Þ ¼ 1=a:

                   ð5Þ Here, D is the gradient operator approximated with a first-order upwind finite-difference scheme. For example, in the x-direction, Dþx and the Dx ij s ¼ ðsiþ1; j  si;j Þ=Dx, ij s ¼ ðsij  si1; j Þ=Dx, “max” (upwinding) ensures that the front propagates monotonically outward. The same holds in the y-direction, Dy ij s ¼ ðsi; j  si; j1 Þ=Dy and Dþy ij s ¼ ðsi; jþ1  si; j Þ=Dy. The discretization results in a local quadratic equation for sij that can be solved very efficiently. The key to the FMM lies in the observation that the upwind approximation possesses a specific causality relationship. The causality relationship means the value of a point depends only on the values of its neighboring points. It made the one-pass calculation possible, which leads to the computational efficiency. Thus, we need to solve Eq. 5 concurrently in the order of increasing values of s (Sethian and Vladimirsky 2000). To illustrate the FMM, let us look at Fig. 1. We first label the well location as the “accepted” point (s ¼ 0) shown in Fig. 1a. Its adjacent nodes are labeled as “neighbor” points shown in Fig. 1b, and the remaining nodes are called “far-away” points. Now, we can iterate to calculate the diffusive TOF at each point. The detailed procedure is as follows. 1. Start from the “accepted” point in black. 2. Calculate the diffusive TOF for all “neighbor” points (A, B, C, D) with the finite-difference approximation. 3. Pick the minimum of the current “neighbor” points. a. Label it as “accepted” (Point A in Fig. 1c). b. Add its neighbors that were in “far-away” as additional “neighbor” points (Points E, F, and G in Fig. 1d). 4. Repeat Steps 2 and 3 until all the points in the domain are labeled as “accepted.” The mathematical proof of the convergence of the FMM was given by Sethian (1999b). In addition, more numerical study about the accuracy of the FMM was conducted (Datta-Gupta et al. 2011), in which the researchers demonstrated that the solution of the FMM coincided with the analytical radius of investigation in homogeneous media. After calculating the diffusive TOF with the FMM, we can convert it to physical time with Eq. 2. To calculate the well-drainage volumes, we apply different time contours and approximate the total drainage volume at cutoff as the pore volume (PV) inside the contour. We can also plot drainage volume Vp ðtÞ as a function of physical time after converting diffusive TOF to physical time.

348

April 2015 SPE Journal ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

Stage:

Page: 349

Total Pages: 13

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 1—Illustration of the FMM.

We illustrate the drainage-volume calculations with a single fracture from a multistage-fractured horizontal well in an unconventional reservoir. In this example, the reservoir is characterized with three permeability regions: the effective fracture permeability, an enhanced-permeability (102 to 103 md) region near the fracture, and the matrix permeability (104 to 105 md), as shown in Fig. 2a. The homogeneous porosity is set at 6%. From the arrival time of the pressure front calculated with the FMM, as shown in Fig. 2b, we can determine the drainage volume at different times (Fig. 2c). At approximately 100 days, the curve flattens out, which indicates that the drainage volume has reached the boundary of the enhanced-permeability zone. The drainage volume at this time is essentially the PV corresponding to the SRV. Subsequently, the drainage volume reaches the total reservoir PV, as shown by the second plateau in Fig. 2c, but only after orders of magnitude of more time. One of the major advantages of our proposed method is that the drainage volume can be calculated in seconds, even for reserlog perm (mD)

voir models with multimillion grid cells. In the following section, we describe how we can use the efficient drainage-volume calculations to estimate well-production rates under constant-flowingbottomhole-pressure (BHP) well conditions Transient-Rate Calculations With a Geometric Approximation. In this subsection, we introduce a geometric approximation with the drainage volume to compute well rates under constant well flowing BHP. The key underlying idea is to relate well-production rate and drainage volume calculated from the FMM. Starting from Darcy’s law, we basically assume pseudosteady state conditions within the drainage volume, and the gas-production rate can be derived as shown in Eq. 6. The detailed derivation is given in Appendix A. ! Tsc p mð piÞ  mð pÞ qg ðtp Þ    : . . . ð6Þ ð Vp ðtp Þ Tpsc lg Z 1 V i 1  dV k/A2 Vp ðtp Þ 0

log arrival time (day) 106

2 4 1

3

Reservoir pore volume

1 –1

0 –1

–2

–2

Drainage volume, ft3

2 0

105

SRV pore volume

104

–3

–3

–4 –4

–5 –6

(a)

(b)

103 10–1

100

101

102 Time, days

103

104

105

(c)

Fig. 2—Top view of a single fracture from a multistage-fractured horizontal well, (a) reservoir-permeability field (md) in log scale; (b) arrival time of pressure front (days) in log scale; and (c) drainage volume vs. time in log-log scale. April 2015 SPE Journal

349 ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

300 Eclipe FMM

Stage:

Page: 350

Total Pages: 13

Sensitivity analysis to identify key parameters

Production rate, Mscf/day

250 Initialize the first GA generation with LHS

200 Incorporate microseismic information

Set up simulation model

150 Simulation with FMM & Evaluate Obj. function

100

Y

50

Converge Check

Uncertainty Analysis & SRV estimation

N

0

0

50

100

150

200

250

300

350

400

Time, days Fig. 3—Comparison of production rate vs. time (Eclipse vs. FMM).

A comparison of the FMM with geometric approximation and finite-difference simulation is now presented. We adopt the same model as shown in Fig. 2, with approximately 60151100.3 million cells. Each grid in this model has the same size of 222 ft. The initial reservoir pressure is set to be 2,500 psi, and the BHP is constrained at 500 psi. One year of production is simulated to compare the gas-production rate and the computational performance with Eclipse, as shown in Fig. 3. We can see that the FMM with geometric approximation compares reasonably well with Eclipse. In this case, it only takes less than 30 seconds for the FMM with geometric approximation compared with more than 400 seconds by Eclipse, more than a tenfold increase in speedup. The benefit of computational efficiency mainly comes from the fact that the FMM involves only local quadratic solution and can be efficiently implemented with min-heap data structure (Sethian 1999b). We will use this geometric approximation as a forward model to estimate the performance of unconventional wells. With the geometric approximation, we are able to calibrate fracture and reservoir parameters, given field-production data, through inverse modeling. In this work, the GA is used for the parameter-estimation problem. Sensitivity Analysis, Experimental Design, and GA. Given a large set of potential parameters on the basis of our concepts of geologic uncertainties, we first perform a sensitivity analysis with H1

Accepted by fitness Generation updating

Fig. 4—Flow chart of GA with FMM.

a Plackett-Burman two-level initial experimental design with high and low values for each parameter (Cheng et al. 2008). Starting with a base case, only one parameter is changed to high- or lowboundary values, keeping the rest of the parameters fixed. With the geometric approximation, the pressure and rate calculations are performed for each of the experiments. Afterward, the effects of each parameter on a predefined objective function are ranked. The parameters with the strongest influence on the objective function will be kept, and other less-sensitive parameters will be discarded. This is a standard procedure for parameter screening (Schulze-Riegert et al. 2003). We adopt a stochastic-search approach to model calibration, with a modified GA (Yin et al. 2010, 2011) to adjust shale-gas reservoir and fracture parameters because the solution space is typically not smooth and GA offers great flexibility in terms of parameter selection. The GA imitates the biological principles of evolution (i.e., survival of the fittest). It has been extensively applied to the history-matching problem (Bittencourt and Horne 1997; Floris et al. 2001; Romero and Carter 2001; SchulzeRiegert et al. 2002; Williams et al. 2004). Because the GA works with a population of model parameters, rather than a single set of model parameters, it simultaneously allows for history matching as well as production forecast through multiple solutions. Experimental design is used to create the first generation of the population of parameters. Latin Hypercube Sampling (LHS) (Iman et al. 1980), with a space-filling design (Yeten et al. 2005), is used to construct an initial generation of selected key parameters. LHS is a stratified-random procedure, and provides an efficient way of sampling variables from their distributions. Unlike simple random sampling, this method ensures a full coverage of the range of each variable by maximally stratifying each marginal distribution. LHS requires fewer experiments compared with a full-factorial design or a D-optimal design. Interested readers can refer to the earlier work of Iman et al. (1980) and Yeten et al. (2005). The entire workflow is demonstrated in Fig. 4. Application In this section, we apply our approach to a horizontal well with multistage hydraulic fractures. This example is designed after a real field case. In addition, synthetic microseismic events are generated, as shown in Fig. 5. In this figure, the four different colors indicate the four fracture stages. Let us first discuss how we incorporate microseismic events into our analysis with a DFN model.

Fig. 5—Generated microseismic events.

Model Setup and Integration of Microseismic Events. Our way to model a complex fracture network is analogous to the UFM and the wire-mesh model described by Cipolla et al. (2011b). In our model, we will have a primary fracture and an

350

April 2015 SPE Journal ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

Stage:

Page: 351

Total Pages: 13

log perm (mD) –2 a

–2.5

XF

–3

b

–3.5

–4 (a)

(b)

(c)

Fig. 6—Workflow to integrate microseismic information, (a) microseismic events (red star); (b) a DFN model; and (c) generated permeability field.

associated stimulated region in each fracture stage. First, we define a primary fracture in each stage following the orientation of the microseismic cloud of that stage. Note that the fracture half-length XF is not fixed but is to be determined during parameter estimation. The stimulated region is defined as an ellipsoid with axis a, b, and c around the primary fracture. To reduce the number of history-matching parameters, we provide some extra constraints to define the ellipsoid. The axis in the z-direction, c, is fixed at twice the reservoir thickness to represent a fully penetrated fracture. The minor axis, b, in the horizontal plane can vary in a certain range whereas the major axis, a, is calculated from the fracture half-length. Further, the major and the minor axes are related as follows: a ¼ XF þ b. With such a relationship, we can ensure that the SRV extends beyond the fracture half-lengths (Song and Ehlig-Economides 2011). More importantly, our complex fracture network inside each stimulated region is represented by a heterogeneous permeability field. Here, we discard microseismic events outside the ellipsoid to account for the fracture closure and unreliable microseismic events far away from the perforation. The DFN model is used to generate a natural-fracture network and integrate microseismic events. In the DFN model, orthogonal fractures are generated around each microseismic event with fixed azimuths (following the primary fracture). The fracture properties (permeability, height, and length) are then P assigned to each simulation grid to calculate a property, I ¼ khl. The h and l in the property are the height and length of natural fracture, respectively, and k is the permeability associated with that natural fracture. This property I is a more proper way to measure the density of natural fractures within simulation grids. The effective permeability in the stimulated region is linearly distributed from the property: Minimum

corresponds to matrix permeability and KM and maximum correspond to enhanced permeability, KE (Al-Harbi et al. 2005). The result is a heterogeneous permeability field constrained to microseismic events inside the ellipsoid. Note that the enhanced permeability is the same for different fracture stages, as is the fracture permeability KF . Fig. 6 shows the workflow to integrate microseismic information. Fig. 7 shows an example of a generated permeability field. If we filter out the matrix part, Fig. 8 gives us an idea of the stimulated regions. It is worth mentioning that our SRV is actually less than the volume of the ellipsoid because part of the ellipsoid is treated as matrix (for some grid cells in the ellipsoid, the permeability equals the matrix permeability). This is clearly shown in Fig. 8. Further, we can incorporate the generated permeability field along with other properties into our forward model (FMM and geometric approximation) to obtain simulation results, as shown in Fig. 9. Part of our goal during parameter estimation is to determine the SRV through the matching of rate transients, as discussed next. Inverse Problem and Parameter Estimation. The parameters to be estimated by means of history matching and the associated ranges are listed in Table 1. We assume a uniform distribution between lower and upper bounds for each of the parameters. The base model indicates the prior knowledge of the most likely parameter values, and the reference/true model is used to generate the observed data for history matching. The horizontal well produces at a constant BHP of 1,000 psi. The first-year gas production will be integrated to predict production for the following 4 years. The geometric approximation and the FMM are used for H1

K K

1

1 0.1 0.1

0.01

0.01

0.001

0.001

0.0001

0.0001

Fig. 7—Generated permeability field (md).

Fig. 8—Visualization of the stimulated regions.

April 2015 SPE Journal

351 ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

Stage:

Arrival time (day)

H1

Page: 352

108

Drainage volume, ft3

Simulation time 100000 10000 1000 100 10 1 0.1 0.01 0.001

107

106

105

104 10–2

102 104 Time, days (b) Drainage volume versus time in log-log scale

(a) Arrival time (days) H1

104 Production rate, Mscf/day

Total Pages: 13

100

Pressure (psi)

P

103

102 0

100 150 200 250 300 350 400 Time, days (c) Production rate versus time in semi-log scale

106

3000 2900 2800 2700 2600 2500 2400 2300 2200 2100 2000

50

(d) reservoir pressure (psi) at the end of the first year

Fig. 9—Simulation results from FMM and geometric approximation, (a) arrival time; (b) drainage volume vs. time; (c) production rate vs. time; and (d) reservoir pressure at the end of first year.

the simulation of pressure and rate response. The misfit for history matching is defined as the sum of squared differences of production rate between simulation results and the reference. To evaluate the impact of various parameters, a sensitivity analysis is performed with 1 year of simulation, with the rapidperformance prediction techniques. Fig. 10 shows a tornado plot of the objective function (logarithm of rate misfit). The line in the middle shows the objective function of the base model. We vary one parameter at a time and rerun the simulation. Low (dark green) and high (yellow) represent the impact of low and high values of each parameter in Table 1. From Fig. 10, we can easily tell that the fracture permeability and enhanced permeability have the most significant influences. It is also worth pointing out the low

impact of the matrix permeability. This is a result of the short simulation time period, during which the pressure disturbance has not yet reached the matrix. The GA is carried out for 15 generations, and each generation has 80 realizations. Fig. 11 shows the gas-production rate in both the semilog plot and log-log plot for the first generation. The reference model is shown in red dots, whereas green lines denote initial models. We can see the wide range of predictions from the initial models indicating the initial parameter ranges. After 15 generations, all simulation misfits vs. generation number are plotted in Fig. 12. We select the best 50 models with a misfit cutoff and plot the gas-production rate for the selected models, as shown in Fig. 13.

TABLE 1—HISTORY-MATCHING PARAMETER RANGES AND THE REFERENCE Model Parameters

Base

Lower

Upper

Reference/True

Matrix permeability (KM) Enhanced permeability (KE) Fracture permeability (KF) Fracture-1 half-length (XF1) Fracture-2 half-length (XF2) Fracture-3 half-length (XF3) Fracture-4 half-length (XF4) Ellipsoid-1 minor axis (b1) Ellipsoid-2 minor axis (b2) Ellipsoid-3 minor axis (b3) Ellipsoid-4 minor axis (b4)

12105md 0.022 md 2.4 md 390 ft 600 ft 600 ft 320 ft 180 ft 260 ft 230 ft 190 ft

8105 md 0.005 md 0.5 md 300 ft 400 ft 350 ft 250 ft 150 ft 200 ft 200 ft 170 ft

15105 md 0.05 md 4.0 md 530 ft 800 ft 700 ft 450 ft 250 ft 300 ft 300 ft 230 ft

10105 md 0.02 md 2.50 md 400 ft 650 ft 550 ft 300 ft 200 ft 250 ft 250 ft 200 ft

352

April 2015 SPE Journal ID: balamurali.l Time: 15:05 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

KF KE XF2 XF3 XF4 XF1 B3 B4 B2 B1 Low High

KM 6.0

6.2

6.4

6.6

6.8 7.0 7.2 7.4 Objective function

7.6

7.8

8.0

Fig. 10—Sensitivity analysis of history-matching parameters.

To further assess the history-matching results, an additional 4-year simulation is carried out to predict the performance (blue curves in Fig. 14). They all give a gas-rate interval of 0.1 to 0.2 MMscf/D at the end of the fifth year. This wide range is dominantly controlled by the spread of the matrix permeability, which is not constrained during history matching because the flow in the first year mainly is inside the stimulated region. As the pressure front reaches the boundary of the stimulated region, matrix permeability plays a more and more important role in predicting the production rate. The spread in the model prediction can also be explained by comparing the initial and the updated fracture and matrix permeability after history matching. These results are summarized in the boxplot of permeability in Fig. 15. Here, three permeability parameters were normalized to fall between zero and unity. The range of model parameters in the population is indicated by the blue box, with reference indicated by triangle and median in red line. After the history matching of the first-year production data, the range of the matrix permeability is not reduced whereas the fracture and enhanced permeability converge to the reference with very narrow ranges. We acknowledge that full uncertainty quantification under a Bayesian framework could be performed to assess the parameter uncertainties. For example, we can use Markov-chain Monte

Gas Rate (Mscf/Day)

obs

initial models

initial models

Gas Rate (Mscf/Day) 150 200 250 300 Simulation Time (Days) (a)

Initial Models

obs

103

100

Total Pages: 13

SRV Estimation. From our drainage-volume calculations, we can estimate the stimulated reservoir PV. As shown in Fig. 2c, the first plateau in the drainage-volume plot indicates the PV corresponding to the SRV. For this application, we plot drainage volume vs. time in log-log scale for the reference and the fifty updated models in Fig. 18. From this figure, we can see the curves bending, instead of forming a plateau, because of the permeability heterogeneity in the stimulated region. The inflection point in Fig. 18 in which we see a reversal in concavity is an indication of the SRV PV. The two dark green lines give us an estimate of the range of the SRV PV. The CDFs of the SRV PVs of the initial and updated models along with that of the reference model are plotted in Fig. 19. The SRV PV of the reference model is 5.534 MMcf. After history matching, the range of the updated SRV PV is considerably reduced compared with that of the initial model, as shown in Fig. 19. The mean and the standard deviation of the SRV PV from the updated models are 5.559 and 0.367 MMcf, respectively. In Fig. 20, the SRV of the reference model is compared with the SRVs of a set of randomly selected models from the updates. We can see that the SRV of each fracture stage can vary substantially, whereas the total SRV PV is constrained, as shown in Fig. 20. The reduced range in the estimated fracture parameters and the

105

104

50

Page: 353

Carlo to sample the posterior distribution on the basis of given prior distribution. However, the limitation is that the prior distribution is unknown to us. So, the topic is beyond the scope of this paper. During history matching, we also estimate the fracture halflengths in each stage. Instead of examining the fracture half-lengths in each individual stage, which are likely to vary considerably, we compare the total-fracture half-lengths from all stages and plot the cumulative distribution functions (CDFs) from the initial and the updated models, as shown in Fig. 16. The total-fracture halflength of the reference model is 1,900 ft, as shown by the red vertical line. From these results, we can clearly see the reduced range in updated models as indicated by the reduced range and spread of the updated distribution compared with the initial distribution of total-fracture half-lengths. In the fifty updated models shown in Fig. 16, the mean and standard deviation of total-fracture halflength are 1,836 and 116 ft, respectively. This compares very favorably to the reference model. In Fig. 17, we have also shown the cumulative distribution of the composite quantity, X pffiffiffiffiffiffi XF KF , which is closely related to the flow across the fracture surfaces. The reduction in range is even more significant here.

Initial Models

105

102 0

Stage:

350

400

104

103

102 10–1

100

101 102 Simulation Time (Days)

103

(b)

Fig. 11—The simulation results with initial 80 models (green) compared with reference model (red), (a) in semilog scale and (b) in log-log scale. April 2015 SPE Journal

353 ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

7.5 7 6.5 Misfit

6 5.5 5 4.5 4 3.5 0

2

4

6 8 10 Generation Number

14

12

16

Fig. 12—The objective function vs. generation number.

Total Pages: 13

Conclusions • We have presented a novel approach to calculate drainage volume and well performance in shale-gas reservoirs with an FMM combined with a geometric pressure approximation. The current work extends our previous work on drainage-volume and pressure calculations that were based on constant well rates to constant well-flowing-pressure conditions. Our approach can fully account for complex-fracture-network geometries associated with multistage hydraulic fractures and their impact on the well pressure and rates. The major advantages of our proposed approach are its simplicity, intuitive appeal, and computational efficiency. The proposed method can be orders of magnitude faster than the conventional numerical reservoir simulator. • The speed and versatility of our proposed approach make it ideally suited for parameter estimation by means of the inverse modeling of shale-gas performance data. We use experimental design to perform sensitivity analysis to identify the “heavy hitters” and a GA to calibrate the relevant fracture and matrix parameters in shalegas reservoirs by the history matching of production data.

Updated Models

105

Page: 354

SRV leads to improved production forecasting and reserves estimation.

GA Evolution 8

3

Stage:

Updated Models

105

obs history matching

Gas Rate (Mscf/Day)

Gas Rate (Mscf/Day)

obs history matching

104

103

102 0

50

100

150

200

250

300

350

400

104

103

102 10–1

Simulation Time (Days)

100 101 102 Simulation Time (Days)

103

Fig. 13—Selected fifty models (green) compared with reference model (red), (a) in semilog scale and (b) in log-log scale.

Prediction with Updated Models

Prediction with Updated Models

obs

obs

history matching prediction

history matching prediction

104

103

102 0

105

Gas Rate (Mscf/Day)

Gas Rate (Mscf/Day)

105

200 400 600 800 1000 1200 1400 1600 1800 2000 Simulation Time (Days) (a)

104

103

102 10–1

100

101 102 Simulation Time (Days) (b)

103

104

Fig. 14—Four-year prediction with updated models, (a) in semilog scale and (b) in log-log scale. 354

April 2015 SPE Journal ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

Box Plot of Permeabilities (updated)

1

1

0.8

0.8

0.6

0.4

0.2 0

KM

KE

KF

1

2 (a)

3

0.6

0.4

0.2 0

Stage:

KM

KE

KF

1

2 (b)

3

Page: 355

Total Pages: 13

1 updated

0.9 Cumulative distribution function

Box Plot of Permeabilities (initial)

Normalized distribution

Normalized distribution

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

initial reference

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 Total fracture half-length, ft Fig. 16—CDF plot of total-fracture half-length.

Fig. 15—Permeability (fracture, enhanced, and matrix) ranges: (a) initial range and (b) updated range.

1

108

updated

reference updated

initial reference

0.8

107

0.7 Drainage volume, ft3

Cumulative distribution function

0.9

0.6 0.5 0.4 0.3 0.2

106 Range of SRV pore volume

105

0.1 0 1000

1500

2000

2500

∑ XF

3000

3500

4000

4500

√KF , ft . mD1/2

Cumulative distribution function

102

104

106

Fig. 18—Drainage volume vs. time in log-log scale.

• The power and utility of our approach was illustrated with the history matching of performance data of a shale-gas well. We have integrated microseismic events into a history-matching workflow to calibrate fracture/matrix parameters for improved production forecast and estimate the SRV regions with the calibrated models. The integration of microseismic events in addition to rate-transient data considerably reduced the parameter ranges and led to better estimation of the SRV.

1 updated initial reference

0.8

100

Time, days

pffiffiffiffiffiffiffi Fig. 17—CDF plot of XF KF .

0.9

104 10–2

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

4

4.5

5

5.5

6

6.5

SRV, ft3 Fig. 19—CDF plot of SRV PV.

7

7.5 × 106

Nomenclature a ¼ major axis in horizontal ellipse, ft A ¼ cross-sectional area, ft2 b ¼ minor axis in horizontal ellipse, ft B ¼ FVF, RB/bbl or res cf/scf c ¼ axis in z-direction, ft ct ¼ total compressibility, psi1 D ¼ spatial-derivative operator h ¼ reservoir thickness, ft i ¼ index in x-direction j ¼ index in y-direction k ¼ permeability, md mðpÞ ¼ pseudopressure, psi

April 2015 SPE Journal

355 ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

H1

Stage:

H1

Page: 356

Total Pages: 13

H1

Simulation time

Simulation time

Simulation time

100 10 1

100 10 1

100 10 1

0.1 0.01 0.001

0.1 0.01 0.001

0.1 0.01 0.001

the reference model H1 Arrival time (days) 100 10 1 0.1

Simulation time

Simulation time

100 10 1

100 10 1

0.1 0.01 0.001

0.1 0.01 0.001

0.01 0.001 Fig. 20—The SRV comparison, the reference vs. four updated models.

p ¼ pressure, psi p ¼ average pressure, psi pwf ¼ flowing BHP, psi q ¼ Darcy flux, RB/D qD ¼ dimensionless rate qwell ¼ well-production rate, RB/D, res cf/D, B/D, scf/D r ¼ distance, ft t ¼ time, days T ¼ temperature,  F V ¼ drainage volume, integral variable, ft3 Vp ðtÞ ¼ drainage volume, ft3 Z ¼ the compressibility factor a ¼ hydraulic diffusivity, ft2/D / ¼ porosity pffiffiffiffiffiffiffi s ¼ diffusive TOF, day l ¼ viscosity, cp Subscripts g ¼ gas i ¼ initial reservoir condition o ¼ oil sc ¼ standard condition Acknowledgments The authors would like to acknowledge financial support from members of the Texas A&M University Joint Industry Project, MCERI (Model Calibration and Efficient Reservoir Imaging). References Albright, J.N. and Pearson, C.F. 1982. Acoustic Emissions as a Tool for Hydraulic Fracture Location: Experience at the Fenton Hill Hot Dry Rock Site. SPE J. 22 (4): 523–530. SPE-9509-PA. http://dx.doi.org/ 10.2118/9509-PA. Al-Harbi, M., Cheng, H., He, Z. et al. 2005. Streamline-Based Production Data Integration in Naturally Fractured Reservoirs. SPE J. 10 (4): 426–439. Bittencourt, A.C. and Horne, R.N. 1997. Reservoir Development and Design Optimization. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 5–8 October. SPE-38895MS. http://dx.doi.org/10.2118/38895-MS. Cheng, H., Dehghani, K., and Billiter, T. 2008. A Structured Approach for Probabilistic-Assisted History Matching Using Evolutionary Algorithms: Tengiz Field Applications. Presented at the SPE Annual Tech-

nical Conference and Exhibition, Denver, Colorado, 21–24 September. SPE-116212-MS. http://dx.doi.org/10.2118/116212-MS. Cipolla, C.L., Lolon, E.P., Erdle, J.C. et al. 2009. Modeling Well Performance in Shale-Gas Reservoirs. Presented at the SPE/EAGE Reservoir Characterization and Simulation Conference, Abu Dhabi, UAE, 19–21 October. SPE-125532-MS. http://dx.doi.org/10.2118/125532-MS. Cipolla, C.L., Fitzpatrick, T., Williams, M.J. et al. 2011a. Seismic-to-Simulation for Unconventional Reservoir Development. Presented at the SPE Reservoir Characterization and Simulation Conference, Abu Dhabi, UAE, 9–11 October. SPE-146876-MS. http://dx.doi.org/ 10.2118/146876-MS. Cipolla, C.L., Weng, X., Mack, M.G. et al. 2011b. Integrating Microseismic Mapping and Complex Fracture Modeling to Characterize Hydraulic Fracture Complexity. Presented at the SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, 24–26 January. SPE-140185-MS. http://dx.doi.org/10.2118/140185-MS. Cipolla, C.L., Maxwell, S., Mack, M.G. et al. 2012. Engineering Guide to the Application of Microseismic Interpretations. Presented at the SPE Hydraulic Fracturing Technology Conference, The Woodlands, Texas, 6–8 February. SPE-152165-MS. http://dx.doi.org/10.2118/152165MS. Clarkson, C.R., Nobakht, M., Kaviani, D. et al. 2012. Production Analysis of Tight-Gas and Shale-Gas Reservoirs Using the Dynamic-Slippage Concept. SPE J. 17 (1): 230–242. SPE-144317-PA. http://dx.doi.org/ 10.2118.144317-PA. Cramer, D.D. 2008. Stimulating Unconventional Reservoirs: Lesson Learned, Successful Practices, Area for Improvement. Presented at the SPE Unconventional Reservoirs Conference, Keystone, Colorado, 10–12 February. SPE-114172-MS. http://dx.doi.org/10.2118/114172MS. Datta-Gupta, A., Xie, J., Gupta, N. et al. 2011. Radius of Investigation and Its Generalization to Unconventional Reservoirs. J Pet Technol 63 (7): 52–55. Fan, L., Thompson, J.W., and Robinson, J.R. 2010. Understanding Gas Production Mechanism and Effectiveness of Well Stimulation in the Haynesville Shale Through Reservoir Simulation. Presented at the Canadian Unconventional Resources and International Petroleum Conference, Calgary, Alberta, Canada, 19–21 October. SPE-136696-MS. http://dx.doi.org/10.2118/136696-MS. Fetkovich, M.J. 1980. Decline Curve Analysis Using Type Curves. J Pet Technol 32 (6): 1065–1077. Fisher, M.K., Wright, C.A., Davidson, B.M. et al. 2005. Integrating Fracture-Mapping Technologies to Improve Stimulations in the Barnett Shale. SPE Prod & Fac 20 (2): 85–93. SPE-77441-PA. http://dx.doi. org/10.2118/77441-PA.

356

April 2015 SPE Journal ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

Floris, F.J.T., Bush, M.D., Cuypers, M. et al. 2001. Methods for Quantifying the Uncertainty of Production Forecast: A Comparative Study. Petroleum Geosci. 7 (S): S87–S96. Ghods, P. and Zhang, D. 2010. Ensemble-Based Characterization and History Matching of Naturally Fractured Tight/Shale Gas Reservoirs. Presented at the SPE Western Regional Meeting, Anaheim, California, 27–29 May. SPE-133606-MS. http://dx.doi.org/10.2118/133606-MS. Holditch, S.A. 2010. Shale Gas Holds Global Opportunities. The American Oil & Gas Reporter, August 2010 Editor’s Choice. Ilk, D., Stotts, G.W.J., Anderson, D.M. et al. 2010. Production Data Analysis—Challenges, Pitfalls, Diagnostics. SPE Res Eval & Eng 13 (3): 538–552. SPE-102048-PA. http://dx.doi.org/10.2118/102048-PA. Iman, R.L., Davenport, J.M., and Zeigler, D.K. 1980. Latin Hypercube Sampling (Program User’s Guide). [LHC, in FORTRAN] Kim, J.U., Datta-Gupta, A., Brouwer, R. et al. 2009. Calibration of HighResolution Reservoir Models Using Transient Pressure Data. Presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana, 4–7 October. SPE-124834-MS. http//dx.doi.org/ 10.2118/124834-MS. King, G.E. 2010. Thirty Years of Gas Shale Fracturing: What Have We Learned? Presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19–22 September. SPE-133456-MS. http:// dx.doi.org/10.2118/133456-MS. Kulkarni, K.N., Datta-Gupta, A., and Vasco, D.W. 2000. A Streamline Approach for Integrating Transient Pressure Data Into High Resolution Reservoir Models. Presented at the SPE European Petroleum Conference, Paris, France, 24–25 October. SPE-65120-MS. http://dx.doi.org/ 10.2118/65120-MS. Lee, W.J. 1982. Well Testing. SPE Textbook Series. Richardson, Texas: Society of Petroleum Engineers. Lee, W.J., Rollins, J.B., and Spivey, J.P. 2003. Pressure-Transient Testing. SPE Textbook Series. Richardson, Texas: Society of Petroleum Engineers. Maxwell, S.C., Urbancic, T.I., Steinsberger, N. et al. 2002. Microseismic Imaging of Hydraulic Fracture Complexity in the Barnett Shale. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 29 September–2 October. SPE-77440-MS. http://dx. doi.org/10.2118/77440-MS. Mayerhofer, M.J., Lolon, E.P., Warpinski, N.R. et al. 2010. What Is Stimulated Reservoir Volume? SPE Prod & Oper 25 (1): 89–98. SPE119890-PA. http://dx.doi.org/10.2118/119890-PA. Mayerhofer, M.J., Lolon, E.P., Youngblood, J.E. et al. 2006. Integration of Microseismic-Fracture-Mapping Results With Numerical Fracture Network Production Modeling in Barnett Shale. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 24–27 September. SPE-102103-MS. http://dx.doi.org/10.2118/102103-MS. Nordbotten, J.M., Celia, M.A., and Bachu, S. 2004. Analytical Solutions for Leakage Rates Through Abandoned Wells. Water Resources Res. 40: W04204. Romero, C.E. and Carter, J.N. 2001. Using Genetic Algorithm for Reservoir Characterization. J. Petrol. Sci. & Eng. 31 (2–4): 113–123. Rutledge, J.T. and Philips, W.S. 2003. Hydraulic Stimulation of Natural Fractures as Revealed by Induced Microearthquakes, Carthage Cotton Valley Gas Field, East Texas. Geophysics 68 (2): 441. Schulze-Riegert, R.W., Axmann, J.K., Haase, O. et al. 2002. Evolutionary Algorithms Applied to History Matching of Complex Reservoirs. SPE Res Eval & Eng 5 (2): 163–173. SPE-77301-PA. http://dx.doi.org/ 10.2118/77301-PA. Schulze-Riegert, R.W., Haase, O., and Nekrassov, A. 2003. Combined Global and Local Optimization Techniques Applied to History Matching. Presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 3–5 February. SPE-79668-MS. http://dx.doi.org/10.2118/79668MS. Sethian, J.A. 1999a. Fast Marching Method. SIAM Rev. 41 (2): 199–235. Sethian, J.A. 1999b. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Material Science. Series: Cambridge Monographs on Applied and Computational Mathematics; 3. Cambridge, UK: Cambridge University Press. Sethian, J.A. and Vladimirsky, A. 2000. Fast Methods for the Eikonal and Related Hamilton-Jacobi Equations on Unstructured Meshes. Proc. Natl. Acad. Sci. USA 97 (11): 5699–5703.

Stage:

Page: 357

Total Pages: 13

Song, B. and Ehlig-Economides, C.A. 2011. Rate-Normalized Pressure Analysis for Determination of Shale Gas Well Performance. Presented at the North American Unconventional Gas Conference and Exhibition, The Woodlands, Texas, 14–16 June. SPE-144031-MS. http:// dx.doi.org/10.2118/144031-MS. Valko, P.P. and Lee, W.J. 2010. A Better Way To Forecast Production From Unconventional Gas Wells. Presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19–22 September. SPE-134231-MS. http://dx.doi.org/10.2118/134231-MS. Vasco, D.W., Keers, H., and Karasaki, K. 2000. Estimation of Reservoir Properties Using Transient Pressure Data: An Asymptotic Approach. Water Resources Res. 36 (12): 3447–3465. Warpinski, N.R., Wolhart, S.L., and Wright, C.A. 2004. Analysis and Prediction of Microseismicity Induced by Hydraulic Fracturing. SPE J. 9 (1): 24–33. SPE-87673-PA. http://dx.doi.org/10.2118/87673-PA. Williams, G.J.J., Mansfield, M., MacDonald, D.G. et al. 2004. Top-Down Reservoir Modelling. Presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, 26–29 September. SPE-89974MS. http://dx.doi.org/10.2118/89974-MS. Winestock, A.G. and Colpitts, G.P. 1965. Advances in Estimating Gas Well Deliverability. J Can Pet Technol 4 (3): 111–119 Xie, J., Gupta, N., King, M.J. et al. 2012. Depth of Investigation and Depletion Behavior in Unconventional Reservoirs Using Fast Marching Methods. Presented at the EAGE Annual Conference and Exhibition Incorporating SPE Europec, Copenhagen, Denmark, 4–7 June. Xu, W., Thiercelin, M., and Walton, I. 2009. Characterization of Hydraulically Induced Shale Fracture Network Using a Semi-Analytical Model. Presented at the Annual Technical Conference and Exhibition, New Orleans, Louisiana, 4–7 October. SPE-124697-MS. http://dx.doi.org/ 10.2118/124697-MS. Yeten, B., Castellini, A., Guyaguler, B. et al. 2005. A Comparison Study on Experimental Design and Response Surface Methodologies. Presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, 31 January–2 February. SPE-93347-MS. http://dx.doi.org/10. 2118/93347-MS. Yin, J., Park, H., Datta-Gupta, A. et al. 2010. A Hierarchical StreamlineAssisted History Matching Approach With Global and Local Parameter Updates. Presented at the SPE Western Regional Meeting, Anaheim, California, 27–29 May. SPE-132642-MS. http://dx.doi.org/10. 2118/132642-MS. Yin, J., Xie, J., Datta-Gupta, A. et al. 2011. Improved Characterization and Performance Assessment of Shale Gas Wells by Integrating Stimulated Reservoir Volume and Production Data. Presented at the SPE Eastern Region Meeting, Columbus, Ohio, 17–19 August. SPE148969-MS. http://dx.doi.org/10.2118/148969-MS.

Appendix A We start with Darcy’s law in radial flow, qðr; tÞ ¼

kAðrÞ @pðr; tÞ : . . . . . . . . . . . . . . . . . . . . ðA-1Þ l @r

We can express Darcy’s law with the reservoir drainage volume as the spatial coordinate, dVp ¼ /AðrÞ  dr: . . . . . . . . . . . . . . . . . . . . . . . . ðA-2Þ From chain rule, we have @ dVp @ @ ¼ ¼ /AðrÞ : . . . . . . . . . . . . . . . . ðA-3Þ  @r @Vp dr @Vp Substituting Eq. A-3 into Eq. A-1 results in qðVp ; tÞ ¼

k/A2 @pðVp ; tÞ : . . . . . . . . . . . . . . . . . . . ðA-4Þ l @Vp

In Eq. A-4, the pressure and Darcy flux are expressed in terms of the drainage volume. Permeability is the averaged permeability over the surface area A along the boundary of the drainage volume and so is the porosity. Notice that in this formulation, the drainage volume is considered as a new spatial variable. We could have equivalently chosen to treat the diffusive TOF as the spatial variable.

April 2015 SPE Journal

357 ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

The Darcy flux in Eq. A-4 can be approximated as the production rate at the well location multiplied by a dimensionless flux along the drainage volume as follows (Winestock and Colpitts 1965): qðVp ; tÞ  qwell ðtÞ  qD ðVp ; tÞ: . . . . . . . . . . . . . . . .

ðA-5Þ

Combining Eqs. A-4 and A-5, we obtain qwell ðtÞ  qD ðVp ; tÞ 

k/A2 @pðVp ; tÞ : . . . . . . . . . . . . ðA-6Þ l @Vp

Rearranging Eq. A-6 and integrating result in ð Vp ðtÞ qD ðV; tÞ dV: . . . . . . . ðA-7Þ Dp ¼ pi  pwf  qwell ðtÞl k/A2 0 In principle, the integral is taken to the outer boundary of the system. However, in practice, the dimensionless flux is sufficiently small beyond the depth of investigation that this integral can be restricted over this finite moving boundary. Because the well is being produced under a constant flowing BHP drop, Dp is considered to be specified, leading to the calculation of the well rate. We now make a few assumptions and observations to quantify the dimensionless flux qD in Eq. A-7 (Nordbotten et al. 2004): • Pseudosteady-state conditions exist inside the drainage volume. We can rewrite the diffusivity equation in terms of drainage volume in the following form: ct

@p @q @qD ðVp Þ ¼ ¼ qwell : . . . . . . . . . . . . . . . . . ðA-8Þ @t @Vp @Vp

Under pseudosteady-state conditions, the left-hand side of Eq. A8 remains a constant. Thus, the dimensionless flux qD must be a linear function of drainage volume: qD ðVÞ ¼ a þ b  V: • Inner-boundary condition: Darcy flux is equal to production rate at the well, qðVp ; tÞ ¼ qwell ðtÞ, which yields qD ð0Þ ¼ 1; thus, a ¼ 1. • Outer-boundary condition: Darcy flux is equal to zero at the depth of investigation, qðVp ; tÞ ¼ 0, which yields qD ðVp Þ ¼ 0; thus, b ¼ 1=Vp ðtÞ. With these assumptions, we can formulate the dimensionless flux as V : . . . . . . . . . . . . . . . . . . . . . . . .ðA-9Þ qD ðVÞ ¼ 1  Vp ðtÞ Now, we substitute Eq. A-9 into Eq. A-7 and rearrange the equation to obtain an expression for the well rate: qwell ðtÞ 

Dp ð l Vp ðtÞ 0

1   : . . . . . . . . .ðA-10Þ 1 V dV 1  k/A2 Vp ðtÞ

In Eq. A-10, the drainage volume Vp ðtÞ was calculated from the fast-marching calculation. The surface area A is the outer surface of the drainage volume and can also be estimated during the drainage-volume calculations. Note that the permeability and porosity are averaged properties along the drainage volume. The production rate is in reservoir conditions, either RB/D or res cf/D. The oil- or gas-formation volume factor (FVF) is needed to convert to standard conditions: Dp qwell ðtÞ  ð Bl Vp ðtÞ 0

1   : . . . . . . . . .ðA-11Þ 1 V 1 dV k/A2 Vp ðtÞ

We can also calculate the pressure at a location within the reservoir or at the wellbore. Combining Eqs. A-8 and A-9 results in ct

@p qwell ðtÞ  : ...................... @t Vp

ðA-12Þ

Stage:

Page: 358

Total Pages: 13

Integrating, we obtain the pressure: ðt qwell ðtÞ pðVp ; tÞ  pi  dt: . . . . . . . . . . . . . . . . ðA-13Þ tðVp Þ ct Vp The integral begins at the time the pressure front reaches that location. The average pressure inside the drainage volume is given as ð 1 Vp ðtÞ p¼ pðVp ; tÞdVp : . . . . . . . . . . . . . . . . . ðA-14Þ Vp ðtÞ 0 In gas reservoirs, the gas viscosity and compressibility can be strongly dependent on the reservoir pressure. To account for these effects, pseudopressure and pseudotime are commonly used in well-test analysis (Lee et al. 2003). Here, we adopt the normalized pseudopressure and pseudotime: ðt dt . . . . . . . . . . . . . . . . . . . . . . ðA-15Þ tp ¼ ðlg ct Þi l 0 g ct and mðpÞ ¼

  ðp lg Z pdp : . . . . . . . . . . . . . . . . . . . . ðA-16Þ p i 0 lg Z

After substituting Eqs. A-15 and A-16, the diffusivity equation takes on the following form (Lee 1982): @mðpÞ k ¼ r2 mðpÞ: . . . . . . . . . . . . . . . . . . ðA-17Þ @tp /ðlg ct Þi Introducing the diffusive TOF, the Eikonal equation can now be rewritten as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k jrsj ¼ 1: . . . . . . . . . . . . . . . . . . . . . . ðA-18Þ /ðlg ct Þi Instead of tracking the pressure front with physical time, we track the pseudopressure front propagation with pseudotime. The diffusive TOF s corresponds to the arrival time of the pseudopressure front. Notice that the gas viscosity and total compressibility are taken at initial reservoir conditions. With a similar derivation as before, we can write the gas-production rate as ! Tsc p mð pi Þ  mð pÞ qg ðtp Þ    ; ð Vp ðtp Þ Tpsc lg Z 1 V i dV 1  k/A2 Vp ðtp Þ 0                    ðA-19Þ where T is the reservoir temperature. Jiang Xie is research scientist at Chevron Energy Technology Company. Xie holds a BS degree in polymer science and engineering from University of Science and Technology of China and MS and PhD degrees in petroleum engineering from Texas A&M University. Changdong Yang is currently a PhD candidate in petroleum engineering at Texas A&M University. Yang holds a BS degree in mechanics from Peking University, China, and an MS degree in petroleum engineering from Texas A&M University. Neha Gupta is currently a reservoir engineer at Occidental Oil and Gas Corporation, Houston. She holds an MS degree in petroleum engineering from Texas A&M University and a BS degree in petroleum engineering from Indian School of Mines. Gupta has almost 4 years of industry experience working as reservoir engineer with Shell and Oxy. She has expertise in the areas of reserves estimation, reservoir simulation, and reservoir management. Michael J. King is a professor and holder of the LeSuer Chair in Reservoir Management in the Petroleum Engineering Department at Texas A&M University. Previously, he was a senior adviser

358

April 2015 SPE Journal ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042

J161357 DOI: 10.2118/161357-PA Date: 10-April-15

in reservoir modeling and simulation at BP America, where he worked in a variety of reservoir-engineering leadership, research, and operational roles beginning in 1982. King holds both MS and PhD degrees in theoretical physics from Syracuse University and a BS degree in mathematics and physics from Cooper Union. Akhil Datta-Gupta is Regents Professor and holder of the L.F. Peterson ’36 Endowed Chair in Petroleum Engineering at Texas A&M University. Previously, he worked for BP Exploration and

Stage:

Page: 359

Total Pages: 13

Research and for the Lawrence Berkeley National Laboratory. Datta-Gupta holds a PhD degree from the University of Texas at Austin. He is the recipient of the 2009 SPE John Franklin Carl Award and the 2003 SPE Lester C. Uren Award. Datta-Gupta also received the SPE Cedric K. Ferguson Certificate twice (2000 and 2006) and the 1992 AIME Rossiter W. Raymond Memorial Award. He is an SPE Distinguished Member and a member of the US National Academy of Engineering.

April 2015 SPE Journal

359 ID: balamurali.l Time: 15:06 I Path: S:/3B2/J###/Vol00000/140042/APPFile/SA-J###140042