INTEGRATED METHOD TO EVALUATE ACID STIMULATION OF

HORIZONTAL WELLS IN CARBONATE RESERVOIR THROUGH TREATMENT

PRESSURE ANALYSIS

A Thesis

by

KENJI UEDA

Submitted to the Office of Graduate and Professional Studies of

Texas A&M University
in partial fulfillment of the requirements for the degree of

MASTER OF SCIENCE

Chair of Committee, Ding Zhu

Committee Members, A. Daniel Hill
Marcelo Sanchez

Head of Department, A. Daniel Hill

August 2015

Major Subject: Petroleum Engineering

Copyright 2015 Kenji Ueda

 ABSTRACT

Unlocking a tight carbonate formation for oil and gas production by multi-stage acid

stimulation is a relatively cost-effective method as an alternative to propped fracturing for

production enhancement. Depending on whether treatment pressure is below or above the

formation closure stress, acid stimulation is basically divided into matrix acidizing and

acid fracturing. In this study, practical methodology to evaluate both matrix acidizing and

acid fracturing through treatment monitoring is presented respectively.

For matrix acidizing, monitoring and optimizing a matrix acidizing has been

achieved by integrating a forward model used in acidizing design for horizontal wells with

a real-time monitoring model for skin evolution during the stimulation. The effect of

acidizing is described as an overall skin factor change, and productivity improvement is

predicted for the treatment. Then the field treatment data monitored on-site was used to

estimate the skin response by treatment injection. History matching procedure of design

and actual treatment data will be carried out to update near-wellbore and key wormholing

parameters. Through sensitivity study, which parameter should be updated is discussed.

Finally optimum rate schedule is identified based on updated parameters.

Meanwhile, for acid fracturing treatment, new method for real-time monitoring of

acid fracturing, the inverse injectivity vs. superposition time function plot is proposed,

subject to the condition that the treatment pressure is above closure pressure after the

breakdown. Combining a linear dual-porosity transient slab model with injectivity concept,

actual growing cross-sectional area induced by acid fracturing treatment can be monitored

ii
in real-time. After production starts, linear flow diagnostic approach with rate-transient

analysis provides cross-sectional area flowing from matrix, which is compared with the

area induced by acid fracturing during the stimulation. The treatment efficiency provides

engineers with additional information as to whether the designed acid fracturing was

performed appropriately under the in-situ closure stress field.

A field case example of both multi-stage matrix acidizing and acid fracturing acid

in horizontal well are also presented respectively in the study to illustrate the application

of the approach developed, and to show the value of the integrated approach to monitor

and diagnose acid stimulation in horizontal wells.

iii
 DEDICATION

To my wife, Risako, little one, Riku and unborn baby.

iv
 ACKNOWLEDGEMENTS

I would like to thank my committee chair, Dr. Zhu, and my committee members, Dr. Hill,

and Dr. Sanchez, for their guidance and support throughout the course of this research.

Thanks also go to my colleagues in room 714 for making my time at Texas A&M

University a great experience. I also want to extend my gratitude to Nozomu Yoshida for

his help with friendship and brain-storming through a coffee break.

Finally, I thank my wife, Risako, for her patience, encouragement and love, and

my little one, Riku, for making our life fun in College Station. I could not have completed

without their support.

v
 NOMENCLATURE

A pipe cross-sectional area, L2, ft2

Ac cross-sectional area to flow, ft2

Acm total matrix surface area draining into fracture system, ft2

B formation volume factor, L3/L3, RB/STB

b intercept of linear function between flow rate, pressure and time

ct total compressibility, M-1L-1T2, psi-1

D tubing inner diameter, L, inch

de,wh effective wormhole radius, L, ft

dperf perforation diameter, L, in.

f(u) relation used in Laplace space to distinguish matrix geometry

ff fanning friction factor, dimensionless

g acceleration of gravitation, LT-2, ft/sec2

gc gravitational dimensional constant

h thickness, L, ft

h reservoir thickness, ft.

hperf perforation spacing, L, ft.

hs coordinate of wellbore location in z-direction

hwh wormhole axial spacing, L, ft.

hx length of the flow field modeled

hz height of the flow field modeled

vi
Iani anisotropy ratio, dimensionless

k permeability, L2, md

K' Fluid-consistency index

Kd Perforation coefficient

kf bulk fracture permeability of dual porosity models, md

kH horizontal permeability, L2, md

km matrix permeability, md

ks damage permeability, L2, md

kV vertical permeability, L2, md

kx permeability in x-direction

ky permeability in y-direction

kz permeability in z-direction

L length, L, ft

l half of fracture spacing, ft.

L general fracture spacing, ft.

l−1 inverse Laplace space operator

Lcore core length, L, inch

Lw length of horizontal well

Lxd left coordinate of wellbore in x-direction

Lxl right coordinate of wellbore in x-direction

m slope of linear function between flow rate, pressure and time

m ̃4 slope of regions 4 (matrix linear flow region)

vii
m(p) pseudopressure (gas), psi-2/cp

mD dimensionless psuedo pressure

mDL dimensionless pressure (rectangular geometry, gas)

mwh number of dominant wormholes per plane

n' flow-behavior index

NAc acid capacity number, dimensionless

Nperf number of perforations

NRe reynolds number, dimensionless

p pressure drop, ML-1T-2, psi

pD dimensionless pressure, dimensionless

pDL dimensionless pressure based on rectangular geometry, linear

pi initial reservoir pressure, ML-1T-2, psi

PVbt,opt optimum pore volume to breakthrough, dimensionless

pWDL dimensionless pressure based on rectangular geometry, linear

pwf bottom flowing pressure, ML-1T-2, psi

q flow rate, L3T-1, ft3/min

qg gas rate, Mscf/day

qw flow rate in the wellbore, L3t-1, bbl/min

rd damage radius, L, ft

rw wellbore radius, L, ft

rw’ equivalent wellbore radius, L, ft

rwh wormhole penetration radius, L, ft

viii
s skin factor, dimensionless

soverall overall skin factor, dimensionless

t time, hours

T absolute temperature, oR

tD dimensionless time

tDAc dimensionless time based on Ac (rectangular geometry)

u laplace space variable

v velocity, LT-1, ft/min

vi interstitial velocity, LT-1, ft/min

vi,opt optimum interstitial velocity, LT-1, cm/min

vi,tip interstitial velocity at the tip of wormholes, LT-1, cm/min

vwh wormhole propagation rate, LT-1, ft/min

w fracture width, ft

xf fracture-half length, ft.

yDe dimensionless reservoir length (rectangular geometry), ft

ye drainage area half-width (rectangular geometry), ft

Z gas compressibility factor

Greek

αz wormhole axial spacing coefficient

γ specific gravity

∆pf frictional pressure drop, ML-1T-2, psi

ix
∆pPE hydrostatic pressure drop, ML-1T-2, psi

∆tsup superposition time function

ε roughness, dimensionless

λAc dimensionless interporosity parameter

μ viscosity, ML-1T-1, cp

ρ density, ML-3, g/cm3

ϕ porosity, fraction

ω dimensionless storativity ratio

Subscripts

i initial

f fracture

m matrix

f+m total system (fracture + matrix)

sc standard condition

sf surface

x
 TABLE OF CONTENTS

Page

ABSTRACT .......................................................................................................................ii

DEDICATION ..................................................................................................................iv

ACKNOWLEDGEMENTS ............................................................................................... v

NOMENCLATURE ..........................................................................................................vi

TABLE OF CONTENTS ..................................................................................................xi

LIST OF FIGURES ........................................................................................................ xiii

LIST OF TABLES ..........................................................................................................xvi

1. INTRODUCTION .......................................................................................................... 1

1.1. Research Background ......................................................................................... 1

1.2. Literature Review ............................................................................................... 6
1.2.1. Matrix Acidizing ........................................................................................... 7
1.2.2. Acid Fracturing ........................................................................................... 18
1.3. Research Objectives ......................................................................................... 21

2. INTEGRATED OPTIMIZATION OF MATRIX ACIDIZING .............................. 24

2.1. Background ...................................................................................................... 24

2.2. Treatment Design ............................................................................................. 25
2.2.1. Horizontal Well Acid Simulator (HWAS) .................................................. 25
2.2.2. Integration of Acid Jetting Function into HWAS ....................................... 26
2.2.3. Model Validation for Jetting Function ........................................................ 28
2.2.4. Integration of Limited Entry Function ........................................................ 30
2.3. Treatment Monitoring ...................................................................................... 32
2.3.1. Real Time Monitoring During Matrix Acidizing ........................................ 32
2.3.2. Integration of Power-Law Fluid Function ................................................... 33
2.3.3. Integration of Gas Reservoir Function ........................................................ 35
2.3.4. Case Study ................................................................................................... 38
2.4. History Matching.............................................................................................. 42
2.4.1. Sensitivity Study ......................................................................................... 43
2.4.2. History Matching Process ........................................................................... 57
xi
 2.4.3. Case Study ................................................................................................... 58
2.5. Integrated Optimization.................................................................................... 60
2.5.1. Optimum Rate Calculation .......................................................................... 60
2.5.2. Case Study ................................................................................................... 62

3. EVALUATION OF ACID FRACTURING THROUGH TREATMENT

MONITORING AND PRODUCTION ANALYSIS ....................................................... 65

3.1. Introduction ...................................................................................................... 65

3.2. Treatment Monitoring in Acid Fracturing in Horizontal Well ......................... 66
3.2.1. Dual-Porosity Transient Slab Model ........................................................... 66
3.2.2. Bilinear Analytical Model ........................................................................... 69
3.2.3. Approach ..................................................................................................... 71
3.2.4. Bilinear Flow Validation by Numerical Simulation ................................... 74
3.2.5. Case Study ................................................................................................... 77
3.3. Production Data Analysis after Acid Stimulation ............................................ 82
3.3.1. RTA Procedure for Gas Well ...................................................................... 82
3.3.2. Case Study ................................................................................................... 84

4. CONCLUSION ........................................................................................................ 89

REFERENCES ................................................................................................................. 91

xii
 LIST OF FIGURES

Page

Fig. 1 Wormholes formed by different injection rates (Fredd et. al, 2000) ..................... 3

Fig. 2 Wormhole Efficiency Curve (Buijse-Glasbergen, 2005) ....................................... 4

Fig. 3 Wormhole growth rate as a function of injection rate ........................................... 8

Fig. 4 Linear core flooding experiments of 1" and 4" in. (Furui at el. 2010)................. 10

Fig. 5 Top and side views from optimal-rate experiments (McDuff et al. 2010) .......... 12

Fig. 6 Top and side views from low-rate experiments (McDuff et al. 2010)................. 13

Fig. 7 Field case results (Paccaloni et al., 1988) ............................................................ 15

Fig. 8 Field case results (Zhu and Hill, 1998) ................................................................ 16

Fig. 9 Pressure and skin evolution match (Buijse et al. 2005) ....................................... 16

Fig. 10 3-D etched profiles for an acid system (Pournik, et. al, 2010) .......................... 19

Fig. 11 Backlit acid-etched width from a laboratory core (Kalfayan, 2007) ................. 20

Fig. 12 Flow regime identification (Samandrli et al. 2014) ........................................... 21

Fig. 13 Square root of time plot (Samandrli et al. 2014)................................................ 21

Fig. 14 Integrated approach for optimum matrix acidizing............................................ 24

Fig. 15 HWAS model scheme (Wu, 2008) .................................................................... 26

Fig. 16 Acid Jetting Assembly (Sasongko et al. 2011) .................................................. 26

Fig. 17 Filter cake removal in segment 2 and wormhole creation in segment 1 ............ 27

Fig. 18 HWAS flow chart with tubing location loop ..................................................... 28

Fig. 19 Wormhole length vs. horizontal section at each treatment ................................ 30

Fig. 20 Spherical flow geometry (Furui at el, 2010) ...................................................... 31

Fig. 21 Schematic of the system for the bottomhole pressure calculation ..................... 33

Fig. 22 Schematic for the horizontal well skin model (Zhu et al., 1999) ....................... 36
xiii
Fig. 23 System schematic of Eq. (2.34) ......................................................................... 39

Fig. 24 Reservoir face pressure with time ...................................................................... 40

Fig. 25 Skin evolution diagnostic plot ........................................................................... 41

Fig. 26 Parametric study on the post-treatment skin effect (Frick et al. 1994) .............. 43

Fig. 27 Wormhole region inside the damaged zone (Tran, 2013).................................. 49

Fig. 28 Pressure-dependent sensitivity on reservoir and near-wellbore parameter ........ 50

Fig. 29 Skin-dependent sensitivity of reservoir and near-wellbore parameter............... 51

Fig. 30 CT scan of wormhole conduit on core-flooding (Furui et al. 2012) .................. 54

Fig. 31 Velocity contour plot on 2D radial flow with mwh = 6. (Furui, 2012) ............... 54

Fig. 32 Pressure-dependent sensitivity on wormhole model parameter......................... 55

Fig. 33 Skin-dependent sensitivity of wormholing model parameter ............................ 56

Fig. 34 Pressure match result ......................................................................................... 59

Fig. 35 Skin match result................................................................................................ 59

Fig. 36 Rate schedule comparison.................................................................................. 63

Fig. 37 Skin evolution comparison ................................................................................ 63

Fig. 38 Schematic of slab matrix linear model well (Bello, 2009) ............................... 67

Fig. 39 Stimulation schematic (Plan view) .................................................................... 71

Fig. 40 Logarithmically spaced LGR grid system ......................................................... 74

Fig. 41 Relative permeability for matrix and fracture .................................................... 76

Fig. 42 Plan view of pressure shows bilinear flow ........................................................ 76

Fig. 43 Normalized rate plot during the injection .......................................................... 77

Fig. 44 Stage 1 treatment pressure & rate ..................................................................... 80

Fig. 45 Diagnostic plot to identify flow regimes. .......................................................... 80

Fig. 46 Evolution of skin and cross-sectional area growth ............................................ 81

xiv
Fig. 47 Production history .............................................................................................. 84

Fig. 48 Recombined gas rate and calculated bottomhole pressure................................. 85

Fig. 49 Normalized rate vs. material balance time ......................................................... 86

Fig. 50 Square root of time plot ..................................................................................... 86

xv
 LIST OF TABLES

Page

Table 1 Data input (Sasongko, 2011) ............................................................................. 29

Table 2 Input data (SPE 134265, Furui at el.) ................................................................ 38

Table 3 Influential parameter list for history matching .................................................. 42

Table 4 Influential parameter list (Buijse-Glasbergen and Furui model) ...................... 44

Table 5 Input data of synthetic well ............................................................................... 45

Table 6 Parameter set-up for sensitivity analysis ........................................................... 45

Table 7 Measurement type of each parameter set .......................................................... 57

Table 8 Wormholing parameter inputs ........................................................................... 58

Table 9 Geometry of the simulation grid data ................................................................ 75

Table 10 Reservoir & simulation data............................................................................ 75

Table 11 Input data for fracture monitoring ................................................................... 78

Table 12 Treatment schedule of stage 1 ......................................................................... 78

Table 13 Wellbore & stimulation data ........................................................................... 87

Table 14 Calculated fracture geometry .......................................................................... 87

Table 15 Comparison summary ..................................................................................... 88

xvi
 1. INTRODUCTION

1.1. Research Background

It was estimated that more than 60% of the world's oil and 40% of the world's gas reserves

are held in carbonate reservoirs before unconventional reservoir came to

commercialization in the past decade. However even though unconventional resource play

has produced more and more oil and gas, it is still important and attractive for the

remaining reserves from carbonate formations to be produced by cost effective manner in

uncertain oil and gas prices. The Middle East, for example, is dominated by carbonate

fields, with around 70% of oil and 90% of gas reserves held within these reservoirs. Some

shale play such as Eagle Ford and Bakken is actually known as carbonate rich shales.

Recent trend shows multi-stage acidizing and acid fracturing in carbonate formations

based on the technology applied from shale play has been recognized and some cases have

been reported in the literature.

Matrix acidizing (called “matrix” because the injection pressure is below the

formation fracture pressure) can be done in either sandstone or carbonate reservoirs, but

since these carbonate formations are highly soluble in acid, matrix acid stimulation is used

as a cost-effective means to enhance well productivity in carbonate reservoirs. The major

goal of acid stimulation under matrix conditions into carbonate formation is to create deep

penetrating conductive flow channels known as wormholes that bypass the damaged near-

wellbore region where there will be virtually no flow in the low-permeability. Effective

1
and highly conductive wormhole penetration beyond the damaged zone should result in a

smaller pressure drop than in the original undamaged formation. Thus, the post-treatment

skin effect could be negative. Actually Furui et al. (2010) showed historical field post-

stimulation buildup test data in various carbonate reservoirs in Middle East and North Sea

fields, where negative -3 to -4 range of skin factors are commonly achieved, which

equivalently means wormhole lengths on the order of 10 to 20 ft. achieved. Very effective

acid stimulation is achievable with a certain completion method.

Carbonate acidizing designs usually consist of 15 - 20 %wt. HCl, which causes a

high surface reaction rate. It is well known that at a given temperature, the ability of a

particular acid to generate wormholes is largely dependent on the acid injection rate. At

low injection rates, acid spends rapidly on the face of the core and no wormholes, or only

short wormholes form. Since, however, at high injection rates, acid flux at the wormhole

is very high, tip splitting and side branching is formed thus wormhole growth rate is

inhibited, as a result wormholes have ramified structures. Dominant wormhole pattern is

obtained at intermediate injection rates. Enough acid reaches the tip to grow a single

wormhole, avoiding excessive side-branching. Compared to the two aforementioned

wormhole patterns, the dominant wormhole structure propagates the longest wormhole

penetration depth with the least amount of acid injected, which is most desired and

efficient injection condition shown in the middle in Fig. 1.

2
 Fig. 1 Wormholes formed by different injection rates (Fredd et. al, 2000)

According to Buijse et al. (2005), under linear core flooding experiment, there

exist three regions characterized by a low, optimum and high interstitial velocity on the

wormhole efficiency curve. Below the optimal flux, dissolution is mostly confined to the

rock face nearest to the acid injection point, and this is called compact dissolution regime,

which should be avoided because longer wormholes will not be formed. Above the optimal

flux, dissolution occur more side branching and called ramified dissolution regime, which

is less efficient compared to the dominant (optimum) wormhole regime. Wormhole

efficiency curve generated through laboratory linear-core flooding experiments as shown

in Fig. 2 indicates that there is a certain optimal flux for which wormholes will most

efficiently propagate along the main axis of the core plug. The injection rate, at which the

dominant wormhole pattern is obtained, is called the optimum injection rate on an acid

treatment. For highly reactive acid/rock systems, the optimum injection rate does exist and

it depends on the rock mineralogy, acid concentration and reaction temperature and other

parameters (Wang et al., 1993).

3
 Fig. 2 Wormhole Efficiency Curve (Buijse-Glasbergen, 2005)

Since wormholes are much larger than the pores in non-vuggy carbonates, pressure

drop through the region penetrated by wormholes is small and can often be neglected.

Thus the wormhole growth is tracked throughout the entire injection period and the

stimulation effect contributed by the wormholes is evaluated by a skin factor.

One of the main challenges in predicting the effectiveness of a carbonate

stimulation treatment is accounting for the wide range of dissolution structures that can be

formed and their impact on skin evolution. Skin varies significantly with dissolution

structure due to changes in the depth of penetration (Fredd, 2000). The costs of matrix

stimulation treatments depend primarily on the volume of injected acid (which influences

treatment time) and the equipment to pump acid into the formation at a certain injection

rate thus optimization should be done by maximizing incremental production obtained by

a matrix acidizing treatment and to balance with cost (Economides at el, 1994). Recent

completion and stimulation methodology for long horizontal wells which are more

4
complicated, affects completion costs (Jackson, 2012). Limited-entry multi-stage

stimulation, pre-perforated liner and isolation packers, acid jetting are the examples.

Although selection of acid placement from various scenarios depends on what the goal of

treatment is, maximum productivity throughout optimum injection condition in field scale

is critical to improve economics.

How to create long wormholes in near-wellbore region through optimum injection

condition (least acid volume) has been challenging but it is still an area of extensive study.

Among lots of wormhole predicting model presented before, it would be key to cover that

a model is simple to apply and allows upscaling from laboratory condition to field scale.

Skin evolution in treatment monitoring can be an important source of information that

helps us characterize the reservoir and understand wormhole growth thus by changing

uncertain parameters, the acidizing forward model used in design phase is able to be

history matched with comparing the real-time observation of simulation results. Practical

guideline to implement optimum acid treatment should be necessary.

On the other hands, acid fracturing is a stimulation conducted above the closure

stress, so that a hydraulic fracture is created. Viscous pad fluid creates an initial crack

geometry, which increases the contact of the well with the reservoir. Then viscous acid is

injected to dissolve the rock along the faces of the fracture and generate etched width

along the fracture. After flushing and injection pressure released, the fracture is allowed

to close by the closure stresses in the formation but it does not close completely due to the

removal of rock at the fracture surface.

5
 Multi-stage acid fracturing in a tight carbonate formation can be an alternative to

propped fracturing as relatively cost-effective stimulation treatment. However the success

of treatments depends on many factors as to whether enough conductivity is secured and

if the conductivity sustains, the selected treatment works under the in in-situ in a specific

geologic environment. Thus, observation and evaluation of past practice is important and

inevitable step to improve stimulation desings. In this study, the methodology to conduct

performance evaluation of acid fracturing treatment using treatment and production

records is described, and a field example is used to illustrate how the procedure works.

Field application of the evaluation procedure shows the effectiveness of the approach.

Actual cross-sectional area growth, using the monitoring program, is shown as fracture

extension and acid etching during the stimulation, which also can be a validation tool for

fracture simulation. With production data analysis, it allows us to compare flowing cross-

sectional area with induced area by the treatment therefore treatment efficiency is

measured. The suggested approach provides engineers with additional information as to

whether the designed acid fracturing was performed appropriately under the in-situ closure

stress field. It is eventually helpful to discuss past practice and improve candidate

selectivity in a company decision making process.

1.2. Literature Review

Literate review for both matrix acidizing and acid fracturing is conducted from the

research perspectives. This literature review summarizes past and current practices and

6
clarifies what was done before. On the basis of the literature review, the research

objectives and approach should be defined in next section.

1.2.1. Matrix Acidizing

The review of matrix acidizing covers laboratory experiment and simulation conducted by

both linear and radial flow geometry, field condition and upscaling from laboratory

condition, treatment monitoring, history matching, and optimal condition.

1.2.1.1. Treatment Design

For treatment design, literature review is categorized into linear core-flooding experiment

and simulation, radial core-flooding experiment and simulation, and field condition and

upscaling. Matrix acidizing has been extensively studied under linear core flooding

condition experimentally. Under laboratory conditions, acid is injected from one end of a

cylindrical core sample at a constant rate and the overall pressure drop is monitored.

Generally, HCl acid is injected axially into carbonate core samples with diameters ranging

from 1in. and lengths ranging from 1in. to 20 in. Acid is pumped at a constant rate until

acid dissolves enough material to break through the opposite end of the core.

1.2.1.1.1. Linear core flooding and simulation

Wang et al. (1993) studied the optimum injection rate in carbonate cores. They found that

there exists optimum rate and reaction pattern in high reactive acid, varying with rock

mineralogy, acid concentration and reaction temperature. Bazin (2001) that reported

optimum injection rate increases with acid concentration, temperature, and limestone
7
permeability. An increase in temperature, concentration and permeability shifts the

optimum flow rate to higher values. The optimum rate became independent of the core

length greater than 20 cm when the wormhole develops in the mass-transport-limited

regime. It was also confirmed that high injection rates are required to increase wormhole

penetration distance thus maximum flow rate provides maximum penetration. Buijse and

Glasbergen (2005) showed in the plot of wormhole growth rate, V wh versus interstitial

velocity, Vi, where optimum injection rate (Vi-opt) associated with the minimum pore

volume to breakthough (PVbt-opt, the smallest volume of acid required for wormhole

breakthrough) does not imply a maximum wormhole growth rate. Rather Vwh is always

an increasing function of Vi, even if Vi is larger than Vi-opt as shown in Fig. 3.

Fig. 3 Wormhole growth rate as a function of injection rate

Cohen et al. (2008) conducted numerical simulation study for different domain

(length is equal to 25cm constant, height varying between 5, 10, 20 and 40 cm) and
8
observed PVbt is in general higher for “confined” conditions (height is equal to 5 and 10

cm) than for “unconfined” domains (height is equal to 40 cm). They found geometry effect

can drastically change the dissolution patterns and optimum conditions, and boundaries of

the core can greatly disturb the wormholing phenomenon by inhibiting the mechanism of

wormhole competition. They demonstrated in the simulation, when shape factor, F (the

ratio of the domain height over the domain length, F equal to 1 would define a square

domain) is more than 1 (unconfined domain, field scale), the flow rate in the wormhole at

optimum conditions becomes independent and stays constant. They concluded

experiments done in confined conditions have probably over-estimated the optimum

injection velocity and the dissolution dynamic is dependent of the core geometry in

unconfined conditions.

Izgec et al. (2008) found that optimal flux in vugular limestone is one to two orders

of magnitude lower than that measured with the same acid formulation in homogeneous

carbonates through the acidizing process. The experiments were conducted with 4-inch

diameter by 20-inch long core samples. This was likely because the acid is flowing through

only a small portion of the rock in the vugular limestone case.

Kalia and Balakotaiah (2009) showed from their simulation study that PVbt

depends on the aspect ratio at low and intermediate injection rate but is independent of the

domain size at high injection rates. They also studied the effect of initial porosity on

pattern formation and reported the amount of acid required to break though is higher for

higher porosity sample because the amount of leak-off from the channel increases as the

local porosity increase. Their study also found that as the heterogeneity magnitude was

9
increased, PVbt decreases, depending on the presence and the connectivity of vugs.

Injection through perforation showed flow expansion in the transverse direction occurs as

soon as it enters the medium, which lead to higher PVbt as compared to the regular

injection case.

Furui at el. (2010) confirmed the consistent trend from the high porosity outcrop

chalk samples that PVbt-opt values in larger core samples are smaller as shown in Fig. 4. In

3D FEM simulation study of comparing 1 in. diameter with 4 in. diameter core effect, the

results showed the interstitial velocity at the tip of the wormhole, Vi,tip-opt increases until

the wormhole penetration length reaches to core diameter and after reaching the length,

the tip velocity of the wormhole becomes relatively constant, which explains why the

apparent wormholing efficiency is greater for larger diameter core samples.

Fig. 4 Linear core flooding experiments of 1" and 4" in. (Furui at el. 2010)

Kai at el. (2013) studied the effects of core length and core diameter on the

optimum condition experimentally and found that a stable value of the optimal acid flux
10
was obtained in 6 inches or longer cores. Because there may be more than one wormhole

formed close to the entry, the wormhole competition effect reduces the tip velocity of the

wormholes in short core length. However, when the core is long enough, a relatively

constant tip velocity of the wormholes can be maintained, and the wormhole propagates

more efficiently.

1.2.1.1.2. Radial core flooding and simulation

According to Buijse et al. (2005), the basic physical and chemical principles underlying

the wormhole growth process are the same in linear and in radial geometry. However a

fundamental difference between linear and radial geometry is the dependence of Vi on

position. In linear geometry, Vi is independent of the position in the core and consequently,

Vwh depends only on the injection rate, and not on the position of the wormhole front in

the core. However, in radial geometry, Vi decreases as the distance R from the wellbore

increases. Thus longer, deeply penetrating wormholes would receive less acid at the tip,

and consequently would grow at a lower rate, compared to short wormholes.

Mostofizadeh and Economides (1994) conducted radial core flood experiments

and reported that reaction pattern was consistent with those obtained by linear core-

flooding conducted by Wang et al. (1998). However, the pattern are far more complicated

and multi-faceted in radial geometry. It was also reported that the optimum acid injection

rate depends on the mineralogy and morphology of the reservoir and on the transition point

from diffusion-limited to fluid loss-limited modes, which is also related to the acid

concentration and temperature. Optimum injection rate highly depends on permeability.

11
For low permeability rock, optimum injection rate should be lower but rock with

increasing permeability would require higher optimum injection rate because the area-to-

volume relation increases and the highly unstable nature of the acid attack on the rock is

enhanced.

McDuff et al (2010) conducted larger-scale. CT scan image of two experiments

conducted at near-optimal rate and at below-optimal rate with Indiana limestone using

15% HCl acid are shown in Fig. 5 and Fig. 6, where the experimental condition

corresponds roughly to pumping 40 bpm into 5000 ft. of completion interval and

generating 10-25 ft. wormholes into the formation. Nearly four times the acid volume was

required in the low-rate condition compare to the near-optimal rate experiments.

Fig. 5 Top and side views from optimal-rate experiments (McDuff et al. 2010)

12
 Fig. 6 Top and side views from low-rate experiments (McDuff et al. 2010)

Cohen et al. (2008) conducted numerical simulation study about radial geometry.

Their main observation was that dissolution close to the optimum injection velocities is

achieved in radial flow at higher injection velocity than in linear flow.

1.2.1.1.3. Field condition and upscaling from lab condition

As evidenced by many researchers describing the optimum injection rate measured in a

linear core-flooding experiment cannot be translated easily to field conditions, where flow

is very likely to radial flow. In radial geometry, optimal wormhole-growth rate in radial

flow does not come from a single pump rate, but changes as the formation is penetrated

by the wormhole front. Longer wormholes require a higher pump rate to grow efficiently

thus pump schedule should start at relatively low pump rates, and increase the pump rate

during the treatment. It actually requires real time knowledge of wormhole growth rate

and penetration depth, which is difficult to obtain.

13
 Buijse et al. (2005) and Glasbergen (2009) recommended the best practice is to

pump at the maximum rate possible below the fracturing pressure. Glasbergen (2009)

reported several other factors will play a role under field conditions. Injection rate

limitations because of fracture pressure or equipment limitations, heterogeneities in

mineralogy, variations in injection temperature, variations in injection rates including

shut-ins, the effect of diverters and wellbore effects such as travel time for the acid from

the heel to the toe is the examples.

Furui et al. (2010) developed the modified wormhole growth model based on

Buijse and Grasbergen’s empirical correlation. The new wormhole model estimates the

wormhole evolution with the consideration of acid flux at the tip of the wormhole,

overcoming core size dependencies, which allows upscaling from laboratory linear core-

flooding on small cores to radial/spherical flow geometries in field-size treatment.

1.2.1.2. Treatment Monitoring

Paccaloni et al. (1988) introduced a method that used instantaneous pressure and rate to

calculate skin factor continuously. This method is based on the concept of finite-radius

“acid bank” and steady-state, radial Darcy’s flow. Fig. 7 illustrates skin evaluation for the

field case using Paccoloni’s method.

14
 Fig. 7 Field case results (Paccaloni et al., 1988)

Another technique by Prouvost and Economides (1989) proposed to analyze skin

evolution by continuous comparison of measured and simulated pressures. The technique

requires simulation of the transient pressure response to the injection of inert fluids. The

flow rates used for the inert fluids follow the exact injection schedule of the acid treatment.

The difference between the actual bottomhole pressure and simulated bottomhole pressure

is utilized to evaluate the changing skin factor. Hill and Zhu (1996, 1998) proposed a

method based on the theory of standard injectivity test using approximate line source

solution for transient flow to monitor changing skin during matrix acidizing treatment (Fig.

8). Furthermore, this approach employs a superposition method to account for the transient

flow effects due to injection of acid at multiple rates and pressures. Therefore, each point

on the inverse injectivity vs. superposition time plot will lie on a straight line having slope,

m, with its intercept depending on the skin factor at time t.

15
 Fig. 8 Field case results (Zhu and Hill, 1998)

For recent practical example, Kent et al.(2014) showed the skin evolution treated

by multi-stage intelligent completion. The evolution of effective wellbore radius through

the chalk formation in North Sea was monitored during injection. Skin was -4.4 at the end

of the treatment, showing a successful stimulation treatment.

1.2.1.3. History Matching

Buijse et al. (2005) conducted a field case study using acid placement model based on

their wormhole model. Fig. 9 shows measured and simulated pressures and skin during

the stimulation, having a good agreement.

Fig. 9 Pressure and skin evolution match (Buijse et al. 2005)

16
 Varun et al. (2007) conducted pressure history matching using their developed acid

placement model. The reservoir permeability and the PVbt in the volumetric model were

adjusted to obtain a match of the actual treating pressure. Furui et al. (2010) showed in a

case study of downhole injection pressure matching results during an acid injection where

measured PVbt,opt and Vi,opt were used based on the linear core-flooding experiments.

Simulation with conventional radial flow showed the downhole injection pressure

becomes too high because Vi,tip is lower and calculated wormhole penetration depth is too

short. In the upscaling method, the wormholing axial spacing, αz was changed until a good

match and final value was αz was 0.75, where dominant wormhole with some spacing

allows much higher Vi,tip at wormhole tip than conventional radial flow condition.

1.2.1.4. Optimum Condition

Glasbergen et al. (2009) described the constraints of optimum injection rate for wormhole

propagation on the basis of excellent literature review and case studies. According to their

study, to be a candidate for applying optimum injection rate, the reservoir pressure should

be fairly uniform, the total zone height should be relatively short, and no major variations

in permeability should be present. Although the optimum injection rate does exist for HCl,

it is necessary to know the injection flow distribution. Real-time measurements should

help determine the flow distribution and can potentially be used to determine the best flow

rate or maximum injection rate. It is recommended to increase flow rate during a treatment

and eventually apply maximum available pump rate below fracture pressure in

combination with diverters when limited information of reservoir is available or reservoir

17
pressure variations are expected. Thus application of the flow rate lower than the

maximum allowed flow rate can be considered only when (1) the optimum flow rate is

fully understood under radial conditions; (2) the flow distribution during a treatment is

known at all times; and (3) the theoretical optimum flow rate is significantly lower than

the maximum allowed flow rate.

1.2.2. Acid Fracturing

1.2.2.1. Treatment Design and Monitoring

The design of an acid fracturing treatment is accomplished by estimating the optimum

conductivity and acid penetration distance that results in maximum benefit of the treatment.

Design parameters include selecting the fluid types, number of stages, pumping rate, and

injection time. Changing these parameters results in different fracture geometry, etching

patterns, and acid-penetration distance. A complete study of fluid properties, mineralogy

and permeability distributions, and formation temperature should be conducted prior to

the stimulation. Simulators are usually used to estimate how these design parameters affect

the stimulation job (Jawad, 2014).

The effects of acid solutions injected into hydraulic fractures created in carbonate

formations can be assessed at the laboratory scale in acid fracture conductivity test that

mimic the conditions in an actual acid fracture treatment (Pournik, et. al, 2010). Fig. 10

shows that 3D images of the etching pattern by a certain acid. The etching profile shows

the difference between the original surface and after the acid etching process. The etched

fracture surface volume was calculated from the difference in surface volume between the
18
before and after acidizing. From the etched surface volume, the etched surface width is

calculated using the cross sectional area of the fracture.

Fig. 10 3-D etched profiles for an acid system (Pournik, et. al, 2010)

1.2.2.2. Production Data Analysis

Evaluations of acid-fracturing treatments are usually based on a production increase or a

comparison with other wells. These evaluations determine the relative success of materials

and techniques compared with other materials and techniques (Elbel, 1994). The analysis

of post-fracture treatment is conducted by pressure-buildup data and/or production data.

Production from acid fractures differ significantly from propped fractures. Although

created open fractures have an almost infinite conductivity, taking into account in-situ

stresses and formation strength affects the post-stimulation performance significantly

(Ben-Naceur, Economides, 1989). If etching is non-uniform, the fracture may close with

conductivity retained, as depicted in Fig. 11 (Kalfayan, 2007).

19
 Fig. 11 Backlit acid-etched width from a laboratory core (Kalfayan, 2007)

For flow regime identification, linear flow can be detected by 1/2 slope line in log-log

plots of either pressure drop or reciprocal of production rate versus time as shown in Fig.

12. Useful plot is the square root of time plot in which p or 1/q data is plotted versus.

√𝑡 as shown in Fig. 13. For gas well, pseudo pressure should be used and a straight line

of Δm(p)/qg vs. √𝑡 plot, either for the constant qg production or the constant pwf production

cases is equivalent to the half-slope period in a log-log diagnostic plot. We can then record

the slope of the straight line and the actual time when the boundary is reached. We use the

slope of the line and end of half-slope time to calculate√𝑘𝐴𝑐 . If matrix permeability is

identified, cross-sectional area flowing from fractured wall, Ac can be calculated.

20
 Fig. 12 Flow regime identification (Samandrli et al. 2014)

Fig. 13 Square root of time plot (Samandrli et al. 2014)

1.3. Research Objectives

Through the literature review in the previous section, it is clear that how to reduce

uncertainty is still the key to achieve optimal treatment. The research goal is to establish

21
practical guideline to evaluate both matrix acidizing and acid fracturing through integrated

approach based on treatment monitoring.

For matrix acidizing, it can be done by combining design, monitoring, history

matching as a closed loop. In this study, the latest Furui’s augmented wormhole model

(2010) on the basis of Buijse and Glasbergen correlation model (2005) is adopted because

according to the literature review, Furui’s model includes the consideration of acid flux at

the tip of the wormhole, overcoming core size dependencies, which allows upscaling from

laboratory linear core-flooding to radial/spherical flow geometries in field-size treatment.

However uncertain parameters still exist in design phase. Skin evolution available during

treatment monitoring is used for history matching process.

For acid fracturing, the method to evaluate acid fracturing is through treatment

monitoring and production data analysis. New method for real-time monitoring of acid

fracturing, the inverse injectivity vs. superposition time function plot is proposed. By

monitoring cross-sectional area from injection flow, real-time area growth induced by acid

fracturing is monitored. Comparing the stimulated area with the result from linear flow

diagnostic approach from production data, the treatment efficiency can be measured. It

provides engineers with additional information as to whether the designed acid fracturing

was performed appropriately under the in-situ closure stress field.

In the summary, the goal of this work is divided into two:

1. Matrix Acidizing:

22
 1) To integrate the related models described in the literature. Existent horizontal

well acid simulator, Acid Jetting model, and skin monitoring program are

consolidated as one integrated tool package.

2) To conduct parameter sensitivity study in order to identify influential

parameters in Furui wormhole model.

3) To develop a consistent procedure that can be used from treatment monitoring

to history matching so that the result is used to optimum condition study.

2. Acid Fracturing:

1) To develop new monitoring method

2) To develop a consistent procedure that can be used from treatment monitoring

to production analysis to measure treatment efficiency.

23
 2. INTEGRATED OPTIMIZATION OF MATRIX ACIDIZING

2.1. Background

In this chapter, all related existent programs and models are integrated to one package for

matrix acidizing so that one can use it 1) to design a treatment based on the well structure

and completion, 2) to monitor stimulation efficiency during treatment and 3) to conduct

history matching the monitored data during treatment and the designed data to update

reservoir and key wormholing parameters and 4) to optimize future well stimulation using

the updated parameters. This procedure is illustrated in Fig. 14.

Fig. 14 Integrated approach for optimum matrix acidizing

24
2.2. Treatment Design

The existing numerical acidizing simulator for horizontal well used in this study was

originally developed by the Mishra and Furui (2007). The modification was made by

Nozaki (2009) for vertical well in gas reservoir where the viscous diversion effect was

included. Tran (2013) added Furui’s wormhole model and apparent skin factor model for

viscous diversion for radial flow in horizontal well in openhole condition. In this section,

further functions have been added to deal with more completion such as acid jetting and

limited entry method.

2.2.1. Horizontal Well Acid Simulator (HWAS)

HWAS is an horizontal well acid stimulation simulator. The simulator shown in Fig. 15

consist of a wellbore flow model, a wellbore fluid interface tracking model, a transient

reservoir outflow model, a wormhole growth model and a skin model. The wellbore flow

model accounts for pressure drop and material balance inside the wellbore. The fluid

interface tracking monitors the interface between the injected fluids in the horizontal

wellbore. The transient reservoir flow model captures the transient effect of varying

injection rates that are often seen in well test problems. The wormhole model predicts the

wormhole penetration in the formation during the entire acid injection period. The

apparent skin model accounts for well completions damaged region, wormholes, reservoir

mobility and injected fluids mobility. Final skin factor is a function of wormhole

penetration depth with the assumption that wormholes extend beyond damage zone at the

25
end of the treatment (Tran, 2013). More details are found in the literature (Mishra and

Furui, 2007, Tran, 2013).

Fig. 15 HWAS model scheme (Wu, 2008)

2.2.2. Integration of Acid Jetting Function into HWAS

Separately, Sasongko (2011) developed an acid placement model by acid jetting out of a

drill pipe. Acid is jetted onto the face of openhole wellbore as the drill pipe is withdrawn

from the well. The acid is pumped through the nozzle holes giving a jetting effect around

the wellbore (Fig. 16).

Fig. 16 Acid Jetting Assembly (Sasongko et al. 2011)

26
The jetting action helps to remove the drilling fluid filter cake and promote the acid to

penetrate into the formation and form wormholes to stimulate the well. The model

simulates the acid jetting process using a comprehensive model of acid placement and

wormhole propagation in a horizontal well. The following steps are taken.

• Openhole horizontal section is divided into a series of segments and each segment

represents one cycle of the jetting treatment.

• Mechanical action by acid jetting removes the filter cake around the wellbore

section by section. After filter cake is removed, acid goes to the segments and wormhole

is created, meanwhile jetting starts to remove the filter cake in the next section as shown

in Fig. 17.

• The injection acid amounts are controlled by the duration of injection in each cycle,

and the pumping rate and the drilling string pulling speed are set close to field operation

conditions.

Fig. 17 Filter cake removal in segment 2 and wormhole creation in segment 1

Fig. 18 shows the flow chart with jetting operation. Notice that jetting is an individual

module added to tubing location step.

27
 Fig. 18 HWAS flow chart with tubing location loop

2.2.3. Model Validation for Jetting Function

Using same example (Sasongko 2011), integrated jetting function in HWAS was run and

validated with the case study result. The input data is shown in Table 1. Acid is injected

at target volume coverage of 0.5 bbl/ft with a pulling speed of the drill pipe is about 50

ft/min. For the wormhole model, the example case uses a value of 0.53 for optimum pore

volume to breakthrough, PVbt-opt, and value of 1.75 cm/min for optimum interstitial

velocity, Vi-opt. These values were obtained from a laboratory experiment using 1” x 6”

core with 15% HCl at 15 F (Furui et al. 2010).

28
 Table 1 Data input (Sasongko, 2011)

As the result of the simulation, it was confirmed that the same plot shown in Fig. 19 was

generated compared with the one generated by original acid jetting model. Besides, the

simulation time was shrunk dramatically. Compared to 3 minutes running time in original

model, integrated function in HWAS took only 5 seconds to finishing calculation. The

integrated approach was successfully completed.

According to Sasongko’s description, temporarily plugging off the toe end of the

well is necessary to increase the flux in the section towards the heel, since the study result

shows the toe section is the most stimulated zone and most of the well (up to 75%) receives

almost no matrix stimulation other than filter cake removal. Falling to maintain sufficient

acid flux into the formation as longer and longer sections are exposed can result in very

inefficient matrix stimulation. To maintain wormhole growth throughout such a treatment,

the acid flux into the formation cannot drop too far below the optimal flux value.

29
 Fig. 19 Wormhole length vs. horizontal section at each treatment

2.2.4. Integration of Limited Entry Function

Limited entry is the process of limiting the number or reducing the entry-hole diameter of

perforations in such a way that significant perforation friction pressure is achieved during

the treatment. Perforation friction establishes a backpressure in the wellbore that tends to

allocate flow among the multiple perforation intervals/clusters, thus improving control of

fracturing process. Spherical flow geometry is important for limited entry, especially in

the formations in the near-wellbore region until pressure interference among the

perforations occurs (Furui at el, 2010). As seen in Fig. 20, spherical flow near perforation

30
region disappears due to no flow boundary at the half-distance between each perforation

and the flow pattern approaches to radial flow.

Fig. 20 Spherical flow geometry (Furui at el, 2010)

From program perspectives, the implementation of the transition has been given as

follows.

h perf
(i) rwh  (Spherical flow)
2

q
vi ,tip  (2.1)
4d e,wh rwh

h perf  h perf 
s  ln   1 (2.2)
2rwh  2rwh 

h perf
(ii) rwh  (Radial flow)
2

q  1  1 
vi ,tip  (1   Z )  Z  
L mwh   d  (2.3)
de, whrwh  e, wh 

31
 r 
s   ln wh  (2.4)
 rw 

2.3. Treatment Monitoring

The concept of inverse injectivity is applied to calculate skin evolution using rate/pressure

data during the matrix acidizing treatment. Pandya (2012) developed skin monitoring

program in his study. Skin monitoring model has been implemented to HWAS.

2.3.1. Real Time Monitoring During Matrix Acidizing

Hill and Zhu (1996) proposed a method based on the theory of standard injectivity test

using approximate line source solution for transient flow to monitor skin during matrix

acidizing treatment. The pressure response to multiple injection rates for a vertical well

is given by:

p wf  pi
 mt sup b (2.5)
qi

where,

162.6 B
m (2.6)
kh

  k  
b  m log 
2 
 3.23  0.87 s  (vertical well) (2.7)
  ct rw  

q j  q j 1
logt N  t j 1 
N
t sup   (2.8)
j 1 qN

The slope, m, in Eq. (2.6) remains constant as reservoir parameters do not change during

acidizing treatment. The only changing parameter is the skin factor, s the equation for
32
intercept, b Eq. (2.7). Furthermore, this approach employs a superposition method (Eq.

(2.8) (Earlougher, 1977) ). Using Eq. (2.6) through Eq. (2.8), skin factor can be calculated

in real-time from measured pressure, injection rate, and time during an acid treatment from

the equation given below:

1  b  k  
s  log   3.23
2 
(2.9)
 
0.868  m  ct rw  

2.3.2. Integration of Power-Law Fluid Function

In order to use the concept of inverse injectivity, it is required to calculate bottomhole

pressure. In most acid treatments, bottomhole pressure is not measured and only surface

pressure is recorded at the injection tubing or the annulus. Fig. 21 shows the schematic of

the system used for the calculation of bottomhole pressure.

Fig. 21 Schematic of the system for the bottomhole pressure calculation

In the system, the surface pressure can be converted to the bottomhole pressure by

p wf  p sf  p PE  p f (2.10)

33
where, psf is the surface pressure, pwf is the bottomhole flowing pressure, ΔpPE is the

hydrostatic pressure drop, and Δpf is the frictional pressure drop. In case the surface

pressure is measured in the annulus of the well, the frictional pressure drop is zero and

only hydrostatic pressure drop is used to calculate the bottomhole pressure.

For single phase liquid, the hydrostatic pressure drop depends only on the density

of the fluid and the height of the fluid column. Therefore, the hydrostatic pressure drop

changes when a fluid with different density is injected into the tubing. The hydrostatic

pressure drop can be calculated by

g cos N
p PE 
gc A
  V
i 1
i i  qt new  (2.11)

where, A is the cross-sectional area of the tubing, q is the injection rate, θ is the average

inclination of the tubing, ρi-1 is the density of the fluid in the tubing, ρi is the density of

the fluid being pumped, Vi is the cumulative injected volume of the i-th fluid, L is the

height of fluid of the tubing, and Δtnew is the time increment after start of pumping the new

fluid.

Considering viscous acid is used, the power-law fluid should be added to the

program. Basically the frictional pressure drop depends on the injection rate, fluid density,

and fluid viscosity, which may vary during an acid treatment. The friction pressure drop

is determined from the Fanning equation,

2
q
2 f f   L
 A (2.12)
Pf 
gc D

34
where, ff is the Fanning friction factor, and D is the diameter of the tubing, In the above

equation, the Fanning friction factor depends on the Reynolds number and is explicitly

calculated by (Chen, 1979)

1 16 if NRE < 2000

 (2.13)
ff N RE

1 
  5.0452   1.1098  7.149 
0.8981
 
 4 log  log      if NRE > 2000 (2.14)
ff 

3.7065 N RE  2.8257  RE 
N  


For Power-law fluid, Metzner-Reed generalized Reynolds number equation is used in

oilfield units as follows:

0.249u 2 n ' D n '

N Re 

96 n ' K ' 3n '  1 / 4n' n' (2.15)

The function of selecting “Power-law” fluid and “K’ & n’ ” in addition to Newtonian fluid

was added to the program to calculate bottomhole pressure from surface pressure in

viscous fluid.

2.3.3. Integration of Gas Reservoir Function

Fig. 22 shows the schematic of the system for the skin model for horizontal gas well. The

system consists of a horizontal wellbore that is bounded in x-direction and z-direction.

However, the system is not bounded in the y-direction. Therefore, Goode and

Thambynayagam (1987) referred to this model as the semi-infinite slab model.

35
 Fig. 22 Schematic for the horizontal well skin model (Zhu et al., 1999)

As for gas reservoir, semi-slab model has been modified for gas case as follows using gas

dimensionless pseudo pressure:

khTsc
mD  (mi  m) (2.17)
50,300qpscT

In injection, since pressure range is normally more than 3,000 psi, the pseudo pressure can

be converted into single pressure:

khTsc  2 pi 
pD  ( p i  p )  (2.18)
50,300qpscT  i Z i 

Now superposition time function is written as usual:

q j  q j 1
 t D, N  t D, j 1 
N
t sup   (2.19)
j 1 qN

Since intercept b is obtained in pressure and time monitoring:

pwf  pi
b  m  t sup (2.20)
qi

Slope m and intercept b is defined:

2  50,300 p scT i Z i rw '

m (2.21)
2 pi Tsc hx hz k y

36
 hx2 N q j q j 1 hxhz N q j q j 1 50,300 p scT i Z i 1 (2.22)
bm
 2vx
j 1 qN
m
Lw v z 
j 1 qN

2 pi Tsc k y k z Lw
s

The remaining parameters used in Eq. (2.22) are defined:

 

1 
    erf vxn t D  t D , j 1  2n  (2.23)
n1  n 

 

1 
    erf v z l t D  t D, j 1  l coslz e  (2.24)
l 1  m 

1   nL xl   nL xd 
n  sin   sin  (2.25)
nLw   h x   hx 

1   l   l 
l  sin hs  2rw '  sin hs  2rw ' (2.26)
4lrw '   h z   hz 

ze 
1
hs  1.47rw ' (2.27)
hz

0.000264k y t
tD  (2.28)
ct rw '

rw ' k x
vx  (2.29)
hx k y

rw ' k z
vz  (2.30)
hz k y

Lw  Lxl  Lxd  (2.31)

1/ 4
k 
rw '  rw  z  (2.32)
k 
 y 

Finally, skin evolution is given as below:

2 pi Tsc k y k z Lw  hx2 N q j q j 1 hxhz N q j q j 1 

s b  m 2
50,300 p scT i Z i   vx
 m
Lw v z 
  (2.33)
j 1 qN j 1 qN 

37
 In the above equations, the geometry functions and dimensionless groups are as

follows. The infinite summation terms are approximated by the first 40 terms for a stable

result. The drawback of this skin model is that it does not consider the effect of reservoir

heterogeneity. The variation in permeability and skin along the horizontal wellbore is

interpreted using three constant values of permeability (kx, ky, and kz) and one skin value.

Thus, the horizontal skin model provides a global estimate of skin evolution and will not

provide the skin profile along the horizontal length.

2.3.4. Case Study

Limited entry case study is shown using same input as Furui et al (2010) and the actual

skin evolution has been generated. Given input parameters are shown in Table 2.

Table 2 Input data (SPE 134265, Furui at el.)

Well & Stimulation Data
Parameters Input Values
Stimulation length, ft (Stage1) 970
Perforation diameter (in) 0.21
Perforation shot density (SPF) 0.1
Perforation interval (ft) 10
Top perf depth (ft) –MD / TVD 15,173 / 9,632
Acid concentration (wt%) 28
Acid volume (bbl /ft) 1.65

In this case, bottomhole pressure were recorded by the downhole gauge. Since reservoir

face pressure is necessary to conduct skin calculation appropriately, hydrostatic, frictional

38
and perforation frictional pressure drop has been subtracted from the measure pressure

(Fig. 23).

presevoir  p gauge  pPE  p f _ pipe  p f _ perf (2.34)

1.98244q 2  (2.35)
p f _ perf  2
N Perf K d2 d Perf
4

where, q is the injection rate in bpm, γ is the fluid specific gravity, Nperf is the number of

perforations, Kd is the discharge coefficient, dperf is the diameter of the openings in inch.

Usually the value of Kd between 0.6 and 0.8 is used.

Fig. 23 System schematic of Eq. (2.34)

With fitted friction reducer multiplier, calculated reservoir face pressure is shown in Fig.

24. When the reservoir face pressure drop below closure pressure, matrix acidizing occurs.

39
 9000 120
Frac Extension Pressure
Closure Pressure
8000 100

Rate (bpm), FR Multiplier (%)

7000 80
Pressure (psi)

6000 60

5000 40

4000 20
Downhole Gauge Pressure Reservoir Face Pressure
Injection Rate Friction Reducer Multiplier
3000 0
0 5 10 15 20 25
Elapsed time (min)

Fig. 24 Reservoir face pressure with time

With reservoir face pressure estimated, skin evolution is calculated. The result is shown

in Fig. 25. Because the reservoir was hydraulically fractured in the beginning the treatment,

initial skin was started by a negative value. Considering limited entry completion, early

flow regime is spherical flow but the effect should be dissipated after created wormhole

length is extended more than 5ft (the half of the perforation spacing). In this case, around

3 minutes are identified as transition point from spherical flow to radial flow since -3.0

skin is equivalent to 5ft wormhole radius. The wormhole continue to extend in the radial

flow condition and the final obtained skin was -4.5 which is equivalent to 20ft wormhole

length created by matrix acidizing. Through this case study, the program successfully

follows the limited entry case and generate reasonable skin evolution result.

40
 0
-0.5
-1
rw
-1.5
rwh
-2 lf
Skin

-2.5
-3
rw
-3.5 rwh
lf
-4
-4.5 Skin = -2.3 Skin = -4.5
-5
0 5 10 15 20 25
Elapsed time (min)

Fig. 25 Skin evolution diagnostic plot

On-site monitoring by skin evolution to evaluate treatment efficiency has been

successfully used during stimulation, which can be an important information that helps us

to characterize the reservoir and understand wormhole growth. Through history matching

process described in next section, uncertain parameters used in the well placement

simulator with a wormhole model is able to be modified by comparing the real-time

observation with the simulation results.

41
2.4. History Matching

Now that actual and simulated pressure and skin with time is obtained, next is to history

match the two sets of data to evaluate the influential parameters. Table 3 shows the

influential parameters list used in Furui wormhole model.

Table 3 Influential parameter list for history matching

Reservoir & Near-Wellbore Stimulation / Wormhole

 Porosity,   Optimum Pore Volume to Breakthrough, PVbt-opt

 Horizontal Permeability, kH  Optimum Interstitial velocity, Vi-opt
 Permeability Impairment ratio, kS/k  Effective wormhole diameter, de,wh
 Permeability Ratio, kH/kV  No. of Dominant Wormholes, mwh
 Damage Penetration Radius, rd  Wormhole axial-spacing coefficient, z

First of all, we need to reduce the number of parameters in the list. Formation

evaluation technique by well-logging analysis and pressure fall-off test are highly

recommended to obtain the best-estimated reservoir and near-wellbore parameters such as

porosity, permeability, damaged zone radius and permeability before designing an acid

treatment. Well completion also matters. Although the flow is converged to radial flow at

some point, whether the well is perforated and/or fractured before matrix acidizing may

give different flow regime rather than radial flow during the early time of acid treatment.

In this study, a sensitivity study using Furui wormhole model was conducted for

reservoir and near-wellbore parameters as wells as the wormhole model parameters. The

understating of each parameter’s effect on pressure and skin evolution behavior during

acid stimulation is crucial to further optimization study.

42
 2.4.1. Sensitivity Study

Through the sensitivity study, the influential parameter and physical meaning of Furui

wormhole model on the bottomhole pressure and damaged skin versus time is described.

Economides et al. (1994) conducted parametric studies on the post-treatment skin

effect. The impact of various factors on the skin effect, including acid volume, ratio of

undamaged to damaged permeability, porosity of the original formation, fractal dimension,

and injection rate was plotted using their wormhole model (Fig. 26).

Impact of permeability ratio on post-treatment skin Impact of porosity on post-treatment skin

Impact of the fractal dimension on post-treatment skin Impact of acid injection rate on post-treatment skin

Fig. 26 Parametric study on the post-treatment skin effect (Frick et al. 1994)

43
Following the literature, the parametric study conducted in this section plots the impact of

various factors on the skin evolution as well as bottomhole pressure over the treatment

time, using Furui wormhole model. Table 4 is the parameter list used in HWAS including

both Buijse-Glasbergen model and Furui wormhole model. It is recognized that 3 more

parameters are used in Furui wormhole model.

Table 4 Influential parameter list (Buijse-Glasbergen and Furui model)

Buijse-Glasbergen Model Furui Model

[Reservoir & Wellbore] [Reservoir & Wellbore]

 Porosity, ϕ  Porosity, ϕ
 Horizontal Permeability, kH  Horizontal Permeability, kH
 Permeability Impairment ratio, kS/k  Permeability Impairment ratio, kS/k
 Permeability Ratio, kH/kV  Permeability Ratio, kH/kV
 Damage Penetration Radius, rd  Damage Penetration Radius, rd

[Stimulation / Wormhole] [Stimulation / Wormhole]

 Optimum Pore Volume to  Optimum Pore Volume to
Breakthrough, PVbt-opt Breakthrough, PVbt-opt
 Optimum Interstitial velocity, Vi-opt  Optimum Interstitial velocity, Vi-opt
 Effective wormhole diameter, de,wh
 No. of Dominant Wormholes, mwh
 Wormhole axial-spacing coefficient, z

In this sensitivity study, synthetic horizontal well in a homogeneous limestone

reservoir under openhole and radial flow geometry condition (Table 5) is used. Each

parameter is varied with 3 values, low, mid, and high. Simulation using HWAS and Furui

wormhole model is performed repeatedly with no change to the other parameters in order

to quantify the sensitivity to the parameter changed.

44
 Table 5 Input data of synthetic well
Summary of Well and Stimulation Data
Parameters Input Values
Initial Reservoir Pressure (psi) 3000
Fluid Viscosity (cP) 0.50
Formation Thickness (ft) 150
Total Compressibility (1/psi) 3.50E-06
Stimulation Length (ft) 1000
Wellbore Radius (ft) 0.328
Tubing OD (in.) 3.00
True Vertical Depth (ft) 7000
Acid Concentration (wt%) 15
Acid Volume (bbl/ft) 0.6
Injection Rate (bpm) 20

Table 6 Parameter set-up for sensitivity analysis

Base values of each parameters
Base
Parameter Low - - High
Values
Porosity, ϕ 0.05 - 0.15 - 0.3
Horizontal Permeability, kH (md) 12 - 15 - 18
Permeability Impairment ratio, kS/k 0.2 - 0.3 - 0.4
Permeability Ratio, kH/kV 5 - 10 - 15
Damage Penetration Radius, rd (ft.) 0.5 - 1.0 - 2.0
Optimum Pore Volume to Breakthrough, PVbt-opt 0.45 - 0.85 - 1.25
Optimum interstitial velocity, Vi-opt (cm/min) 0.5 - 1.0 - 1.5
Effective wormhole diameter, de,wh (in.) 0.1 - 0.25 - 0.5
No. of Dominant Wormholes, mwh 3 - 6 - 12
Wormhole axial-spacing coefficient, z 0 - 0.5 - 1

2.4.1.1. Reservoir And Near-Wellbore Parametric Sensitivity Study

The results of pressure-dependent sensitivity for reservoir and near-wellbore parameters

are summarized in Fig. 28. It is observed that formation damage related parameter such

45
as permeability impairment ratio (ks/k) and damage radius (rd) are sensitive during early

injection time until wormhole radius exceeds damaged zone. As an overall trend, porosity,

horizontal permeability and permeability ratio (kh/kv) are sensitive during entire

stimulation period.

On the other hand, skin sensitivity for reservoir and near-wellbore parameters are

shown in Fig. 29. It is observed that the overall behavior of skin sensitivity is same as

pressure sensitivity. Because the nature of skin is an additional pressure drop in radial flow

condition to count for production efficiency deviated from the ideal condition (Furui at el

2003), it is not surprising the parameters affect pressure and skin in a similar trend.

However horizontal permeability and permeability ratio (kh/kv) are not sensitive parameter

because of the way that skin is calculated. The description of each parameter effect on

pressure and skin is as follows.

(a) Porosity:

With the superposition principle, the reservoir flow estimation to include the transient

effects during acid injection process can be estimated from Lee et al. (2003);

2kl
 p R  pw    q j p D (t n  t j 1 ) D   qn sn
n
 (2.36)
 j 1

where,

q j  q j  q j 1 n (2.37)

4.395 106 kt
tD  (2.38)
ct rw 2

46
 1
pD  (ln t D  0.80907) (2.39)
2

In HWAS, reservoir flow rate and pressure at each time step in radial geometry is

solved simultaneously because wellbore material balance is coupled with reservoir flow

model. Eq. (2.36), (2.38) and (2.39) indicates wellbore pressure decreases as porosity

increases for injection. Skin deceases because of pressure drop. However, in the wormhole

velocity described in radial geometry in Eq. (2.40) and (2.41) in Furui’s model, Vi,tip

decreases, resulting in a slower Vwh and skin reduction as porosity becomes larger. Thus

the difference of 3 cases are minimal.

q  1  1 
vi ,tip  (1   Z )  Z  
L mwh   d  (2.40)
de, whrwh  e, wh 

2

   v PV 2
 (2.41)
 v i ,tip PVbt ,opt N AC   
v wh  v i ,tip N AC    1  exp 4
i ,tip bt , opt N AC Lcore
  
 v i ,opt     v i ,opt rwh  
     

(b) Horizontal permeability:

Transient reservoir flow in Equation (2.36) shows drawdown is lower as permeability

increases. Because permeability effect is included in PVbt,opt and interstitial velocity

measured in linear core flooding experiment, wormhole velocity does not change if

permeability changes in HWAS thus skin is not changeable.

(c) Permeability ratio:

The distribution of damage around a horizontal well is likely to be highly non-uniform

and reservoir anisotropy may lead to an elliptically shaped damage zone perpendicular to

47
the well, depending on the ratio of the vertical to horizontal permeability (Furui et al.

2002). Iani is the index of anisotropy defined by permeability ratio as:

I ani  k h / k v (2.42)

Overall skin factor (Furui et al., 2003) over the horizontal lengths is given as:

L  I ani h 
seq   ln  
 rw I ani  1
1
L
  I ani h   (2.43)
0 ln rw I ani  1  sx 
 

As anisotropy is higher, equivalent skin generates higher thus the bottomhole

pressure is calculated higher. For skin calculation, skin govern by wormhole radius thus

this effect is not included.

(d) Permeability impairment ratio (ks/k):

If the wormhole region is still inside the damaged zone, damage skin impacts the flow.

After penetrating outside the damage zone, skin factor become negative. (Fig. 27).

For rwh<rd:

k  r ( x)   r ( x) 
s ( x)  ln d   ln d  (2.44)
k d ( x)  rwh ( x)   rw 

And for rwh>rd:

 r ( x) 
s( x)   ln wh  (2.45)
 w 
r

where rwh is radius of region penetrated by wormholes at that particular point, which is to

be calculated from the wormhole model.

48
 Fig. 27 Wormhole region inside the damaged zone (Tran, 2013)

The plot (d) in Fig. 29 indicates wormhole created inside of damaged region during

early time thus skin was still huge effect thus bottomhole pressure is still high but the

values of both pressure and skin converges to lower values after wormhole penetrates

longer than damaged region and not sensitive anymore.

(e) Damage penetration radius (rd):

As same reason with (d), Eq. (2.45) shows deeper damage radius provides higher skin.

After wormhole penetrate damage radius, bottomhole pressure and skin is converged.

49
Bottom hole pressure (psi) 7000 7000

Bottom hole pressure (psi)

6500 6500
6000 6000
5500 5500
5000 5000
4500 4500
4000 0.05 0.15 0.3 4000 12 15 18
3500 3500
3000 3000
0 10 20 30 0 10 20 30
Time Time
(a) Porosity (b) Horizontal Permeability

7000 7000
Bottom hole pressure (psi)

Bottom hole pressure (psi)

6500 6500
6000 6000
5500 5500
5000 5000
4500 4500
4000 4000 0.2 0.3 0.4
3500 5 10 15 3500
3000 3000
0 10 20 30 0 10 20 30
Time Time
(c) Permeability Ratio (kh/kv) (d) Permeability Impairment Ratio (ks/k)
7000
Bottom hole pressure (psi)

6500
6000
5500
5000
4500
4000 0.5 1 1.5
3500
3000
0 10 20 30
Time
(e) Damage Penetration Radius
Fig. 28 Pressure-dependent sensitivity on reservoir and near-wellbore parameter

50
 4 4
3 3 12 15 18
2 2
1 1
Skin

Skin
0 0
-1 -1
-2 -2
0.05 0.15 0.3
-3 -3
0 10 20 30 0 10 20 30
Time Time
(a) Porosity (b) Horizontal Permeability
4 7
3 5 10 15 0.2 0.3 0.4
5
2
1 3
Skin

Skin
0 1
-1
-1
-2
-3 -3
0 10 20 30 0 10 20 30
Time Time
(c) Permeability Ratio (kh/kv) (d) Permeability Impairment Ratio (ks/k)
5
0.5 1 1.5
4
3
2
Skin

1
0
-1
-2
0 10 20 30
Time
(e) Damage Penetration Radius
Fig. 29 Skin-dependent sensitivity of reservoir and near-wellbore parameter

51
2.4.1.2. Wormhole Model Parametric Sensitivity Study

The results of pressure and skin dependent sensitivity for Furui’s wormhole model

parameters are summarize in Fig. 30 and Fig. 31. It is observed that the overall behavior

of skin-dependent sensitivity is same as pressure-dependent sensitivity. The description of

each parameter effect on pressure and skin is as follows.

(a) Optimum pore volume to breakthrough (PVbt-opt):

The relation of the average growth rate of the wormhole front, Vwh and PVbt is given by:

Vi
PVbt  (2.46)
V wh

After each data point of PVbt and Vi is obtained, the minimum point is obtained by fitting

the wormhole model to the results of laboratory tests. This lowest point defines the

optimum condition (optimal pore volume to breakthrough, PVbt-opt, and optimal interstitial

velocity, Vi-opt) . It should be noted that PVbt-opt is only an optimum in volume and smallest

volume of acid required for wormhole breakthrough.

Considering Eq.(2.46) an (2.47), smaller PVbt-opt generates higher wormhole velocity at

the tip thus skin drops faster and botomhole pressure deceases faster as well.

vi ,opt rwh
vi ,tip,opt  (2.47)
PVbt ,opt N AC Lcore

(b) Optimum interstitial velocity (Vi,opt):

Eq.(2.47) shows higher Vi-opt generates higher Vi,tip,op. According to the Eq.(2.46), Vwh is

always an increasing function of Vi. Therefore, as Vi-opt is higher, wormhole propagate

faster, reducing skin faster and botomhole pressure deceases faster as well.

52
(c) Wormhole axial-spacing coefficient (z):

In Furui wormhole model, the interstitial velocity at wormhole tip under radial flow

condition is approximated as Eq. (2.48), which assumes that in-situ injection velocity at

the tip of wormholes would be much higher than that predicted by conventional radial

flow.

q  1  1 
vi ,tip  (1   Z )   Z  
 (2.48)
L mwh  d e, wh rwh  d e , wh 

where, z is defined by the FEM simulation result as follows

0.7
h 
 Z   wh  (2.49)
 rwh 

z is set to be 0 for wormholes closely spaced in the axial direction, while z is set

to be 1 for wormholes sparsely distributed in the axial direction. As z is larger, dominant

wormhole structure is assumed and Vi,tip should be faster compared with conventional

radial flow. Therefore, wormhole propagation would be faster and pressure and skin

decreased faster.

(d) Effective wormhole diameter (de,wh):

de,wh is actual wormhole diameter inside of wormhole conduit as seen in Fig. 30. de,wh has

significant impact on Vi,tip in radial geometry as seen in Eq.(2.48). Smaller diameter gives

higher velocity. Therefore, as diameter becomes smaller and wormhole propagates faster,

so that the pressure and skin decline faster.

53
 Fig. 30 CT scan of wormhole conduit on core-flooding (Furui et al. 2012)

(e) Number of dominant wormhole per 2D or 3D plane (mwh):

Furui’s model assumes mwh is symmetry as shown in Fig. 31. It has been observed in the

literature that a few dominant wormholes grow and are spaced apart in both angular and

axial directions in a certain pattern under radial flow conditions. Eq. (2.48) implies that

Vi,tip decreases as mwh increases because of injection-fluid competition mechanism.

Fig. 31 Velocity contour plot on 2D radial flow with mwh = 6. (Furui, 2012)
54
Bottom hole pressure (psi) 7000 7000

Bottom hole pressure (psi)

6500 6500
6000 6000
5500 5500
5000 5000
4500 4500
4000 0.45 0.85 1.25 4000 0.5 1 1.5
3500 3500
3000 3000
0 10 20 30 0 10 20 30
Time Time
(a) Optimum pore volume to breakthrough (b) Optimum interstitial velocity

7000 7000
Bottom hole pressure (psi)

Bottom hole pressure (psi)

6500 6500
6000 6000
5500 5500
5000 5000
4500 4500
4000 0 0.5 1 4000 0.1 0.25 0.5
3500 3500
3000 3000
0 10 20 30 0 10 20 30
Time Time
(c) Wormhole axial-spacing coefficient (d) Effective wormhole diameter
7000
Bottom hole pressure (psi)

6500
6000
5500
5000
4500
4000 3 6 12
3500
3000
0 10 20 30
Time
(e) No. of Dominant wormholes
Fig. 32 Pressure-dependent sensitivity on wormhole model parameter

55
 4 4
3 3 0.5 1 1.5

2 2
1 1
Skin

Skin
0 0
-1 -1
-2 0.45 0.85 1.25 -2
-3 -3
0 10 20 30 0 10 20 30
Time Time
(a) Optimum pore volume to breakthrough (b) Optimum interstitial velocity

4 4
0 0.5 1 0.1 0.25 0.5
3 3
2 2
1 1
Skin

0 Skin 0
-1 -1
-2 -2
-3 -3
0 10 20 30 0 10 20 30
Time Time
(c) Wormhole axial-spacing coefficient (d) Effective wormhole diameter
4
3 6 12
3
2
1
Skin

0
-1
-2
-3
0 10 20 30
Time
(e) No. of Dominant wormholes
Fig. 33 Skin-dependent sensitivity of wormholing model parameter

56
2.4.2. History Matching Process

The result of sensitivity analysis shows that in openhole condition, damage related

parameters such as permeability impairment ratio (ks/k) and damage radius (rd) are

sensitive during early time period until wormhole penetrates the damage region. Reservoir

parameters such as porosity and permeability are sensitive during entire stimulation period.

On the other hands, among wormholing parameters used in the Furui model, wormhole

axial-spacing coefficient (z) and the effective wormhole diameter (dewh) are recognized

as most uncertain parameters.

Table 7 shows how to measure each parameter before a stimulation. As recognized,

z is inherently not given parameter in lab condition but should be obtained by history

matching process under the radial flow condition in field scale. Number of Dominant

Wormholes, mwh cannot be measured in the linear core-flooding experiment, but it does not

affect the result significantly so that it is assumed to use a suggested number.

Table 7 Measurement type of each parameter set

Well Core Fall-off
Parameter
Logging flooding test
Porosity, ϕ   -
Horizontal Permeability, kH (md)   
Permeability Impairment ratio, kS/k  - -
Permeability Ratio, kH/kV  - -
Damage Penetration Radius, rd (ft.)  - -
Optimum Pore Volume to Breakthrough, PVbt-opt -  -
Optimum interstitial velocity, Vi-opt (cm/min) -  -
Effective wormhole diameter, de,wh (in.) -  -
No. of Dominant Wormholes, mwh - Radial -
Wormhole axial-spacing coefficient, z - - -
57
2.4.3. Case Study

The limited entry case study shown in 2.2.4 is used as a history matching case study. As

described before, wormhole axial-spacing coefficient (z) is the parameter that is

optimized in the study. In this history matching, the same values of all parameters

presented by Furui at el (2010) were used. Note that this case is limited entry and the well

is hydraulically fractured before matrix acidizing, and it is assumed that damage

parameters does not affect the pressure behavior. The input parameters are shown in Table

8.

Table 8 Wormholing parameter inputs

Inputs Value
Optimum Pore Volume to Breakthrough, PVbt-opt 0.393
Optimum interstitial velocity, Vi-opt (cm/min) 1.468
Effective wormhole diameter, de,wh (in.) 0.1
No. of Dominant Wormholes, mwh 6

For this case study, integrated limited entry function introduced in 2.2.3 was used. The

result of pressure history matching is shown in Fig. 34. This match was obtained with z

set to 0.75. The correspondent result of skin history match is in Fig. 35. Due to the friction

reducer effect, actual reservoir face pressure is fluctuated but the entire match trend is

acceptable.

58
 8000 120
7000
Reservoir Face Pressure (psi) 100

Injection Rate (bpm)

6000
80
5000
4000 60
3000
40
2000
Simulated Pressure Actual Pressure
20
1000 Simulated Rate Actual Rate

0 0
0 5 10 15 20 25
Elapsed time (min)

Fig. 34 Pressure match result

0
-0.5
-1
Simulated Skin Actual Skin
-1.5
-2
Skin

-2.5
-3
-3.5
-4
-4.5
-5
0 5 10 15 20 25
Elapsed time (min)

Fig. 35 Skin match result

59
2.5. Integrated Optimization

After finishing the history match process, wormhole axial-spacing coefficient fitted by the

actual field condition is obtained. Tran (2013) studied the optimum wormhole propagation

conditions of matrix acidizing treatments in carbonate formations and compared acid

distribution and wormhole penetration along the wellbore for the increasing rate injection

to those of a constant optimum rate (determined in the parametric study), and to a constant

maximum allowable rate below the fracture gradient of the formation. He showed that the

rate needs to be increasing as more acid is injected into the formation. It actually followed

the study of Furui (2010) that showed the wormhole growth rate decreases as the acid

spends in the formation, yielding lower efficient wormhole propagation over time. If we

inject acid at a constant (even the optimum) injection rate, the interstitial velocity at the

tip of wormholes decreases as injection time increases. To keep the wormhole growth

process at the optimum conditions, the interstitial velocity needs to be maintained at or

close to its optimum as the wormhole length increases.

2.5.1. Optimum Rate Calculation

For Furui’s model, the growth rate of a wormhole extending around the wellbore is given

as,
2
 v i ,tip PVbt ,opt N AC 

   v PV 
2

 bt , opt N AC Lcore

 v i ,tip N AC    1  exp 4 
i ,tip
v wh
   (2.50)
 v i ,opt     v i ,opt rwh 
 
  

Assuming the diffusion limited reaction is the dominant process, the optimum tip velocity

is,
60
 v i ,opt rwh
v i ,tip,opt  (2.51)
PVbt ,opt N AC Lcore

The goal of the optimum acidizing treatment design is to maintain the wormhole

propagation at its optimum conditions. Especially for Furui’s model, wormhole tip is

satisfied with optimum condition all the time. Thus,

vi ,tip  vi ,itpopt (2.52)

Replacing Eq. (2.50) with Eq. (2.52) and (2.53) when the fluid loss limited wormholing,

γ = 1/ 3,
1 2
 vi ,tip,opt PVbt ,opt N AC 

  v 
2

3
 i ,tip, opt PVbt , opt N AC Lcore
 
v wh  vi ,tip,opt N AC    1  exp 4  (2.53)
 vi ,opt    vi ,opt rwh  
     

Then substituting Eq. (2.51) in to Eq. (2.53),

2
v i ,opt  rwh 3
v wh     1  exp 42 (2.54)
PVbt ,opt  Lcore 

When second term of Eq. (2.54) is simplified,

2
 v  r  3
vwh  0.96 i ,opt  wh  (2.55)
 PV 
 bt ,opt  Lcore 

From the Eq. (2.55), it is noticed that the wormhole growth rate depends on the wormhole

penetration radius and the length of the core. Assuming in a time increment Δt during

which the wormhole growth rate Vwh is constant, the wormhole penetration can be

calculated by,

n 1
rwh  rwh
n
 vwh t (2.56)

61
From the wormhole penetration we can calculate the optimal injection rate as a function

of time. Finally, a theoretical incremental rate schedule in radial flow regime is calculated

by,

vi ,tip, optL mwh

q n 1 
 1  1  (2.57)
(1   Z )  Z  
  d 
d e, whrwh  e , wh 

Practically speaking, it was noted in the literature (Glasbergen et al.,2009) that the

reservoir pressure should be fairly uniform, the total zone height should be relatively short,

and no major variations in permeability should be present in order to be a candidate for

applying the optimum injection rate. Otherwise, maximum rate should be used.

2.5.2. Case Study

Using same the limited entry case as treatment monitoring and history matching, the

optimum rate schedule is formulated based on the derivation in section 2.5.1. In this case

study, the effectiveness of calculated increasing rate schedule has been compared with the

case of both the maximum injection rate available in operated condition and the actual

treated injection rate. For the comparison, three injection rate schedule are shown in Fig.

36. It is noted that increasing rate is going to be flat because of the constraint of maximum

rate possible after 9 minutes. For the increasing rate schedule, the optimum rate is

increasing in a fast pace if one want to keep the optimum wormhole tip condition.

62
 60

Injection rate (bpm)

50

40
Actual Rate
Increasing Rate
30
Maximum Rate
20
0 5 10 15 20 25
Elapsed time (min)

Fig. 36 Rate schedule comparison

Fig. 37 shows the comparison of skin evolution based on the rate schedule. It is

observed that final skin evolutions are converged in 26 minutes because all the rate

reached the maximum rate. This means that the total acid volume used in increasing rate

schedule may be the least volume among the three rate patterns.

0
Skin (Actual)
-1 Skin (Increasing Rate)
Skin (Maximum rate)
-2
Skin

-3

-4

-5
0 5 10 15 20 25
Elapsed time (min)

Fig. 37 Skin evolution comparison

63
 However, it should be noted that the increasing rate schedule was estimated based

on the assumption that flow pattern is always radial flow and reservoir pressure is fairly

uniform, the total zone height is relatively short, and no major variations in the

permeability.

There may exist deviations in the reality from these assumptions. Since this case

was limited entry and short fractures were created before matrix acidizing, higher rate to

satisfy optimum spherical flow may be necessary.

64
 3. EVALUATION OF ACID FRACTURING THROUGH TREATMENT

MONITORING AND PRODUCTION ANALYSIS

3.1. Introduction

Unlocking tight carbonate formations for oil and gas production by multi-stage acid

fracturing in horizontal wells is relatively cost-effective as alternative to the use of

proppants. Real-time assessment of skin evolution in carbonate reservoir during matrix

acidizing treatment has been successfully applied in the industry in vertical and horizontal

wells however an extended applicability of skin evolution analysis from matrix acidizing

to acid fracturing (above closure pressure) is limited. Since fracturing can be approximated

move to a linear flow patter rather than radial flow for matrix acidizing, another method

should be developed to monitor in real-time.

In this chapter, new method for real-time monitoring of acid fracturing, the inverse

injectivity vs. superposition time function plot is presented, subject to the condition that

the treatment pressure is above closure pressure after the breakdown. A linear dual-

porosity transient slab model to simulate rate transient response in hydraulically fractured

horizontal wells normally used in production analysis has been applied for horizontal well

stimulation performance that can be implemented for field evaluation of acid fracturing

treatment. Using that model with injectivity methodology, actual growing “cross-sectional

area” by acid fracturing treatment can be monitored in real-time. The treatment is subject

to be transient flow occurred in early time thus cross-sectional area from injection flow

65
has been identified as real-time monitoring parameter, assuming bilinear flow regime

occurred during injection.

By monitoring cross-sectional area from injection flow, its real-time area growth

induced by acid fracturing is monitored during acid fracturing treatment. Even if the

condition such that pressure fluctuates around closure pressure, which makes it difficult

to identify whether fracture is closed or open, monitoring both existent skin evolution

under matrix acidizing and acid fracturing effect allows us to observe the effect of acid

stimulation comprehensively. Another benefit of selecting cross-sectional area as a

monitoring parameter is that one can examine the stimulation effectiveness by conducting

production data analysis after production starts. Linear flow diagnostic approach with

Rate-Transient Analysis provides cross-sectional area flowing from matrix, which is

compared with the area induced by acid fracturing in real-time. The treatment efficiency

provides engineers with additional information as to whether the designed acid fracturing

was performed appropriately under the in-situ closure stress field.

3.2. Treatment Monitoring in Acid Fracturing in Horizontal Well

3.2.1. Dual-Porosity Transient Slab Model

The purpose of this section is to introduce analyzing a fracturing treatment in horizontal

well with same monitoring method approach shown in matrix acidizing. Using the

transient dual porosity solution which is most commonly used in the literature in rate

transient analysis for a linear model originally presented by El-Banbi and extended by

66
Bello and Wattenbarger as shown in Fig. 38, necessary flow regime occurring during

fracturing treatment was incorporated to existing method.

However skin over time during acid fracturing treatment is difficult to be generated

conceptually because fracturing skin is given by pseudo radial flow. Our treatment is

subject to be transient flow occurred in early time thus cross-sectional area from injection

flow has been identified as real-time monitoring parameter. Monitoring cross-sectional

area from injection flow, its real-time growth is monitored during acid fracturing treatment.

If we assume rectangular shape of fracture, approximate fracture half-length is monitored

in real-time. Using a field example, we demonstrate the augmented method be consistent

and useful.

matrix
fractures
ye

xe
h

Fig. 38 Schematic of slab matrix linear model well (Bello, 2009)

The transient dual porosity solutions presented here for a linear model was seen in

the literature (El-Banbi, 1998). Bello and Wattenbarger extended the model towards an

actual multi-fractured horizontal well, which has the advantage of being simpler than other

horizontal well models and allows relaxing five flow-regime. The diffusivity equations

67
for the matrix and fracture are solved in Laplace space. Although mathematical details of

the linear dual porosity model (slab matrix) are found in the literature (Bello, 2009), the

constant rate inner boundary and closed outer boundary reservoir (slab matrix) solution in

Laplace space is given in Eq. (3.1).

2 1  e 2 uf ( u ) y De 
p wDL    (3.1)
u uf (u ) 1  e  2 uf ( u ) y De


where, dual porosity parameters for the slab matrix case,

 Ac 31   u
f (u )    1    tanh (3.2)
3u  Ac

12 k m
 Ac  Ac (3.3)
L2 k f

( ct ) f
 (3.4)
( c t ) f  ( ct ) m

It is defined that f(u) is a Laplace space function used in transient dual porosity

model. λAC is the Warren and Root dimensionless interporosity flow parameter and

controls how fast the fluid drains from the matrix to the fractures and ω is dimensionless

storativity ratio and controls how much fluid is initially in the fracture system. In our

analytical model, flow-regime assumes to be bilinear flow which has simultaneous

transient flow in the fracture system and matrix.

The following dimensionless variables for gas reservoir are defined to convert

dimensionless pressure and time into pressure and time in real space. Note that Ac is

defined as cross-sectional area to flow. It is noted that unit of t is hours when 0.000264.

68
 0.000264k f t
t DAc  (3.5)
(  ct ) f  m Ac

kf Ac [m( p i )  m( p wf )]
m DL  (3.6)
1422q g T

ye
y De  (3.7)
Ac

Since injection pressure is normally more than 3,000 psi, the pseudo pressure in Eq. (3.6)

can be converted into single pressure as follows:

kf Ac  p i  p wf   2 p i 
p DL    (3.8)
1422q g T  i zi 

3.2.2. Bilinear Analytical Model

In a fracturing treatment with injected liquid, bilinear flow caused by simultaneous

transient flow into the fracture system and matrix is assumed.

The derivation of constant rate solution for bilinear flow regime are shown as

below. Starting with Eq. (3.1) again,

2 1  e 2 uf ( u ) y De 
p wDL    (3.9)
u uf (u ) 1  e 2 uf ( u ) y De


This can be shown to be the same as:

pwDL 
2
u uf (u )

Coth uf (u ) y De  (3.10)

where,

e2x  1
Coth x   (3.11)
e2x  1

69
With approximately for x>3, Coth ( x)  1 . Eq. (3.10) becomes,

2
p wDL  (3.12)
u uf (u )

For slab matrix:

 Ac 31   u
f (u )    1    tanh (3.13)
3u  Ac

 Ac 31   u
Assuming   1    tanh then, 1    1
3u  Ac

3 1    u
Approximately for x >3, tanh  x   1 . Since 31   u  3, tanh 1.
 Ac Ac

Therefore substituting those assumptions into Eq. (3.13), the simplified form is,

 Ac
f (u )  (3.14)
3u

Then substituting Eq. (3.14) into Eq. (3.12):

2 (2 )(3) 0.25

p wDL  
  u 1.25  Ac
0.25 0.25
(3.15)
uu 0.5  Ac 
 3u 

Mathematically,

1.251
 1  t t 0.25
 1  1.25    (3.16)
 u  1.25 0.9064

1  (3.17)
 1    1
u 

Inverting from Laplace to real space:

(2 )(3) 0.25 0.25 9.12305 0.25

pwDL  t  t DAc (3.18)
0.9064 Ac  Ac 0.25
0.25 DAc

70
Eq. (3.18) is dimensionless constant rate solution for bilinear flow model.

3.2.3. Approach

Derived Eq. (3.18) is applied to the treatment monitoring methodology by Zhu and Hill

(1996, 1998). One big difference with production analysis, however, is that fracture keeps

growing, thus the area of fracture, Ac is not a constant as depicted in Fig. 39.

Fig. 39 Stimulation schematic (Plan view)

Using Eq. (3.5) and (3.8), Eq. (3.18) is re-written in linear relationship, noting that

quadratic root of time is used as superposition time function because of bilinear flow.

Applying the dimensionless pressure definition to convert to real pressure:

kf Ac  p wf  p i   2 p i  9.12305 0.25
   t DAc (3.19)
 i zi   Ac
0.25
1422q T

71
where, kf is adopted as effective fracture permeability (Kazemi, 1969) in fracture spacing,

L.

w k f  (l  w) k m wk f
kf   (3.20)
L L

Since in the direction perpendicular to the transverse fracture, the flow from high

conductive fracture region is restricted to low matrix permeability region, weighted

average permeability is reasonable in the dual porosity model.

Rearranging Eq. (3.19),

p wf  pi 1422T   i z i  9.12305 0.25

   t DAc
k f Acw  2 pi
(3.21)
  Ac
0.25
q

Substituting dimensionless time of Eq. (3.5) into Eq. (3.21):

0.25
p wf  p i 1422T   i z i  9.12305  0.000264k f 
   0.25   t 0.25 (3.22)
q k f Ac  2 p i   Ac  (  c t ) f  m Ac 

where,
12 k m
 Ac  Ac (3.23)
L2 k f

Rearranging the Eq. (3.22) for Ac,

0.25
pwf  pi 1 1422T  i zi  9.12305  0.000264k f 
     t 0.25
 (  ct ) f  m 
0.75 0.25
q Ac kf  2 pi   12 km  (3.24)
 A 
 L2 k c 
 f 

p wf  p i 1 1422T   i z i  9.12305 L  0.000264 

0.25
(3.25)
     t 0.25
k f  2 p i  12k m   (  c t ) f  m
0.25
q Ac 

Now superposition time function is written as usual,

72
 q j  q j 1
t  t j 1 
N
tsup  
0.25
N (3.26)
j 1 qN

Slope m is defined by Eq. (3.26).

0.25
1422T   i z i  9.12305 L  0.000264 
m     (3.27)
 12k m   (  c t ) f m
0.25
k f  2 pi 

Rearranging again using superposition time, the following equation is obtained:

pwf  pi 1
 mtsup (3.28)
q Ac

Finally cross sectional area to flow is obtained with reforming Eq. (3.28),
mtsup
Ac (t ) 
pwf  pi (3.29)
q

Assuming rectangular shape of fracture, cross sectional area to flow in bilinear flow
regime is given as follows:
Ac (t )  4wh  4hx f n f (3.30)

Since area of injection flow inside of fracture is much narrower than fracture-half length,
first term of Eq. (3.30) should be minimized. Thus, Ac is approximated into:
Ac  4wh  4hx f n f  4hx f n f (3.31)

Thus, fracture half-length is obtained in each time step as:

Ac (t )
xf  (3.32)
4hn f

Summarizing the procedure to calculate real-time evolution of cross-sectional area from

injection flow as a function of injection time by continuous monitoring of injection rates

and pressures;

73
1. Calculate slope, m, using Eq. (3.27) for fractured horizontal wells in gas reservoir. The

slope, m, is constant and depends on the reservoir parameters and dimension such as

formation volume factor, permeability, and viscosity of the reservoir fluid.

2. From the measured pressure and injection rate at a desired time interval, calculate

superposition time using Eq. (3.26) for horizontal wells.

3. Calculate cross-sectional area from injection flow using Eq. (3.29) at each time step.

3.2.4. Bilinear Flow Validation by Numerical Simulation

In this section, the assumed bilinear environment during injection treatment in

hydraulically fractured reservoir is validated by the commercial numerical simulator. Fig.

40 shows a finite-difference grid system employed in this work, where three transverse

fractures with fracture half-length 250 ft. are set up but only toe stage is opened and

simulated for the injection.

Fig. 40 Logarithmically spaced LGR grid system

74
 To capture transient behavior occurred during injection, near-induced fracture

region by acid fracturing is expressed by the logarithmically spaced local grid refinement

(LGR). In this LGR model, the global grid which contains hydraulic fractures is divided

into 21 different widths according to a logarithmically distribution in x direction. The table

of input summary used in this work is shown below.

Table 9 Geometry of the simulation grid data

Parameters Input Values
Dimensions (nx,ny,nz) = (7,3,1)
Grid block size (dx,dy,dz) = (400, 500, 300) ft
Matrix permeability kx,ky,kz = 0.0044 md
Fracture permeability kx,ky,kz = 2000 darcy
Fracture width 0.001 ft.

Table 10 Reservoir & simulation data

Reservoir Properties Well and Stimulation Data
Input
Parameters Input Values Parameters
Values
Initial Reservoir Pressure 9837 psi Wellbore Radius 3.307 in.
Formation Thickness 300 ft Tubing Diameter 2.992 in.
Total Porosity 0.03 Stimulation Interval 400 ft
Total Compressibility 5.78E-05 1/psi Reservoir top 20,000 ft
Injection fluid
Initial Sw (matrix) 0.1 1 cp
viscosity
Initial Sw (Fracture) 1 Injection rate 15 bpm

Based on Brooks and Corey model, two relative permeability curves are set up for

matrix and fracture respectively as shown in Fig. 41.

75
 Gas-Water Relative Permeability Gas-Water Relative Permeability
1 (Matrix) 1 (Fracture)
0.8 0.8

0.6 0.6 krg

Kr
Kr

0.4 0.4 krw

krg
0.2 krw 0.2

0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1
Sg Sg

Fig. 41 Relative permeability for matrix and fracture

Fig. 42 is the plan view of pressure profiles around the wellbore from early time.

It is seen that after the injected flow reached fracture-half length (250 ft.), fluid was

expanded to x axis. Thus the assumption made in analytical model is validated.

Fig. 42 Plan view of pressure shows bilinear flow

76
 As the simulation result, Fig. 43 shows that normalized rate vs. time plot in log-

log scale. The diagnostic plot identifies ¼ slope which means that the flow regime during

injection was bilinear. Thus, linear flow inside of the fracture and linear flow from the

fracture to matrix occurs simultaneously. The deviation of early time is considered as

infinite acting radial flow from the center of wellbore. After reaching the boundary, the

flow was converged to bilinear flow.

10000

¼ slope
(Pwf-pi)/qw (psi/bpm)

1000

100

10
0.1 1 10 100 1000
Time (min)

Fig. 43 Normalized rate plot during the injection

3.2.5. Case Study

3.2.5.1. Treatment Description

The acid stimulation in toe stage is analyzed and monitored. The input data for the

reservoir, the wellbore for the analysis is summarized in Table 11.

77
 Table 11 Input data for fracture monitoring

Reservoir Properties Well and Stimulation Data

Input
Parameters Input Values Parameters Values
Initial Reservoir Pressure 9837 psi Wellbore Radius 3.307 in.
Formation Volume Factor 0.003 ft3/scf Tubing Diameter 2.992 in.
Total Porosity 0.03 Stimulation Interval 1,357 ft.
Total Compressibility 5.78E-05 1/psi Vertical Depth 19,880 ft.
Formation Thickness 250 ft. Measured Depth 23,343 ft.
Reservoir Fluid Viscosity 0.0315 cp
Reservoir Temperature 280 F (*) Effective permeability:
Z factor 1.36
matrix permeability 0.0044 md wk f 2000[md  ft]
kf    1.5 [md]
fracture permeability (*) 2000 D L 1357[ft]
fracture width 0.001 ft.

3.2.5.2. Monitoring Result

Fluid injection schedule from surface to measured depth at stage 1 is in Table 12.

Table 12 Treatment schedule of stage 1

Injection Schedule
Stage Name Viscosity, cp Volume (gal)
1 Gelled Acid 39 18492
2 X-linked acid 30 31701
3 Non-X-linked frac fluid 56 52834
4 X-linked acid 30 36984
5 Non-X-linked frac fluid 56 52834
6 Gelled acid 39 18492
7 Water Flush 1 5283

78
 The surface tubing pressure and injection rate were measured on-site during the

injection. Based on hydrostatic and friction pressure calculation by frictional pressure drop

analysis made beforehand, the bottomhole pressure was generated in Fig. 44. As seen,

each circle number corresponds to the fluid number in Table 12 and each dashed line also

corresponds to each fluid arrival at bottomhole.

It is noted that calculated reservoir face pressure fluctuates between the formation

breakdown pressure and the closure pressure during the stimulation. It is also noted that

the reservoir face pressure often decreased below the closure pressure but stayed around

the closure pressure. The decision as to whether acid fracturing or matrix acidizing was

happening strongly the calculated bottomhole pressure, breakdown, closure pressure and

acid etched conductivity. Thus both matrix acidizing skin and cross-sectional growth

induced by acid fracturing treatment were calculated in the monitoring process. Based on

Fig. 44 data, rate normalized plot in log-log scale is plotted to diagnose related flow regime.

Fig. 45 shows that after formation breakdown, bilinear flow occurred during acid injection.

Fracture linear flow is identified when frac fluid used, which means fracture extension.

79
 Surface pressure Reservoir Face Pressure Slurry Rate
18000 100
Breakdown Pressure 90
16000
Closure Pressure 80
14000
frac sleeve open 70

Slurry Rate, bpm

Pressure, psi

12000 60
10000 50
2 3 5
8000 40
Completion 30
fluid
6000 6
4 20
1 7
4000 10
2000 0
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170
Time, min
Fig. 44 Stage 1 treatment pressure & rate

Fig. 45 Diagnostic plot to identify flow regimes.

80
Using treatment data in Fig. 44, the skin evolution when matrix acidizing is supposed

and cross-sectional area growth when acid fracturing is supposed were computed as

shown in Fig. 46.

Fig. 46 Evolution of skin and cross-sectional area growth

From Fig. 46, it is observed that consistent cross-sectional area growth during the

treatment. Strictly speaking from the diagnostic plot in Fig. 45, the area growth during

xanthan frac fluid injection is not valid because the flow regime indicated fracture linear

flow, ½ slope and further fracture extension. However, after fracture linear flow, it is noted

that bilinear flow, ¼ slope followed the fracture linear flow and the cross-sectional area

growth was monitored. Final cross-sectional area growth and skin before water flush are

81
100,000 ft2 and -1.8 respectively. The equivalent fracture half-length is 140 ft, if fracture

is assumed to be rectangular shape and the fracture height is 250 ft. uniformly.

3.3. Production Data Analysis after Acid Stimulation

When production data is available, it can be used to conduct a rate transient analysis (RTA).

In our purposes, this analysis can be used to calculate cross-sectional area flowing from

fractured wall, Ac and so that we can compare the calculated area between the stimulation

and production so that the efficiency can be shown quantitatively in a field scale.

3.3.1. RTA Procedure for Gas Well

3.3.1.1. Generate the Specialized Plot

For dry gas reservoirs, the specific gravity of the gas sample obtained at the primary-

separator conditions equals the value of a sample from the reservoir. However, for gas-

condensate reservoir fluids, the specific gravities of gas samples taken from the reservoir

and at the surface differ. By using recombination calculation, dry gas production was

corrected to account for the condensate in the gas phase (Lee and Wattenbarger, 1996)

Qg ,wet  Qg ,dry  133.316 oQo / M o  (3.33)

From PVT report, oil in molecular weight is given. Thus, API and oil specific gravity is

calculated.

5954
API   8.811 (3.34)
Mo

82
 141.5 (3.35)
o 
131.5  API

Recombined wet gas volume is used as production rate in this analysis. On the other hand,

calculated bottomhole pressure in commercial RTA software is converted to real gas

pseudo pressure and finally, Δm(p)/qg is obtained the square root of time plot is generated.
p
p
m( p)  2  dp (3.36)
po
z

The production data, petrophysics, PVT, wellbore schematics, deviation survey etc.

was loaded into commercial RTA software (fekete HARMONY) to obtain the bottomhole

flowing pressure From the PVT report, reservoir is identified as gas condensate reservoir

in this field case and recombined gas rate was calculated. Calculated bottomhole pressure

was converted to pseudo pressure m(p) and finally Δm(p)/qg vs. √𝑡 plot was generated.

For rate fluctuations due to operating changes, superposition time should be used.

m pi   m p wf  ~
m 4 t b (3.37)
qg

3.3.1.2. Diagnose the Specialized Plot

Read straight line slope from Δm(p)/qg vs. √𝑡 plot and identify the time to the end of linear

flow. Finally using Eq. (3.38) and (3.39), stimulated rock area and fracture half-length (xf)

is estimated from this plot and compared with real-time monitoring results.

803.2T 1
k m Acm 
ct  f  m ~
m4 (3.38)

Assuming stimulated rock volume to be uniform, xf is calculated as follows:

83
 Acm
xf  (3.39)
4h n f

3.3.2. Case Study

3.3.2.1. Production History

The well in monitoring case study is used. The summary of production history is shown

in Fig. 47.

Surface Pressure Casing pressure

3000 Oil rate water rate 350 1400
gas rate
2500 300 1200

250 1000

Oil, water rate (boepd)

Pressre (psi)

Gas rate (mcfd)

2000
200 800
1500
150 600
1000
100 400
500 50 200

0 0 0
0 20 40 60 80 100 120
Days

Fig. 47 Production history

As seen, in 113 days production, there was no shut-in during the period. Condensate, water

as well as gas rate and surface pressure was recorded in surface. Choke size when well

produced was 5 mm as constant during those periods. Along with methodology developed

before, recombined gas rate and bottomhole pressure are calculated and shown in Fig. 48.
84
 Calculated Bottomhole Pressure Surface Pressure
5000 Recombined Gas Volume gas rate 1400

1200
4000
Pressre (psi)

1000

Gas rate (mcf/d)

3000 800

600
2000

400
1000
200

0 0
0 20 40 60 80 100 120
Days
Fig. 48 Recombined gas rate and calculated bottomhole pressure

3.3.2.2. Production Analysis Result

Based on recombined gas rate and calculated bottomhole pressure, pressure normalized

rate vs. material balance time in log-log scale is generated and plotted in Fig. 49.

85
 ½ slope

Fig. 49 Normalized rate vs. material balance time

Due to the low quality of data, it’s not clear trend but half slope which indicates matrix

linear flow is shown. Very early time indicates clean-up and the skin effects and it may

mask early time behavior. As next step, square root time plot is generated in Fig. 50.

Corresponding
Straight line

Fig. 50 Square root of time plot

86
The slop m4 was read as 700000 (700 × 103 ). Now that the slope is known, √k m Ac

should be calculated using Eq. (3.38). The constant values of Table 13 are used for the

calculation.

Table 13 Wellbore & stimulation data

Wellbore and stimulation assumptoin
Course length # of frac Effective spacing Thickness porosity visicosity ct
ft ft ft cp psi^-1
3832 1 1357 250 0.03 0.0315 5.78E-05

Assuming the rectangular shape of fracture and 0.0044 md for matrix permeability that

was obtained from the well log interpretation, cross-sectional area to flow in formation

linear flow and uniform fracture half-length of each stage is calculated as below.

Table 14 Calculated fracture geometry

Summary of caluculation
Ack^1/2 Acm xf
md^0.5.ft^2 ft^2 ft
3581 53,990 54

Finally this area is compared with monitoring result shown in 3.2.5 and the summary is

obtained in Table 15. It shows actual contributed area from flow is approximately half

of stimulated area. This information may be helpful for optimization of next stimulation

job.

87
 Table 15 Comparison summary
Monitoring analysis
Production analysis Efficiency

Stimulated area (ft2) Flowing area (ft2) (Production/Stimulation)

100,000 53,990 54%

88
 4. CONCLUSION

In this thesis, integrated approach to evaluate both matrix acidizing and acid fracturing

(“acid stimulation”) has been implemented. Since, in many cases, once the acid reaches

the formation, injection pressure cannot be kept above the formation closure stress, the

fracture eventually closes, and the treatment becomes matrix acidizing rather than acid

fracturing. Appropriate design and diagnostic method should be necessary to evaluate the

treatment.

For matrix acidizing, using the existent numerical simulator and skin monitoring

program, consolidating them into one package has been done. The package allows us to

conduct acidizing design, monitoring the pressure and skin evolution and history match

them and update uncertain parameters to use future treatment. Optimum rate schedule is

also calculated based on the integrated approach.

For acid fracturing, the success of the treatment depends on many factors as to

whether enough conductivity is secured, selected treatment works well in in-situ under

specific geologic environment. Thus, observation and evaluation of past practice is

critically important and inevitable step to develop further optimal stimulation procedures.

We developed the methodology to conduct performance evaluation of acid fracturing

treatment using treatment and production records.

The suggested integrated approach provides engineers with additional information

as to whether the designed acid stimulation was performed appropriately under the in-situ

89
closure stress field. It is eventually helpful to look back past practice and enhance future

stimulation.

90
 REFERENCES

Agnia, A., Alkouh, A., and Wattenbarger, R. A. 2012. Bias in Rate-Transient Analysis
Methods: Shale Gas Wells. Paper SPE 159710 presented at the SPE Annual
Technical Conference and Exhibition, San Antonio, Texas, 8-10 October.
Al Ahmadi, H. A., Almarzooq, A. M., and Wattenbarger, R. A. 2010. Application of
Linear Flow Analysis to Shale Gas Wells - Field Cases. Paper SPE 130370
presented at the SPE Unconventional Gas Conference, Pittsburgh, Pennsylvania,
23-25 February.
Bazin, B. 2001. From Matrix Acidizing to Acid Fracturing: A Laboratory Evaluation of
Acid/Rock Interactions. SPE Production & Facilities 16 (1):22-29, SPE-66566-PA.
Bello, R.O. 2009. Rate Transient Analysis in Shale Gas Reservoirs with Transient Linear
Behavior. PhD dissertation, Texas A&M University, College Station, Texas.
Bello, R. O. and Wattenbarger, R. A. 2010. Modelling and Analysis of Shale Gas
Production With a Skin Effect. Journal of Canadian Petroleum Technology, 49(12),
37-48.
Ben-Naceur, K. and Economides, M. J. 1989. Design and evaluation of acid fracturing
treatments. Paper SPE 18978presented at the SPE Low Permeability Reservoirs
Symposium, Denver, Colorado, 6-8 March.
Buijse, M.A. and Glasbergen, G. 2005. A Semiempirical Model to Calculate Wormhole
Growth in Carbonate Acidizing. Paper SPE 96892 presented at the SPE Annual
Technical Conference and Exhibition, Dallas, 9-12 October.
Cohen, C.E., Ding, D., Quintard, M., Bazin, B. 2008. From porescale to wellbore scale:
Impact of geometry on wormhole growth in carbonate acidization. Chem. Eng. Sci.
63 (12): 3088-3099.
Dong, K., Jin, X., Zhu, D., Hill, A.D. 2014. The Effect of Core Dimensions on the Optimal
Acid Flux in Carbonate Acidizing. Paper SPE 168146 presented at the SPE
International Symposium and Exhibition on Formation Damage Control, Lafayette,
Louisiana, 26-28 February.
Economides, M.J., and Frick, T.P. 1994. Optimization of horizontal well matrix
stimulation treatments. SPE Production & Facilities, 9(02), 93-99. SPE-22334-PA.
Economides, M.J., Hill, A.D., and Ehlig-Economides, C. 1994. PetroleumProduction
Systems. Upper Saddle River, New Jersey: Prentice-Hall, Inc.

91
El-Banbi,A.H. 1998. Analysis of Tight Gas Wells. PhD Dissertation, Texas A&M
University, College Station, Texas.
Elbel, J. L. 1994. Field Evaluation of Acid-Fracturing Treatments With Geometry
Simulation, Buildup, and Production Data. SPE Production & Facilities, 9(01), 43-
46, SPE-19773-PA.

Fredd, C.N. and Miller, M.J. 2000. Validation of Carbonate Matrix Stimulation Models.
Paper SPE 58713 presented at the SPE International Symposium on Formation
Damage Control, Lafayette, Louisiana, 23-24 February.

Furui, K., Burton, R.C., Burkhead, D.W., Abdelmalek, N.A., Hill, A.D., Zhu, D., and
Nozaki, M. 2010. A Comprehensive Model of High-Rate Matrix Acid Stimulation
for Long Horizontal Wells in Carbonate Reservoirs: Paper SPE 134265 presented
at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19-22
September.

Goode, P. A. 1987. Pressure drawdown and buildup analysis of horizontal wells in

anisotropic media. SPE Formation Evaluation, 2(04), 683-697. SPE-14250-PA.

Hill, A.D. and Zhu, D. 1996. Real-Time Monitoring of Matrix Acidizing Including the
Effects of Diverting Agents. SPEPF11 (2): 95–101.

Izgec, O., Key, R., Zhu, D., and Hill, A.D. 2008. An Integrated Theoretical and
Experimental Study on the Effects of Multiscale Heterogeneities in Matrix
Acidizing of Carbonates. Paper SPE 115143 presented at the SPE Annual
Technical Conference and Exhibition, Denver, 21-24 September.

Jackson, A.M., Azizi, B.A., Kofoed, C.W., Shuchart, C.E., Keller, S.R., Sau,R., Grubert,
M.A., Phi, M.V. Completion and Stimulation Methodology for Long Horizontal
Wells in Lower Permeability Carbonate Reservoirs. Paper SPE 161527 presented
at Abu Dhabi International Petroleum Conference and Exhibition, Abu Dhabi,
UAE, 11-14 November

Jawad, M.S. Al. 2014. Modeling of Acid Fracturing in Carbonate Reservoirs. Master of
Science. Texas A&M University. College Station, Texas.

Kalfayan, L.J. 2007. Fracture Acidizing: History, Present State, and Future. Paper SPE
106371 presented at the SPE Hydraulic Fracturing Technology Conference,
College Station, Texas, 29-31 January.

Kalia, N., and Balakotaiah, V. 2009. Effect of medium heterogeneities on reactive

dissolution of carbonates. Chemical Engineering Science, 64(2), 376-390.

Kazemi, H. 1969. Pressure transient analysis of naturally fractured reservoirs with uniform
fracture distribution. SPE J. 9(04): 451-462.. SPE-2156-A.
92
Kent, A.W., Burkhead, D.W., Burton, R.C., Furui,K., Actis,S.C., Bjornen,K.,
Constantine,J.J., Gilbert,W.W., Hodge,R.M., Ledlow, L.B., Nozaki, M., Vasshus,
A., and Zhang, T. 2014. Intelligent Completion Inside Uncemented Liner For
Selective High-Rate Carbonate Matrix Acidizing. SPE Drilling & Completion,
29(02), 165-181. SPE-166209-PA.
Lee, J., Wattenbarger, R.A. 1996. Gas Reservoir Engineering. SPE Textbook series Vol.
5.
McDuff, D., Jackson, S., Schuchart, C., and Postl, D. 2010. Understanding Wormholes in
Carbonates: Unprecedented Experimental Scale and 3D Visualization. J Pet
Technol 62 (10): 78-81. SPE-129329-MS.
Mishra, V., Zhu, D., Hill, A.D., and Furui, K. 2007. An Acid-Placement Model for Long
Horizontal Wells in Carbonate Reservoirs. Paper SPE 107780 presented at the
European Formation Damage Conference, Scheveningen, The Netherlands, 30
May-1 June.
Mostofizadeh, B. and Economides, M.J. 1994. Optimum Injection Rate From Radial
Acidizing Experiments. Paper SPE 28547 presented at the SPE Annual Technical
Conference and Exhibition, New Orleans, 25-28 September.
Paccaloni, G., Tambini, M., and Galoppini, M. 1988. Key Factors for Enhanced Results
of Matrix Stimulation Treatments. Paper SPE 17154 presented at the SPE
Formation Damage Control Symposium, Bakersfield, California.
Pandya, N.D. 2012. A Systematic Study of Matrix Acidizing Treatments Using Skin
Monitoring Method. Master of Science. Texas A&M University. College Station,
Texas.
Pournik, M., and Nasr-El-Din, H. A. 2010. Influence of Acid-Fracture Fluid Properties on
Acid-Etched Surfaces and Resulting Fracture Conductivity. Paper SPE 128070
presented at the SPE International Symposium and Exhibition on Formation
Damage Control, Lafayette, Louisiana, 10-12 February.
Prouvost, L.P. and Economides, M.J. 1989. Applications of Real-Time Matrix-Acidizing
Evaluation Method. SPE Production Engineering 4 (4): 401-407. SPE-17156-PA.
Samandarli, O., McDonald, B., Barzola, G., Murray, M., and Richmond, P. 2014.
Understanding Shale Performance: Performance Analysis Workflow With
Analytical Models in Eagle Ford Shale Play. Paper SPE 169004 presented at the
SPE Unconventional Resources Conference, The Woodlands, Texas, 1-3 April.
Sasongko, H., Zhu, D., and Hill, A. D. 2011. Simulation of Acid Jetting Treatments in
Long Horizontal Wells. Paper SPE 144200 at the SPE European Formation
Damage Conference, Noordwijk, The Netherlands, 7-10 June.

93
Tran, H.T. 2013. Modeling and Optimization of Matrix Acidizing in Horizontal Wells in
Carbonate Reservoirs. Master of Science. Texas A&M University. College Station,
Texas
Warren, J.E. and Root, P.J. 1963. The Behavior of Naturally Fractured Reservoirs. SPE J.
3 (3): 245-255. SPE-426-PA.
Zhu, D. and Hill, A.D. 1998. Field Results Demonstrate Enhanced Matrix Acidizing
through Real-Time Monitoring. SPE Production & Operations 13 (4): 279-284.
SPE-52400-PA
Zhu, D., Hill, A. D., and Looney, M. D. 1999. Evaluation of acid treatments in horizontal
wells. Paper SPE 56782 presented at the SPE Annual Technical Conference and
Exhibition, Houston, Texas, 3-6 October.
Zhu, D., Hill, A. D., and da Motta, E. P. 1998. On-site evaluation of acidizing treatment
of a gas reservoir. Paper SPE 39421 presented at the SPE Formation Damage
Control Conference, Lafayette, Louisiana, 18-19 February.

94