University of Southern California Dissertations and Theses

Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

(USC Thesis Other)

Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

Waterflood Induced Formation Particle Transport and Evolution of
Thief Zones in Unconsolidated Geologic Layers

By

Qianru Qi

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)

December 2020

Copyright 2020 Qianru Qi
ii
Dedication

To my beloved mother, Bei Qi, whose sacrifices made it possible for me to
succeed in life.

iii
Acknowledgements
I am deeply grateful to my advisor Dr. Iraj Ershaghi for his continuous inspiration,
encouragement, and support throughout my years at USC. His generosity for sharing his time
and insights have made my graduate study and this dissertation possible. His phenomenal
mentoring has made me a different individual, as he consistently helps me discover myself
academically and professionally. He had done so from the very beginning when I was first
exposed to the subject of petroleum engineering until today as I complete my Ph.D. dissertation.
I give him my most profound respect and appreciation.
I also sincerely thank Dr. Donald Paul, Dr. Jincai Chang, Dr. Richard Robert, and Dr.
Katherine Shing, who served on my Ph.D. committee. They kindly offered their time, advice, and
suggestions to assist me in the completion of the study from a multitude of perspectives. They
have been considerate of accommodating my qualifying exam and final defense processes. I want
to particularly thank the great professors at USC PTE, especially Dr. Donald Hill, Dr. Donald
Gautier, and Mr. Brad Pierce, for enriching me with their knowledge, advice, and experiences.
I am thankful to Ms. Ebonie Hawthorne, Mr. Andy Chen, Ms. Cory Reano, Mr. Sorin
Marghitoiu, Ms. Idania Takimoto, Ms. Juli Legat, Ms. Nicole Kern, Ms. Aimee Barnard, the DEN
staff, and many other people who offered their care, assistance, and support during my joyful and
memorable time at USC.
I especially thank Miss Sophie Pepin, Miss Katarina Bethel, Miss Leilani Rebolledo, and Mr.
Abdulla Al-Jazzaf, for working together with me on various exciting research topics. My sincere
thanks also go to all the students I served as a teaching assistant (TA), as they also provided me
a context to refresh my understanding of petroleum reservoirs.
I thank my elders Mr. Xiaorong Zhang, Dr. Hanming Cui, Dr. Zhengming Yang, and my dear
friends, Ms. Qiuyan Zhang, Ms. Huafeng Zhu, Miss Na Zhao, and Miss Yuanchun Li, who provided
their encouragement along my Ph.D. journey.
iv
Finally, I appreciate the support received from Computer Modelling Group (CMG) for the
academic license provided to conduct reservoir simulation studies in completing part of the
research project, and California Resources Corporation (CRC), who provided data to complete
this research project and dissertation.

v
Table of Contents
Dedication ................................................................................................................................... ii
Acknowledgements ................................................................................................................... iii
List of Tables ............................................................................................................................ viii
List of Figures ............................................................................................................................. x
Abstract .................................................................................................................................... xiv
Introduction ................................................................................................................................. 1
Excess Water Production Issues and Sand Production Issues .......................... 1
1.1: Environmental aspects of excess water production ........................................................ 1
1.2: Sand production issues and remediation measures ....................................................... 6
1.3: Solutions to mitigate sand productions and associated costs ........................................ 6
1.3.1: Gravel pack ........................................................................................................... 7
1.3.2: Frac pack .............................................................................................................. 9
1.3.3: Plastic consolidation ............................................................................................. 9
1.3.4: High-energy resin placement .............................................................................. 10
1.3.5: Resin-coated gravel pack ................................................................................... 10
1.3.6: Slotted liners or unit filter .................................................................................... 11
1.3.7: Rate control ........................................................................................................ 12
1.3.8: Other measures .................................................................................................. 13
Causes of Sand Production .................................................................................. 14
2.1: The behaviors of unconsolidated formations ................................................................ 14
2.1.1: Dynamic reservoir properties as a function of pore pressure,
permeability, skin, and compressibility ......................................................................... 14
2.1.2: Subsidence ......................................................................................................... 15
2.1.3: Grain movement ................................................................................................. 16
2.1.4: Grain breakage ................................................................................................... 17
2.1.5: Thief zone ........................................................................................................... 17
2.2: Rock characteristics of unconsolidated formations ....................................................... 28
2.2.1: Formation strength, mineral composition of cementation and
formation fines .............................................................................................................. 28
2.2.2: Examples of unconsolidated sands with grain size distributions ........................ 31
vi
2.3: Current understandings on sand particle movements and injectivity issues
for injection wells ................................................................................................................. 33
Injectivity Maintenance and Review of Step Rate Test ...................................... 39
3.1: Injectivity monitoring and permeability variations ......................................................... 39
3.2: Fracturing signals from SRT vs. grain movement in multi-layered sands ..................... 40
3.2.1: A critical review of finding formation breakdown pressure in multi-
layered reservoirs ......................................................................................................... 42
3.2.2: Particle migration ................................................................................................ 44
Modeling of Sand Grain Movement ...................................................................... 45
4.1: Analytics ....................................................................................................................... 47
4.1.1: Assumptions ....................................................................................................... 47
4.1.2: Pressure distribution and apparent fluid flow velocity ......................................... 47
4.1.3: Critical shear velocity for particle migration ........................................................ 48
4.1.4: Grain size distribution and grain population density ........................................... 51
4.1.5: Pore size variation .............................................................................................. 53
4.1.6: Particle movement and migration ....................................................................... 55
4.1.7: Update in particle population density, porosity, and permeability ....................... 56
4.2: Stochastic modeling using the Markov Chain Probability Transition
process ................................................................................................................................ 58
4.3: Inverse modeling .......................................................................................................... 61
4.4: Algorithm ....................................................................................................................... 62
Validation approach ............................................................................................... 65
Results .................................................................................................................... 71
Conclusions and Discussions .............................................................................. 87
7.1: Conclusions .................................................................................................................. 87
7.2: Discussions ................................................................................................................... 89
7.2.1: Completion method and formation damage affecting sand
migration in unconsolidated formations ........................................................................ 89
7.2.2: Low pore volume injected in developing a thief zone ......................................... 99
7.3: Highlights .................................................................................................................... 101
vii
Suggested Future Work ...................................................................................... 102
8.1: Extensions of modeling in multi-phase flow with different grain wettability,
and consideration of formation strength ............................................................................ 102
8.2: Other aspects of work to be clarified related to the sand particle movement
process .............................................................................................................................. 106
Nomenclature .......................................................................................................................... 109
References .............................................................................................................................. 112

viii
List of Tables

Table 1.1: Estimated Cost For Typically Handling One Barrel Of Water. ([1].
Courtesy SPE) .............................................................................................................................. 3
Table 1.2: Reasons And Remediation Measures Associated With Types Of Above
Diagnostic Plots. ([1] [9]. Courtesy SPE.) ..................................................................................... 5
Table 1.3: Summary of A Sand Control Program. ([15] Courtesy SPE.) .................................... 12
Table 2.1: Water Profile Survey Results For Vertical Injector From Kuwait. ([40]
Courtesy SPE.) ........................................................................................................................... 20
Table 2.2: Water Profile Survey Results For The Horizontal Injector From Kuwait.
([40] Courtesy SPE.) ................................................................................................................... 22
Table 2.3: Injectivity Increase Associated With Well Profile In Figure 2.1. ([41],
Courtesy SPE.) ........................................................................................................................... 24
Table 2.4: X-Ray Analysis Of GOM Unconsolidated Formation Fines. ([49].
Courtesy SPE.) ........................................................................................................................... 30
Table 3.1: Data From A Multi-Profile Test Before Polymer Treatment Associated
With Figure 2.6. ([42]. Courtesy SPE.) ....................................................................................... 42
Table 4.1: Water Viscosity And Density. ..................................................................................... 49
Table 5.1: Reservoir Properties For Base Case. ........................................................................ 66
Table 5.2: Sensitivity Parameterizations. .................................................................................... 66
Table 5.3: Hand Extracted And Extrapolated Grain Size Distribution Values From
Figure 5.2. ([16] Courtesy SPE.) ................................................................................................. 69
Table 6.1: Permeability Increase Ratio For 10 Ft. Thick Sand With 800 Md Initial
Permeability. ............................................................................................................................... 78
ix
Table 6.2: Permeability Increase Ratio For 10 Ft. Thick Sand With 1500 Md Initial
Permeability. ............................................................................................................................... 78
Table 6.3: Permeability Increase Ratio For 10 Ft. Thick Sand With 2000 Md Initial
Permeability. ............................................................................................................................... 79
Table 6.4: Permeability Increase Ratio For 10 Ft. Thick Sand With 3000 Md Initial
Permeability. ............................................................................................................................... 79
Table 6.5: Permeability Increase Ratio For 30 Ft. Thick Sand With 800 Md Initial
Permeability. ............................................................................................................................... 79
Table 6.6: Permeability Increase Ratio For 30 Ft. Thick Sand With 1500 Md Initial
Permeability. ............................................................................................................................... 79
Table 6.7: Permeability Increase Ratio For 30 Ft. Thick Sand With 2000 Md Initial
Permeability. ............................................................................................................................... 80
Table 6.8: Permeability Increase Ratio For 30 Ft. Thick Sand With 3000 Md Initial
Permeability. ............................................................................................................................... 80
Table 6.9: Permeability Increase Ratio For 50 Ft. Thick Sand With 800 Md Initial
Permeability. ............................................................................................................................... 80
Table 6.10: Permeability Increase Ratio For 50 Ft. Thick Sand With 1500 Md Initial
Permeability. ............................................................................................................................... 80
Table 6.11: Permeability Increase Ratio For 50 Ft. Thick Sand With 2000 Md Initial
Permeability. ............................................................................................................................... 81
Table 6.12: Permeability Increase Ratio For 50 Ft. Thick Sand With 3000 Md Initial
Permeability. ............................................................................................................................... 81
Table 6.13: Critical Injection Gradient To Incur Sand Migration And Erosion In
Open-Hole Injection Wells .......................................................................................................... 86
Table 7.1: Thickness of Thief Zone Compared to the Formation Thickness ............................ 100
Table 7.2: Assumed Properties For Computing The Case Shown In Figure 2.9. ..................... 100
x
List of Figures

Figure 1.1: Summarized Diagnostic Plots in Identifying Water Production Reasons.
([9]. Courtesy SPE.) ...................................................................................................................... 4
Figure 2.1: Hall Plot For The Vertical Injector From Kuwait. ([40] Courtesy SPE.) ..................... 20
Figure 2.2: Water Profile Survey For The Vertical Injector From Kuwait. ([40]
Courtesy SPE.) ........................................................................................................................... 21
Figure 2.3: Hall Plot Of The Horizontal Injector From Kuwait. ([40] Courtesy SPE.)
.................................................................................................................................................... 22
Figure 2.4: Water Profile Survey For The Horizontal Injector From Kuwait. ([40]
Courtesy SPE.) ........................................................................................................................... 23
Figure 2.5: Effect Of A Shock Treatment On Profile And Injectivity Of A Well. ([41],
Courtesy SPE.) ........................................................................................................................... 24
Figure 2.6: Injection Profile Behavior At Various Rates And Pressures, Lower
Jones Injector A, Huntington Beach Offshore Field. ([42]. Courtesy SPE.) ................................ 25
Figure 2.7: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well A. ([47]
Courtesy SPE.) ........................................................................................................................... 25
Figure 2.8: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well B. ([47]
Courtesy SPE.) ........................................................................................................................... 26
Figure 2.9: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well C. ([47]
Courtesy SPE.) ........................................................................................................................... 26
Figure 2.10: Schematic Of A Reservoir With Horizontal Fracture (A) Vs. With A
Horizontal Washed-Out Zone In Unconsolidated Formation (B). ([47] Courtesy
SPE.) .......................................................................................................................................... 27
xi
Figure 2.11: Dynamic Combination Modulus Curve As A Function Of Depth. ([22]
Courtesy SPE.) ........................................................................................................................... 29
Figure 2.12: Examples Of Fine Particle Located On The Surface Of Larger
Formation Sand Grains. ([49] Courtesy SPE.) ............................................................................ 30
Figure 2.13: Uniform Sized, Unconsolidated California Sand. ([16] Courtesy SPE.)
.................................................................................................................................................... 31
Figure 2.14: Nonuniform Sized, Unconsolidated California Sand. ([16] Courtesy
SPE.) .......................................................................................................................................... 32
Figure 2.15: Typical Grain Size Distribution Of GOM Unconsolidated Sand ([51]). ................... 32
Figure 3.1: Channeling Shown From A Temperature Log In An Injection Well. ([68].
Courtesy SPE.) ........................................................................................................................... 40
Figure 3.2: Illustration Of Stabilized Pressure – Plot From SRT. ([68]. Courtesy
SPE.) .......................................................................................................................................... 41
Figure 4.1: Apparent Flow Velocity At Various Completion Conditions. ..................................... 48
Figure 4.2: Settling Velocity 𝑤𝑠 Calculated With Van Rijn (1993) Correlations. ......................... 50
Figure 4.3: Number Of Realizations Generated In Stochastic Process. ..................................... 54
Figure 4.4: Workflow And Inverse Modeling Of The Particle Migration Algorithm. ..................... 64
Figure 5.1: Schematic Of The Slab Reservoir Assumed In Simulation. ..................................... 66
Figure 5.2: Cumulative Grain Size Distribution Of Sampled California
Unconsolidated Sand. ([16] Courtesy SPE.) ............................................................................... 67
Figure 5.3: Incremental Grain Size Distribution Of Sampled California
Unconsolidated Sand: Data Extracted From Figure 5.2 by Schwartz (1969). ([16]
Courtesy SPE.) ........................................................................................................................... 67
Figure 5.4: Initial Grain Size Distribution Extrapolated upon Figure 5.2. .................................... 68
Figure 5.5: The Assumed Initial Pore Size Distribution. ............................................................. 68
xii
Figure 5.6: Critical Shear Velocity At Various Depths And Formation Temperatures.
.................................................................................................................................................... 70
Figure 6.1: Initial Status, As Well As All Realizations For The First, Second, And
The Third Steps. ......................................................................................................................... 74
Figure 6.2: Averaged Porosity Evolution For The Base Case Described In Table
5.1. .............................................................................................................................................. 75
Figure 6.3: Averaged Permeability Evolution For The Base Case Described In ........................ 75
Figure 6.4: Harmonically Averaged Permeability Increase Vs. Injection Duration
For The Base Case Described In Table 5.1. .............................................................................. 76
Figure 6.5: Porosity Profile Evolutions For The Base Case Mentioned In Table 5.1
Except For A Change In Injector-Producer Distance. ................................................................. 83
Figure 6.6: Permeability Profile Evolutions For The Base Case Mentioned In Table
5.1 Except For A Change In Injector-Producer Distance ............................................................ 85
Figure 7.1: Fluid Drag Forces Vs. Completion Conditions. ([15] Courtesy SPE.) ....................... 90
Figure 7.2: Evolution Of Grain Size Distribution On Initial Status, As Well As All
Realizations For The First, Second, And The Third Steps. ........................................................ 94
Figure 7.3: Porosity Profile Evolutions For The Base Case Mentioned In Table 5.1
Except For A Change In Injector-Producer Distance And A Change From Open-
Hole To Perforated Well. ............................................................................................................ 96
Figure 7.4: Permeability Profile Evolutions For The Base Case Mentioned In Table
5.1 Except For A Change In Injector-Producer Distance And A Change From
Open-Hole To Perforated Well. .................................................................................................. 98
Figure 7.5: A Comparison Among The Cases With Various Injector-Producer
Distances In A Perforated Injection Well .................................................................................... 99
Figure 8.1: Observations On Wettability And Multi-Phase Affecting Particle
Movement. ................................................................................................................................ 105
xiii
Figure 8.2: Formation Strength Affected By Saturation In Multi-Phase Condition.
([15]. Courtesy SPE.) ................................................................................................................ 106
xiv
Abstract
Recovery of hydrocarbons from unconsolidated sands under waterflooding frequently results
in excess water and sand production. Such operational practices may cause reservoir subsidence
and damage to production facilities with associated environmental and safety hazards. During the
lifespan of a reservoir, these two issues can severely and negatively affect operations by
overloading subsurface and surface facilities and pipelines, thereby decreasing hydrocarbon
productivity and substantially increasing operational cost. These often happen during water
injection into unconsolidated formations, where unbonded grains in some layers can migrate
along with fluid flow under injection pressure. Left behind are unswept zones and layers, as more
in-situ grains are dislodged and carried away from near-wellbore area, causing erosion and
development of high permeability zones. Associated changes in reservoir properties, such as
porosity and permeability, are not negligible and must be modeled. Modeling of sand migration,
as a predictive tool for identification of thief zones in a timely manner, should be incorporated as
a critical part of reservoir simulation studies.
In this study, a new algorithm modeling sand particle migration is constructed assuming
single-phase flow in unconsolidated sands composed of unbonded, spherical grains following
typical grain size distributions. The formation pore size distribution assumed is also in the range
reported for these types of formations. The algorithm developed is capable of providing statuses
of grain size distribution at any stage of waterflooding. As an indication of a washed-out zone or
a thief zone being developed, the computed sequential grain size distributions model the
migration of various-sized particles, which are further used in computing the dynamics of reservoir
permeability.
Given adequate pore sizes allowing grain migration, our method calculates the movement of
grains of various sizes when local apparent fluid velocity exceeds critical shear velocity of grains.
xv
In this way, the physics of grain mobilization is better honored as one can anticipate that grains
are subject to gravitational forces and viscous forces in a porous system.
As possibilities exist for both pore throat bridging and pore throat enlargement due to the
local particle movement, variations in pore size distribution are modeled stochastically by
incorporating the Markov Chain probability transition process. During the computation, it is noted
that detailed measurements of grain size distribution, especially for the smaller-sized end, are
needed. Monitoring of formation temperature, which affects fluid viscosity and particle mobilization,
is also important. The algorithm developed requires minimal measurements of reservoir
characteristics, where only general reservoir properties, pore size distribution, and grain size
distribution are required. It also allows validation and re-selections of particle loss factor and
transition matrices through comparison to periodically measured profile surveys.
The study includes a computational example using a typical grain size distribution for
unconsolidated sands presented in the literature with the assumption of other reservoir properties.
The results show that at each location and time, several realizations of grain size distribution are
mapped. Permeability is calculated based on counting the particle population density for each
realization of grain size distribution. This way, with a profile survey, realizations that match the
real permeability increase shown on the survey can be found, as the possible reflections of
reservoir dynamics. In this study, the computed permeability increase is tabulated for scenarios
considering formation depth, initial permeability, injection pressure gradient, and layer thickness.
Moreover, factors including pressure drawdown per unit length, well completion type, and partial
formation damage are emphasized.
With improvement in accuracy of representing formation strength under varying water
saturation and grain wettability conditions, the sand migration model presented here can guide
the prediction and timely remediation measures from an early stage of reservoir development for
recovery under water injection. More importantly, it must be included into reservoir simulations to
history match the waterflooding data in unconsolidated sands.
1
Introduction
In all subsurface injection of fluids using water or gas, when dealing with unconsolidated
formations, sand production can become an operational issue. One side effect of sand grain
movement can be a gradual increase in permeability of some individual layers, causing thief
zones and also sand production. The dissertation focuses on addressing the causes of sand
production and includes new mathematical modeling describing sand migration in unconsolidated
geologic formations. The primary focus is on water injection.
Chapter 1 introduces the issues of excess water production and associated sand production.
Besides its economic burdens on oil and gas operations, excess water production brings short-
term and long-term harm to environmental safety. For unconsolidated formations, sand production
can be a major issue in damaging hydrocarbon recovery process and causing potential formation
subsidence. Some remediation measures are systematically reviewed, with a summary of
application conditions and key experiences for each of the actions included.
Chapter 2 describes the causes of sand production from perspectives of formation
characteristics and flow conditions. Current techniques used to classify those conditions were
also examined. Several laboratory tests, mathematical models, and interpretations of the
processes are reviewed. The correct way of identifying thief zones is discussed with examples.
An analytical solution for detecting a thief zone during its gradual development using a series of
fall-off tests are also discussed.
Chapter 3 briefly reviewed the significance of maintaining injectivity in waterflooding
operations. In addition to profile survey, temperature logging is another method in monitoring the
changes in permeability around injection wells. As concerns for fracturing of unconsolidated
formations exist, a case study is reviewed where pressure response was interpreted as a
fracturing signal, while the test was done on multi-layered unconsolidated sands. It brings up a
very important point that, for a well with commingled flow, using Step Rate Test data to find
2
fracturing gradient is questionable. Fine migration could be a possible reason for slope decrease
on the pressure-rate plot.
In Chapter 4, a new algorithm is proposed to model sand grain movement by using critical
velocity criterion for sand grain dislodge. The new algorithm computes changes in grain
population density following the stochastic Markov Chain probability transition process, which
allows various possibilities to variation in pore size distribution. The multiple realizations
generated can be incorporated into inverse modeling by comparison with profile surveys. If
available, profile survey may also suggest finer discretization on grain size distribution curve, and
re-selection of probability transition matrices taken in the Markov chain process. As indicated by
changes in grain size distribution at each location and time, resultant updates in porosity and
permeability are calculated. The algorithm is based on the physical nature of particle migration
process considering gravitational and viscous forces, while traditional solutions are primarily
based on modeling of local concentration gradient.
In Chapters 5 and 6, validation of the algorithm with a computational example and its results
are presented. The computation process highlights the benefit of a more accurate measurement
of grain size distribution especially on the finer-size end. Results are shown for changed grain
size distribution at each location and time. Along the process, an increase in permeability due to
voidage created by particle loss can be calculated, for a sand layer with specific thickness, depth,
initial permeability, subject to a specific injection pressure gradient. The results show that an
unconsolidated formation with a higher initial permeability, or subject to a higher injection gradient
is prone to more severe particle migration and sand erosion. Thickness may not affect the sand
migration much, if assuming the sand face is clean around the open-hole injection well. It is also
noted that, distance between injection well and producing well is a key factor as it controls the
pressure change per distance to trigger grain mobilization. Sensitivity of the injector-producer
distance is tested, in addition to considerations of injection gradients, formation depth, initial
3
permeability, and sand thickness. ‘Safe’ injection gradients are calculated with the assumptions
of open hole injector around which the sand face is clean.
Chapter 7 includes two key aspects of the discussions, in addition to the conclusions reached
with the computations. First, for perforated injection wells and the injection wells around which
the sand face is partially damaged, sand migration is much more severe than in clean sand face
around the open-hole injectors. Second, in reality, thief zone(s) could be very thin to jeopardize
the profile survey. As a result, on a field scale, we may observe a thief zone being created with
much less fluid injected than seen from lab measurements or numerical simulations. Both aspects
of the discussions are illustrated with examples.
Chapter 8 outlines the extension of the research work to model particle migration using the
proposed algorithm. Additional work suggested includes modeling under multi-phase condition
considering the grain particle wettability, and modeling with different conditions such as
completion type, initial skin around wellbores, grain angularities, and clay type.
1
Excess Water Production
Issues and Sand Production Issues
1.1: Environmental aspects of excess water production
It is estimated that the average ratio of water-to-oil production volume is 3: 1 worldwide. The
treatment, reinjection, as well as disposal of produced water can cost $40 – 50 billion USD
annually [1] [2]. Excess water production issues become more severe when oilfields approach
maturity. Water, as a by-product of oil and gas operations, includes formation water and injected
water. Formation water can come from bottom or adjacent aquifers supporting reservoirs, or
aquifers connected to reservoirs through faults, natural fractures, and natural or artificially-induced
high permeability zones [3]. More significantly, the huge volume of produced water comes from
injected water used for waterflooding, steam flooding, or any other water-based displacement
methods. Connections between producing wells and injection wells via faults, natural fractures,
natural or artificially created high permeability channels become the source of the issue [1].
Produced water is classified as 'acceptable water' and 'excess water' with respect to technical
capabilities and economic tolerance. 'Acceptable water' refers to produced water that can be
economically justified by recovered hydrocarbon, which also cannot be terminated without
shutting off oil production. Also, 'acceptable' water production can be minimized by enhanced
operational practice [1].
'Excess water' refers to unwanted water produced for various reasons. Because of its
chemical composition and disposal issues, excess water production poses negative effects on
onshore and offshore environments. Produced water contains a substantial number of
compounds that vary with geological and geochemical conditions of reservoirs and operations on
each well. The compounds typically include:
2
1. dispersed oil droplets, water-soluble low molecular organic acids, alkylphenols, biocides,
aromatic molecules, and other organic materials;
2. heavy metals, naturally occurring radioactive material (NORM), solid particles, sulfur,
sulfide, and other inorganic materials; and,
3. chemicals added, including corrosion inhibitors, scale inhibitors, clay stabilizers, and
emulsion breakers, etc. [4].
These compounds not only affect the quality of produced water for subsurface re-injection
and ocean disposal, but also leave residual harm to ecology due to bioaccumulation and toxicity,
which especially comes from chemicals of aromatic hydrocarbons and heavy metals [4] [5].
Produced water-to-oil ratio for a well could vary from zero to upwards of 99% as a reservoir
gets more depleted [6]. For natural gas operations, produced water accounts for 80% of total
wastes [7]. Produced water is also considered as waste of the largest volume for oil and gas
operation on most offshore platforms [8].
Excess water production also burdens oil and gas operators as it causes corrosion to
downhole and surface equipment and pipelines. Other side effects may include scale, wax,
asphaltene deposition, and gas hydrates. Besides indirect damages to reservoirs, just handling
produced water itself, the cost can range from $0.05 to $0.50 USD per barrel. Sometimes, the
cost related to water can be high as $4 USD for each barrel of oil produced in fields with a water
cut being higher than 80% [1]. The cost significantly increases as the water cut of produced fluid
increases along the life of a reservoir. Bailey et al. (2000) [1] statistically showed the breakdown
cost of handling water in typical situations on lifting, separation, pumping, reinjecting, and
treatment, as shown in Table 1.1. Note that the price of oil in that year was $27.39 USD (nominal).

3
Table 1.1: Estimated Cost For Typically Handling One Barrel Of Water. ([1]. Courtesy SPE)

20,000
B/D
50,000
B/D
100,000
B/D
200,000
B/D
Average
Lifting
Capex/Opex

$0.044
5.28%
$0.044
7.95%
$0.044
9.29%
$0.044
10.25%
$0.044
7.69%
Utilities
$0.050
6.38%
$0.054
9.62%
$0.054
11.24%
$0.054
12.40%
$0.054
9.30%

Separation
Capex/Opex

$0.087
10.36%
$0.046
8.27%
$0.035
7.24%
$0.030
6.82%
$0.049
8.55%
Utilities
$0.002
0.30%
$0.003
0.45%
$0.003
0.52%
$0.003
0.58%
$0.003
0.43%
Chemical
$0.034
4.09%
$0.034
6.16%
$0.034
7.20%
$0.034
7.94%
$0.034
5.95%

De-oiling
Capex/Opex
$0.147
17.56%
$0.073
12.99%
$0.056
11.64%
$0.046
10.58%
$0.081
13.92%
Chemicals
$0.040
4.81%
$0.041
7.25%
$0.041
8.47%
$0.041
9.34%
$0.041
7.00%

Filtering
Capex/Opex

$0.147
17.47%
$0.068
12.18%
$0.047
9.85%
$0.030
6.87%
$0.073
12.63%
Utilities
$0.012
1.48%
$0.010
1.79%
$0.010
2.09%
$0.010
2.31%
$0.011
1.84%

Pumping
Capex/Opex

$0.207
24.66%
$0.122
21.89%
$0.091
19.06%
$0.079
18.15%
$0.125
21.61%
Utilities
$0.033
3.99%
$0.034
6.01%
$0.034
7.03%
$0.034
7.75%
$0.034
5.81%
Injecting Capex/Opex
$0.030
3.62%
$0.030
5.45%
$0.030
6.37%
$0.030
7.02%
$0.030
5.27%

Total cost/bbl

$0.842
100%

$0.559
100%
$0.478
100%
$0.434
100%
$0.578
100%
Total chemicals
$0.074
8.90%
$0.075
13.41%
$0.075
15.67%
$0.075
17.28%
$0.075
12.96%
Total utilities
$0.102
12.16%
$0.010
17.87%
$0.100
20.88%
$0.100
23.03%
$0.101
17.38%
Total wells
$0.074
8.89%
$0.075
13.40%
$0.075
15.66%
$0.075
17.27%
$0.075
12.95%

Surface
facilities
$0.589
70.05%
$0.309
55.33%
$0.227
47.80%
$0.184
42.41%
$0.328
56.71%

As such, best practices include avoiding excess water and other water control measures that
can significantly benefit oil and gas operations economically and environmentally. It is of key
4
importance that the application of water control measures requires clear identification of the
source(s) of produced water. Tools and methods used on water problem analysis include density
logging, temperature logging, and spinner logging, running profile surveys, multi-phase fluid
logging, installation and usage of downhole sensors, reservoir characterizations, application of
inter-well tracers, and pressure transient tests [1]. Guidance has been offered on classifying
water sources based upon the evolvement of water-cut with time [9]. The diagnostic curves are
summarized in Figure 1.1, where the recorded curves of water-oil-ratio (WOR) with respect to
time provide key information of various situations. Possible reasons for each of the situations are
explained in the first two columns of Table 1.2.

Figure 1.1: Summarized Diagnostic Plots In Identifying Water Production Reasons. ([9].
Courtesy SPE.)

It is recognized that the optimized practice for water control is to separate water as early as
possible. Techniques for eliminating water can include downhole separation, injection, chemical
and mechanical shutoff [1], or keeping production rate below the economic limit of water-oil-ratio.
Typical methods also include re-completion, cement plugging and bridging, drilling lateral drain
holes for vertical wells or sidetrack for horizontal wells, and so on.
Measures on treatments inside reservoirs, such as using gels and foams, can be effective in
control water influx. However, their economic and material costs can be substantial, considering
1
2
3
4
WOR
Time
5
a large formation volume to be treated. Those well-practiced methods on water control are
classified for various types of water problems, as co-tabulated in the last column of Table 1.2.
Another more recent method for water control is to artificially cause formation damage and
permeability reduction by injecting low salinity water, through mobilization of formation fines [3].

Table 1.2: Reasons And Remediation Measures Associated With Types Of Above Diagnostic
Plots. ([1] [9]. Courtesy SPE.)
Reasons for Water Production
Curve
Type
Methods
Leaks from casing, tubing or packers 1, 2
Squeezing shutoff fluids (gel-based
or resin-based [10]),
shutoff using plugs, cement, packer
and patches

Leaks behind casing 1
High-strength squeeze cement,
placement of resin-based fluids or
lower-strength gel-based fluids in the
annulus

Coning or cresting,
during or cusping for horizontal wells
3, 4
Cement plug, bridge plug, gel
injection or drill lateral drain holes
near formation top for vertical wells
Sidetrack for horizontal wells

Channeling from injector 4 Gel injection, Water shutoff

High permeability streak 4

In-fill well drilling, drilling lateral drain
holes for vertical well, near-wellbore
shutoff for horizontal well

Completion near water zone 1, 2
Cement plug, bridge plug, gel
injection or drill lateral drain holes
near formation top

Fracturing propagated to water, faults and
fractures connecting injection well and
producing well, or connecting to water
zone
1, 2 Gel injection, water shutoff
Watered-out zone 3, 4
Carefully drill lateral drain holes
around undrained layer, mechanical
shutoff of injector or producer

6
1.2: Sand production issues and remediation measures
Sand production issues have been observed from field operations and documented in the
literature for many decades [11] [12]. It not only concerns oil and gas operators but also affects
disposal and geothermal projects [13] [14].
In petroleum operations, sand production is an issue for both producing wells and injection
wells. It reduces well productivities and accessible hydrocarbon reserves, which may also cause
formation collapse [15]. For injection wells, it is an issue for achieving the goals of pressure
maintenance in sand layers, which further jeopardizes profile surveys, vertical sweep efficiency,
and the ultimate recovery factors of entire fields. Besides its negative effects on the fulfillment of
a field’s potential, sand production also brings other issues, such as:
1. erosion of downhole and surface equipment, valves, pipes, artificial lifts, etc.;
2. frequent workovers and associated shutoffs of production and injection; and,
3. unexpected leaks and failures of equipment, implying safety and environmental hazards.
All of the problems mentioned above are associated with a substantial economic cost. For
instance, the total cost of sand control actions taken on the wells of the Gulf Coast owned by
Amoco alone was reported as $10 - 15 million USD annually [15]. As a reference, the crude oil
price was $14.40 USD (nominal) in that year.

1.3: Solutions to mitigate sand productions and associated
costs
Many types of actions have been taken for sand control purposes. They includes gravel
packing, frac packing, utilization of slotted liners, plastic consolidation, resin-coated gravel
packing, rate control, or other passive methods. Passive methods include periodically and
routinely workovers and equipment washing when sand accumulates inside wellbores and near-
wellbore regions of formations. They also include selective completion, meaning not producing
7
from the sand layers that indicate a high potential for sand production. For each case, the
selection of sand control measures depends on specific reservoir condition, mechanical condition,
logistic condition, and cost, etc. Details on selection among those measures and the key
experiences are discussed as following.

1.3.1: Gravel pack
Gravel packing is by far the most successful sand control method [15] in either a perforated
or open-hole well. Successful application of gravel packs demands a series of actions to be done
properly. These include good drilling and perforations, careful measurements of formation grain
size distribution, wise selection of gravel sizes, tight packing, as well as cleaning and maintenance
of pipes and equipment.
Gravel pack works well for thick intervals, which may or may not have high clay content.
Generally, its cost is the lowest compared to other methods. Placement of gravel pack works for
both open hole and perforated casing. It is harder to place gravel pack tightly in a well with
perforated casing, and removal of which can be expensive if erosion happens. When the objective
layer is deep, the selected size of gravel packed screen or liner become small, accordingly, the
outcome of sand control become less than ideal.
Studies have been done to address the key points of successful placement of gravel packed
screens in detail. Schwartz (1969) [16] proposed a standardized workflow to design gravel packs
in response to sand production issues in wells operated under high rates in waterflooding
operations. The study demonstrated that, in addition to the proper design of gravel size and
uniformity of those gravels, a successful design of gravel pack also requires accurate analyses of
formation sand size distribution and fluid velocity through slots. Further improvements can be
made by optimization of gravel pack thickness and appropriate selection of completion fluids.
A program on sand control [15] showed that, to achieve the best sand control objective,
primary steps include selection of drilling fluid, treated with salts and clay stabilizers to avoid clay
8
dispersion. Lowering mud density and controlling tripping rates to minimize mud invasion are
required as well. This is important, especially for perforated wells, where flow area is limited so
that apparent flow velocity becomes substantially higher. If some of the perforations are plugged
by the damaging particles caused by drilling, flow velocity will be even higher to trigger sand
mobilization and erosion. In the same principle, adding chemicals also plays a role, in order to
remove oxygen and avoid corrosion.
To keep near-wellbore formation clean to the extent possible, actions such as continuous
filtering of particles bigger than 2 microns for non-solid brine systems is another step. Performing
perforations in an underbalanced condition, as well as conducting perforation washing and back
surging, are emphasized. Wellbore cleanout fluid and gravel pack placement fluids are suggested
to be brine to ensure a clean wellbore for the best gravel packing operation [17]. It was also noted
that perforation of no less than 4 or 6 shots per foot is preferred for optimum sand control results
using gravel-pack filter/screens. Sufficient perforations help avoid too-high flow velocities through
limited flow area [15] [17].
Jennings Jr. (1997) [18] suggested the selection of bigger sized gravels and more spherical
grains (less angular) in designing gravel pack based on the experimental study. The study showed
that using too-small gravels hurts well productivity, while selecting too-large gravels could not
meet the expectation on sand control [15].
Besides the sizing of gravels, tight packing of gravel-packed screens is another key factor for
sand control and prevention of corrosive wear. This is not only noticed from field practice but also
indirectly verified from lab experiments. Measured pore throat diameters from mercury injection
through uniform-sized gravels in rhombic packing are smaller than the theoretical pore size
calculated, indicating certain possibility of poorly packed gravels, which hurts the result of sand
control [15].
9
As shown from the statistical data collected from operators, the average amount of gravels
per foot is substantially higher in the successful gravel packing jobs than in the failed cases, for
either open holes or perforated holes [17].

1.3.2: Frac pack
Bruno et al. (2001) [19] modeled the failure process of unconsolidated granular materials with
a coupled particle and fluid flow approach, which can be useful in guiding the design of frac packs
and prediction of their performance. Hwang and Sharma (2014) [20] experimentally investigated
the fluid velocity during filtration via frac packs, filled with proppants of various sizes. When bigger
proppant is selected, the extent of particle removal significantly decreases as fluid flow rate
increases. However, when smaller proppant is used in frac packs, their particle removal capability
does not change as much. Note that both scenarios require an accurate understanding of water
quality.

1.3.3: Plastic consolidation
Plastic consolidation, or plastic injection, had attracted some attention but later became less
used. Through the injection of plastic at a temperature of 225 ℉ or above, loose formation sand
particles, within a radius of about 3 ft. around the wellbore, get consolidated after plastic curing.
As a result, sand control is achieved at the price of losing some formation permeability. In
comparison to gravel packing, plastic consolidation can work for deeper zones as it does not leave
obstruction inside the narrow wellbore.
Even coverage of plastic over entire formation interval is essential to control sand production,
which implies a need for less permeability variation. Therefore, the best formation candidates are
no more than 12 [17] or 15 ft. thick [15]. Shaly or silty intervals add permeability variation.
Moreover, shale particles have bigger surface area, which hinders the consolidation of injected
plastic. As such, plastic consolidation does not work well for formations with high clay contents.
10
Repeated consolidation jobs are not suggested, as continuous reduction of formation
permeability offsets its sand control benefits. With all that restrictions, a once-and-forever job of
plastic consolidation on carefully selected candidates needs to be guaranteed. Also, those cases
with its plastic consolidation conducted prior to initial production show a higher successful rate
[17].

1.3.4: High-energy resin placement
Resin injection at a high rate and pressure has been used in practice to overcome uneven
placement of resin over target sand intervals. It delivers a more even coverage of target interval
at such a rapid rate, before formation has a chance to fail. There are several ways to achieve a
high-rate injection, such as high-pressure surging, overbalanced perforating, and proppant gas
fracturing.
For either typical plastic consolidation, or high-rate resin injection, their applications are
limited due to permitting issues as chemicals are involved in resin formula.

1.3.5: Resin-coated gravel pack
Resin-coated gravel packing refers to applying dry gravels coated with about 5% (in weight)
resin around all perforations. It is followed by steam injection, so that consolidated gravels stay
inside the perforations when temperature in the formation decreases. The parts left inside the
casing can be drilled out of the wellbore so that only perforations are filled with resin-coated
gravels.
The method works for deep formation where the wellbore is narrow, which also leaves the
wellbore obstruction-free, similar as the plastic consolidation from this perspective. It appears that
Resin-coated gravel packing overcomes the disadvantages and limitations of both gravel packing
and plastic consolidation; however, success of applying resin-coated gravel pack is not easy to
reach.
11
It requires all perforations to be fully filled, which becomes much more difficult when interval
thickness increases. It is also hard to be applied in slanted holes. Another factor controlling sand
control result is resin curing inside perforations. Among the successful cases, not all resin-coated
gravel packings last long, which limits its application.

1.3.6: Slotted liners or unit filter
Unit filters include plastic liner/screens, wire-wrapped screens, or filter sections filled with
beads/agates connected to a slotted-liner [21]. A Survey [17] showed that, due to the fewer slot
plugging by formation clays which left the pathways for fluid flow, unit filters are reported to be
more successful for cleaner and more homogeneous sands with fewer shale breakers around
perforations. It also showed that fluid velocity is higher in filter-installed sections with a perforation
density less than four shots per foot. Applications in such situations lead to more sand production
and a higher possibility for filter cut out.
Wire-wrapped screens are less used compared to gravel packs. They suits for situations
when formation sloughs shortly after initial production of the well [15]. With such reservoir
condition, screens need to be packed completely to avoid erosion under high flow velocity and
sand influx. It holds true, especially for perforated wells where flow is restricted, and apparent flow
velocity is even much higher to mobilize sand particles.
Table 1.3 summarizes the applied completion methods in dealing with sand production issues
using several different methods. The success rate was improved to 92 percent of all 182 applied
jobs [15].

12
Table 1.3: Summary of A Sand Control Program. ([15] Courtesy SPE.)
Type Job Total
Initial
Success
Initial
Failures
Subsequent
Failures
A gravel-packed screen in
perforated casing
164 157 7 5
A gravel-packed screen in
open hole
5 5 0 0
Plastic consolidation 8 6 2 0
Triple-wrapped screen, not
gravel-packed
5 5 0 0
Total 182 173 9 5

1.3.7: Rate control
The goal of rate control is to find out a maximum flowing rate at which the amount of sand
production is acceptable, without involving other sand control methods. Calculation of the
maximum allowable rate needs to be done repetitively, as the water cut and reservoir pressure,
changes during production.
The critical pressure drawdown, ∆𝑃
!
, at a formation of depth 𝐷, can be estimated from the
following equation:
[∆𝑃
!
]
"
= [∆𝑃
!
]
#

(%
"
)
#
(%
"
)
$
(Eq. 1.1, [22])
where subscript 𝑇 refers to the zone tested for reference and 𝐷 refers to the zone of interest,
(𝐸
%
)
&
as shear modulus, can be estimated from density log and acoustic velocity log.
Assuming critical pressure drawdown is proportional to shear modulus, considered as an
index for formation strength, critical rate, (𝑞
!
)
&
, is,
(𝑞
!
)
"
=
'.)*×,)
'(
-. [∆1
)
]
#
34(56
*
+
*
,
8 9 : ).';)
(Eq. 1.2, [22])
Details are shown in Section 2.2.1. Parameters ℎ,𝑘,𝐵,𝜇,𝑟
-
,𝑟
.
and 𝑠 refer to formation
thickness, permeability, fluid formation volume factor, fluid viscosity, draining radius, wellbore
radius, and skin factor, respectively.
13
In practice, instead of having a single value for critical rate or pressure drawdown, a safety
margin needs to be considered. Flowing at a rate below the lower curve of the risk region are
considered ‘safe’, while above the upper curve of the risk region are considered as a high
possibility to ‘fail’ in sand control. The test is conducted with a series of increased flow rates. 'Safe'
rate refers to observation of no apparent sand accumulation, and 'fail' rate refers to the point
where excessive sand production enters the wellbore.
The method is less preferred for the economic reason that the maximum allowable rate is far
less than the reservoir potential. It is chosen mostly when applications of other sand control
methods cannot be economically justified by hydrocarbon production, or, when mechanical
condition of the well does not allow other methods to be applied.

1.3.8: Other measures
Many other kinds of remediation measures used in practice have been documented in
response to sand production issue. When multi-layered unconsolidated sands are under
waterflooding, sub-zoning ensures that each layer receives sufficient water injection for both
pressure maintenance purposes and prevention of potential subsidence, as mentioned in
Robertson et al. (1987) [23]. In deepwater reservoirs with severe sand production issues, where
the high cost, high risk, and logistic difficulty simultaneously challenge the development of the
reservoirs, a combination of frac-pack in intelligent wells with multi-completions has been applied
[24].

14
Causes of Sand Production
2.1: The behaviors of unconsolidated formations
Both consolidated and unconsolidated formations produce sand [25], but most of the sand
production issues relate to unconsolidated formations. For the rest of the dissertation, focus will
be on unconsolidated formations where little to no inter-granular cement exists. Note that in the
dissertation, terminologies of fine migration, particle migration, and grain movement, etc., are
inter-changeable, as the sizes of moving particles are determined by apparent fluid velocity. As
long as the local apparent fluid velocity is high enough, not only the formation fines, but also the
uncemented, loose grain particles of bigger sizes can be mobilized.
Recovery of hydrocarbon from unconsolidated sands poses operational challenges such as
sand erosion, water channeling, development of thief zones, and potential formation subsidence
as a result of its inherent stress-sensitivity. Some of the observations are discussed below.

2.1.1: Dynamic reservoir properties as a function of pore pressure, permeability, skin,
and compressibility
Adams (1983) [26] analyzed the data from the Kuparuk field and conducted pressure buildup
and drawdown analyses. The study also documented core analyses with flow experiments, and
history match to the field data. The author demonstrated that, in high permeability sandstone
reservoirs, permeability behaves as a function of pore pressure. He indicated that permeability
estimations from buildup and drawdown tests could be 2 to 7 times different from each other.
Once correlated, the permeability - pore pressure curves matched the simulation results, along
with the pressure drawdown and buildup data. Moreover, he indicated that the skin estimation
assuming a constant permeability could be very misleading, which often results in unnecessary
acidizing jobs at high costs.
15
It should be noted that in soft and unconsolidated sands, the nature of dynamic permeability
may also invalidate many calculations which assume a constant permeability. We will discuss the
problem in additional publications, besides the dissertation which focus on sand migration induced
permeability changes, which affect history matching of waterflooding performance in
unconsolidated formations.
Yang et al. (1998) [27] pointed out the difficulty in history matching the data from
unconsolidated sands. The dynamics of the compaction-rebound effect in stress-sensitive
formations must be incorporated into simulation so that satisfactory matching can be achieved.
Here, rebound refers to the elastic behavior of high-compressibility formations. Porosity,
temperature (especially for heavy oil recovery), and principal stresses all play their roles in
affecting the compressibility. Reservoir pressure affects formation compressibility as well,
especially when conditions of reservoir depletion and water supporting change with time.

2.1.2: Subsidence
In unconsolidated formations, pore pressure balances overburden stress and confining
stresses from other directions. If pore pressure decreases, a high-permeability, high
compressibility formation may subside [28]. In application of waterflooding or other pressure
maintenance methods, one needs to constantly manage the voidage replacement ratio (VRR). It
helps avoid production of excessive solution gas and prevent reservoir subsidence in poorly
consolidated systems. Therefore, constant monitoring of injection wells is a key to keep the VRR
as needed. It is important to note that injectivity loss may happen due to reservoir pressure buildup,
or skin development caused by issues such as water incompatibility, clay swelling, high total
dissolved solids (TDS) in injection fluid, and fine migration.

16
2.1.3: Grain movement
Grain movement has been verified based on lab measurements [29] and observations of
sand productions in the field [12]. An operational practice has been to periodically backflush or
backflow the producing wells where one can temporarily remove the near-wellbore damage and
enhance the fluid flow. Through backflush, pore plugging in near-wellbore region gets reduced as
damaging particle transports with fluid flow and travels farther into formation. Grain movement
can also be reflected in Drill Stem Test, as when undergoing a period of fluid flow, the formation
gets cleaner.
As water injectivity decline is analogous to filtration process, laboratory studies on filtration
are of great value to modeling of grain movement. More specifically, experimental studies on
gravel packs and frac packs, provide opportunities to understand the development of formation
damage under waterflooding. Several approaches and data sets have been published on the
subject.
Gravel packs and frac packs act as a filter, helping hydrocarbon recovery from an
unconsolidated formation from two perspectives: preventing sand production, fast abrasion of
tubular and equipment, as well as avoiding subsidence. The principle is that the gravels selected
of certain sizes do not form a pathway for majority of moving sand particles. If no sand control
method is taken, formation particles in unconsolidated layers will transport with fluid flow when
flowing rate is high enough.
Therefore, clear evidence of grain movement can be obtained from the studies on sand
control in unconsolidated formations using gravel-packs and frac-packs subject to hydraulic forces.
Experiences have been documented in regard to applying frac-packs.
The study by Monus et al. [30] described the use of tip screen out technique in overcoming
the difficulty of fracturing soft rock, which also suggested on the selection of proppant in fracturing
the unconsolidated formation of the Gulf of Mexico. 40/60 gravels were selected for the gravel
pack placement. Shumbera et al. [21] documented the utilization of frac-pack method in the
17
injection well in the Offshore Gulf of Mexico to improve its performance. Moreover, the study
mentioned self-cleaning capability of screen/filter systems by backflushing, which, besides frack-
packing, stands as another indication of grain movement. Norman (2004) [31] highlighted the
longer lifespan of frac-pack over other sand control methods.

2.1.4: Grain breakage
Unconsolidated formations are naturally stress-sensitive as reflected on their mechanical and
fluid properties. When pore pressure changes in such soft rocks, the net stress change is drastic,
which can lead to several physical processes. In addition to grain movement, grain rotations and
grain breakage may also happen. Davies and Davies (2001) [32] discussed the fracturing of
individual grain on several mineral samples such as feldspar and multi-mineral grains. Note that
for sand grains, breakage of a single particle happens only at high-stress levels.

2.1.5: Thief zone
The thief zone issue has been one of the major concerns in thermal stimulation of heavy oil
sands and waterflooding operation. In the dissertation, by thief zone, we mean the injector-
producer communication as a result of gradual sand erosion, which develops into washed-out
zones in unconsolidated sandstones subject to water injection. Also, the concept should not be
mixed with the modeling of a penny-shaped artificial fracture in low-permeability rock, which has
been extensively presented in literature, such as Sneddon, I. (1946) [33]; Biot, M. A., Masse, L.,
& Medlin, W. L. (1986) [34]; Abousleiman et al. (1994) [35]; Zhang et al. (2002) [36], etc.
Unconsolidated formations, due to the natural lack of cementation, may develop into thief
zones when undergoing to continuous waterflooding. Thief zones could form in either single or
multi-layered unconsolidated sands. Especially in multi-layered reservoirs, development of
washed-out zone(s) can lead to uneven injectivities among all layers. Issues associated include
uneven injectivity profiles, poor vertical sweep efficiency, and also sand production-related issues.
18
One remediation measure is selective plugging using gels or polymer-based materials [37] [38].
Other measures such as mechanical shutoff, are discussed in Section 1.1.
Ideally, thief zones or the layer with gradual increase in permeability among stacked sands
should be recognized early. The Hall Plot method [39], based on the Darcy's equation, can be
used to monitor injectivity changes during waterflooding. However, the use of the Hall plot solely
loses its validity for multi-layer sands if the layers are undergoing commingled injection. One has
to consult profile surveys periodically to monitor the dynamics in injectivity of each layer, as
proposed by Qi et al. (2017) [40]. The cases from California and Kuwait demonstrated the
necessity of incorporating profile surveys into the interpretation of injectivity variations for multi-
layered systems.
Figure 2.1 and Figure 2.3 [40] show two Hall plots from a Kuwait field, with one case for a
vertical injector and the other for a horizontal injector, respectively. Although the slope changes
are associated with changes in injectivity, the plot itself has a limited resolution to identify the
troubled zones among all the layers under flow. Fortunately, profile surveys were available for
both cases, as shown in Figure 2.2 and Figure 2.4.
For the case of the vertical injector, Table 2.1 and Figure 2.2 show that, the two sections
received more than 99 percent of the fluid, , although they only represent about 50% thickness of
the total interval. The injectivity losses were much more affected by those high input layers. It
became clear that, in order to take remedial measures properly, one must have access to profile
survey data. Figure 2.1 shows the overall changes of injectivity but fails to deliver more detailed
information on the injectivity of each layer.
Similar observations were made for the case of the horizontal injector, as shown in Table 2.2
and Figure 2.3. Again, it is the profile survey in Figure 2.4, indicated that less than 23% of the
entire interval is taking more than 86% of the injected water, which could not be identified just
through the Hall Plot shown in Figure 2.3.
19
In absence of a reservoir pressure buildup, accumulation of fluids around an injector requires
an increase in injection pressure. Derivation of a new Hall plot equation for multilayer formations
and the simulation work highlighted limitations of the diagnostic functionality the traditional Hall
plot formulation used for single-layered reservoirs [40]. Profile survey, as a key technique to
disclose information of injectivity variations among all layers, can be run periodically to
demonstrate the injectivity changes of each layer over time.
Although the life of a reservoir can last for decades, obvious thief zone(s) in unconsolidated
formations can be developed in a short period of time, such as a few months. Figure 2.5 and
Table 2.3 [41] demonstrate how rapidly the injectivity profiles changed in unconsolidated sands
in the Huntington Beach oilfield in California. Figure 2.6 [42] is another example of the changes
noticed by comparisons among profile surveys.
As illustrated by the figures, certain layers in unconsolidated formations may receive a higher
fraction of water injection, which lead to insufficient pressure maintenance in other layers, and
water circulation between surface and subsurface through the thief zones. The challenge is to
detect eroded intervals before one of them becomes a thief zone. Early detection of thief zones
requires checking profile surveys periodically in multi-layered formations.
Related diagnostic approaches have been constructed for various conditions. Calhoun (1953)
[43] and Watkins (1973) [44] proposed mathematical methods to diagnose a pancake-shaped
thief zone and its areal extension. However, the techniques assumed steady-state conditions in
the calculations, after water breakthrough into producing wells.
Pressure transient analyses done by Gringarten and Ramey (1974) [45], Valko and
Economides (1997) [46], etc., were proposed in modeling the system where a high conductivity
layer of a certain radius and thickness is embedded inside the reservoir.
See examples shown in Figure 2.7, Figure 2.8, and Figure 2.9, profile surveys of the three
wells demonstrated gradual development of a washed-out zone. More informatively, pressure fall-
off data is available, as shown on the log-log derivative plots. Due to the nature of unconsolidated
20
formations, pressure fall-off tests for the situation, where a thief zone has been developed in the
high-permeability sand matrix, indicate radial flow, illustrated by the zero-slope on pressure
derivative plots [47].

Figure 2.1: Hall Plot For The Vertical Injector From Kuwait. ([40] Courtesy SPE.)

Table 2.1: Water Profile Survey Results For Vertical Injector From Kuwait. ([40] Courtesy SPE.)
Section
Perforation Zone,
ft.
Thickness,
ft.
Water Injected,
bbl./day
Water Injected,
%
1 7930 - 7956 26 13,676 38.2
2
7973 - 7983 10 204 0.6
7983 - 8010 27 21,954 61.3
8010 - 8019 9 0 0.0
3 8025 - 8040 15 0 0.0
4 8059 - 8066 7 0 0.0
Cumulative 94 35,834 100.0
0
500,000
1,000,000
1,500,000
2,000,000
0 2,000,0004,000,0006,000,0008,000,00010,000,000 12,000,000
Cumulative Pressure, psi*days
Cumulative Injected Water Volume, bbl
Hall Plot of the Vertical Injector
21

Figure 2.2: Water Profile Survey For The Vertical Injector From Kuwait. ([40] Courtesy SPE.)
Water Profile Survey for the Verical Sample well showing the majority of the water injected is being sent to
the first 2 sections while the remaining two sections did not receive water.
22

Figure 2.3: Hall Plot Of The Horizontal Injector From Kuwait. ([40] Courtesy SPE.)

Table 2.2: Water Profile Survey Results For The Horizontal Injector From Kuwait. ([40] Courtesy
SPE.)
Section
Perforation Zone,
ft.
Thickness,
ft.
Water Injected,
bbl./d
Water Injected,
%
1 8858 - 9100 242 3,498 27.2
2 9100 - 11300 2200 140 1.1
3 11300 - 11450 150 1,610 12.5
4 11450 - 11914 464 7,618 59.2
Cumulative 3056 12,866 100.0
0
200,000
400,000
600,000
800,000
1,000,000
1,200,000
1,400,000
1,600,000
0 500,000 1,000,000 1,500,000 2,000,000 2,500,000
Cumulative Pressure, psi*days
Cumulative Injected Water Volume, bbl
Hall Plot of the Horizontal Injector
23

Figure 2.4: Water Profile Survey For The Horizontal Injector From Kuwait. ([40] Courtesy SPE.)

Water Profile Survey for the Horizontal Injection well showing the majority of the water injected is being sent
to the 1st, 3rd, and 4th sections while the second section (being the longest interval) received only 1.1% of the
total water injected.
24

Figure 2.5: Effect Of A Shock Treatment On Profile And Injectivity Of A Well. ([41], Courtesy
SPE.)

Table 2.3: Injectivity Increase Associated With Well Profile In Figure 2.1. ([41], Courtesy SPE.)
Run No. 1 2
Date, mm/dd/yy 09/30/66 04/03/67
Rate, B/D 2880 4800
Pressure, psi 1000 1000
Cum. Injection, Million Bbl. 0.13 0.7

Washed-out
zone
develops
0 50 100
Percentage of
Total Injection
5000
5200
5400
5600
Depth
25

Figure 2.6: Injection Profile Behavior At Various Rates And Pressures, Lower Jones Injector A,
Huntington Beach Offshore Field. ([42]. Courtesy SPE.)

Figure 2.7: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well A. ([47] Courtesy
SPE.)

26

Figure 2.8: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well B. ([47] Courtesy
SPE.)

Figure 2.9: Profile Survey And Derivative Of Pressure Fall-Off Data Of Well C. ([47] Courtesy
SPE.)

The derivative plots in Figure 2.7, Figure 2.8, and Figure 2.8 differ from the solutions
proposed by Gringarten and Ramey (1974) [45], Valko and Economides (1997) [46], etc. The
latter ones describe the system where a high conductivity layer, of a certain radius and thickness,
is embedded inside the reservoir. Although Gringarten and Ramey (1974) [45] assumed uniform

27
flux distribution in the fracture, Valko and Economides (1997) [46] further considered finite
conductivity and infinite conductivity of the fracture, both studies assumed that the host rock
contributes little to fluid flow. This is unlikely to be the case for unconsolidated sand, as both the
thief zone and the sand matrix are of high permeability which contribute to fluid flow. The
difference between two systems is illustrated in Figure 2.10.
Satman (1981) [48] proposed a transient pressure solution in modeling uneven water fronts
in multi-layered formations, with the assumptions not limited to that host rock does not contribute
much to fluid flow. While in detecting gradual development of a thief zone, periodic use of pressure
fall-off tests fits better in modeling the dynamics, as proposed by Qi and Ershaghi (2019) [47].

Figure 2.10: Schematic Of A Reservoir With Horizontal Fracture (A) Vs. With A Horizontal
Washed-Out Zone In Unconsolidated Formation (B). ([47] Courtesy SPE.)

In brief, unconsolidated formations under waterflooding for pressure maintenance purposes
may go through sand erosion and development to thief zones, where high permeability intervals
with little to no cementation can be washed out gradually. Before immediate remediation
k
eroded
= k
2
h
eroded
= h
2
r
e
r
eroded
a. No flow entry
into host rock
b. Flow entries both host
rock and eroded zone
k
hp
= k
2
h
hp
= h
2
r
e
r
hp
k
r
k
r
28
measures such as sub-zoning or profile correction using chemicals, one of the ways to identify
thief zone(s) is to consult periodically measured profile surveys. Also, repeated fall-off tests can
help in early detection of eroded zones at a lower cost [47]. It is noted that log-log derivative plots
of the fall-off test data show a radial flow, which is different from the several pressure transient
solutions presented in the literature, addressing a similar situation except for the hard rock
considered.

2.2: Rock characteristics of unconsolidated formations
2.2.1: Formation strength, mineral composition of cementation and formation fines
Several factors that control formation strength include inter-granular cementation, inter-
locking of sand grains due to vertical and horizontal confining stresses, as well as interfacial
tensions among formation fluids, and wetting-phase saturations. What differentiates
unconsolidated formations from consolidated formations is the significantly different levels of
cementation precipitated around sand grains. The cementation materials involve silt and clay, and
sometimes calcite, that have precipitated over a geologic time under formation temperature, which
bond sand grains together. Mineral analysis of sidewall cores and well cleanouts out of
unconsolidated formations show poor cementation among sand particles [15].
Formation strength can be estimated quantitatively from density log and acoustic velocity log
[22]. More precisely, the estimation done for the formation sand near wellbore matters the most
for fields with sand production issues and associated sand control measures. Shear modulus,
(𝐸
/
)
0
as a reflection of relative strength of a formation at certain depth 𝑑, can be obtained by,
*
<
=
𝐸
9
,
>
= *𝐸
3
+
<
=
𝐸
9
,
>
− *𝐸
3
+
<
=
𝐸
9
,
?
(Eq. 2.1, [22])
where,
*𝐸
3
+
<
=
𝐸
9
,
?
≈ (𝐸
3
)
>
(Eq. 2.2, [22])
29
where, subscript A refers to asymptotic conditions at surface or a shallow depth where the
shear modulus of the formation is negligibly small. Illustration of combined modulus 3𝐸
%
+
1
2
𝐸
/
5
as a function of depth, 𝑑 is shown in Figure 2.11.
Estimated value of (𝐸
/
)
3
is also useful in predicting sand influx or in obtaining a critical flow
rate for the passive sand control method, as described in Section 1.3.7.

Figure 2.11: Dynamic Combination Modulus Curve As A Function Of Depth. ([22] Courtesy
SPE.)

In naturally occurring formations, fines cover a broad range of clay particles composed of
montmorillonite, illite, kaolinite, and chlorite, but also a wide variety of particles including, but not
limited to quartz, feldspar, muscovite, sodium chloride, calcite, dolomite, barite, and amorphous
materials. Clays of different types may not be the dominating type of fine particles, see Table 2.4.
Samples analyzed are particles that can pass through a 400-mesh sieve, with a pore size of 37
microns.
Fines are located in interior space of porous matrix. Unlike larger grains, which are bonded
by cementation of different levels, fines are not connected to bigger grains. Fine particles can
move subject to fluid flow through porous media. Figure 2.12 shows the relative sizes of the fines
and typical grain particles observed under the Scanning Electron Microscope (SEM) [49].
Test Zone
Depth
Zone Of Interest
Depth
Relative Strength Of
Zone Of Interest
Relative Strength Of
Test Zone
Asymptotic
Modulus
INCREASING DEPTH
Increasing Modulus
E
B
+ 4/3 E
S
30
Table 2.4: X-Ray Analysis Of GOM Unconsolidated Formation Fines. ([49]. Courtesy SPE.)
Well A Well B Well C Well D Well E
Clays
Montmorillonite 5.5 13.4 2.2 1.4 -
Illite 6.2 9.1 3.3 1.7 -
Kaolinite 0.8 4.2 1.3 0.7 -
Chlorite 3.9 - - - -
Quartz 36.7 24.0 47.3 17.0 68.3
Other Minerals
Feldspar 8.6 5.7 9.1 5.4 11.4
Muscovite 1.6 - 1.6 1.0 -
Sodium chloride 1.1 1.3 7.8 5.0 1.5
Calcite - 1.6 - - 1.5
Dolomite - - 1.8 2.8 -
Barite - - - 22.1 -
Amorphous Materials 35.6 40.7 25.9 42.9 17.3
Total 100.0 100.0 100.0 100.0 100.0

Figure 2.12: Examples Of Fine Particle Located On The Surface Of Larger Formation Sand
Grains. ([49] Courtesy SPE.)

31
2.2.2: Examples of unconsolidated sands with grain size distributions
In loosely cemented sands, unbonded fines, when carried by flowing fluid, move away from
the near-wellbore regions of injectors into a deeper reservoir toward producing wells, as described
in Civan (2017) [50]. Formation fines come from several sources, including deposition during
geologic time, drilling and completion fluids, or parts of total dissolved solids in injected water. For
describing grain size distributions of formations, coefficient C is defined as:
𝐶 =
"
45
"
65
(Eq. 2.3, [16])
where 𝐷
17
is grain diameter at 40 percentiles on cumulative percentage plot, and 𝐷
87
is grain
diameter correspond to 90 percentiles. For example, the grain size distributions of California
unconsolidated sand has a C value ranging between 2 and 8, as shown in Figure 2.13 and Figure
2.14, which represents uniform sands (with a C value less than 3) and nonuniform sands (C
greater than 5).

Figure 2.13: Uniform Sized, Unconsolidated California Sand. ([16] Courtesy SPE.)

0
20
40
60
80
100%
1.0
0.5
0.1
0.05
0.01
0.005
Cumulative Percentage By Weight
Grain Diameter, inches
Gravel Critical
Grain Size
0.102"
Formation Critical
Grain Size
0.017"
Formation Sand
Profile
Gravel
Profile
Slot Size
0.050”
0.001
C = 1.42
32

Figure 2.14: Nonuniform Sized, Unconsolidated California Sand. ([16] Courtesy SPE.)

Figure 2.15: Typical Grain Size Distribution Of GOM Unconsolidated Sand ([51]).

As another example, Figure 2.15 shows a grain size distribution for the U.S. Gulf of Mexico
offshore unconsolidated sand.
It should be noted that the above grain size distributions for unconsolidated sands may not
fully represent the fine particles and their amounts. Although the fine particles seem small and
sparse from Figure 2.12, total fines can actually weigh significantly. For instance, lab
0
20
40
60
80
100%
1.0
0.5
0.1
0.05
0.01
0.005
Cumulative Percentage By Weight
Grain Diameter, inches
Gravel Critical
Grain Size
0.096"
Formation Critical
Grain Size
0.016"
Formation Sand
Profile
Gravel
Profile
Slot Size
0.060”
0.001
C = 1.43
Cumulative weight, %
Grain Diameter, inch
0.1 0.01 0.001 1.0
50
100
0
Typical U.S.
Gulf Coast
Sand
33
measurements on the GOM unconsolidated sand samples show that the weight percentage for
grains smaller than 37 microns can be 2 – 15% through sieve analysis [49]. And the grains smaller
than such value is not clearly reflected in Figure 2.15, possibly due to a lack of data density
achieved in sieve analysis. Another possibility would be unrepresentative sampling done prior to
measurements [25].

2.3: Current understandings on sand particle movements and
injectivity issues for injection wells
Milton et al. (1961) [11] environmentally investigate the effect of the sizes and shapes of the
grains composing an artificial 'formation' sample, as well as its (swelling) clay content on
waterflooding improved oil recovery. All three variables were quantified using the synthetic
materials. It was observed that the optimized grain size leads to the greatest recovery
improvement result. The selection was tested among fine, intermediate, coarse glass beads, and
sand grains. The optimum values differ among the clean glass bead system, clean angular sand
system, and angular sand system with clay content.
Barkman and Davidson (1972) [12] mathematically showed that the bigger the contrast
between formation permeability and cake (damaged zone) permeability, the worse the injectivity
decline will be. It indicates that damage brought by water injection in high permeability,
unconsolidated sand formations is more severe than in formations with lower permeability.
More importantly, the study systematically proposed four conceptual mechanisms of good
impairment subject to waterflooding. Those include wellbore narrowing as a result of cake
developed on sand face; invasion in near-wellbore region of formation, as a result of fine-grain
movements carried by injected water; wellbore fill-up when solids deposit on bottom of a well; and
perforation plugging in a cased and perforated well. For each one of the mechanisms, a
34
mathematical model was proposed. As of today, Barkman and Davidson (1972) ([12]) remains to
be one of the most fundamental frameworks in modeling the process.
In fields, fines migration is noticed when the injection of seemingly non-reactive fluids (water,
oil, gas, solvent/surfactant solutions) can remove formation damage around the wellbore. Another
evidence is the operational practice of periodical backflows in both injectors and producers. The
process is interpreted as that the local pressure changes are high enough to mobilize formation
fine particles and break pore plugging.
It should be noted that not all formation fines migrate with fluid. Reasons for immobilization
of formation fines include:
1. plugging of pore throats by fines or accumulation of fines;
2. local low flow rate; and,
3. multi-phase presence, where the flowing phase has not touched the relatively steady
phase, where the latter wets and coats fine particles and matrix grains [49].
As such, the dynamic balance between pore throat plugging (bridging) and fine migration are
of equal importance in modeling thief zone development. Work on fine-migration related
processes has been addressed in several types of studies. This includes field measurements of
the grain size distribution of the produced sand from Huntington Beach [16], lab measurements
of migrating fines during continued injection of single and multi-phase flow on samples cored from
the U.S. Gulf Coast [49], and explainations of fine migration mechanisms under various wettability
and miscibility conditions.
Davidson (1979) [13] also mentioned that the tests with a condition of constant pressure drop
might not be able to model particle sizes and their distribution, which implies the need for
additional experimental work and appropriate design of a representative sample being flooded, in
addition to the SEM examination of particle composition suspended in water.
The rate-controlled flooding done on the field, which is at neither the condition of constant
pressure drop, nor the condition of constant rate. Ershaghi et al. (1986) ([53]) stated that the
35
action of increasing injection pressure, which happens in field operations when injectivity
decreases, might cause further damage to formation if invasion is the dominating mechanism.
The study also noted the extra core damaging effect if adding another chemical component in
injected water such as corrosion inhibitor.
Hofsaess and Kleinitz (2003) [14] summarized the field data for over 30 years. Several key
points were emphasized, such as that the completion type of a wellbore significantly affects the
extent of injectivity decline, especially at the beginning of a well’s injection life. However, the
calculated injector's half-life can be different from the models presented by Barkman and
Davidson (1972) [12] for orders of magnitude. The errors emanate for reasons include the usage
of membrane filtration in the lab to justify the water quality and the simplification that Barkman
and Davidson (1972) ([12]) made using a less commonly applied condition of constant pressure
applied to the four mechanisms.
Besides the theory that the filter-cake developed is the main cause of injectivity loss, analogy
to deep filtration process has been mentioned in many studies. The phenomena that several
mechanisms of particle capturing happen simultaneously was discussed [54] [55], while the
mechanisms were often correlated separately. Ali et al. (2009) [56] experimentally studied the
residual oil affecting deep filtration process. They noted that the invasion into the residual oil-
saturated formation is deeper compared to the system saturated with pure brine. The study also
acknowledged the possibility for particles attached to oil droplets being left in pore throats.
Guedes et al. (2009) ([57]) proposed an analytical solution based upon the classical
governing equation in mass transport to model deep filtration. Instead of quantifying the several
mechanisms of particle movement with individual coefficients, the solution determined one
combined coefficient to fit the field data into a correlation.
Vaz et al. (2010) [58] upscaled the analytical solution for deep filtration problems from the
core level to the level covering near-wellbore region. Moreover, the study considered
36
heterogeneity, which has been recognized to affect well impairment significantly ([12], [14]). They
presented the solutions for in-parallel, in-series, and combined types of heterogeneity models.
Numerical models have been presented more recently and progressed rapidly. Bruno et al.
(2001) [19] simulated fracturing with coupled particle and fluid flow approach in less-consolidated
formations. The study provided useful information on the mechanical properties of unconsolidated
sands. The conclusions mentioned a transition zone that exists from brittle fracturing of weakly
cemented sands to dilation and inelastic parting. Shear bond in such sands, while being very
weak, controls the transition zone.
In more recent work, dynamics in rock properties have been experimentally tested and
correlated. Davies and Davies (2001) ([32]) tested the change in stress-permeability curve for
three rock types, and each has its own cement level, mineralogy, grain packing type, and grain
angularity. Moreover, the less discussed physics of intra-grain cleavage has been mentioned. In
comparison to the terminology of ‘failure’, which refers to breakage of inter-grain cement,
‘cleavage’ refers to the fracturing of a single grain. Although the high stress level is unlikely to
happen in typical waterflooding operations, the dynamics of properties needed to be included in
numerical modeling.
The grain movement process has been numerically modeled by Barake (2015) [59], where
the grain depositing, grain being eroded, and grain staying with matrix are considered as different
phases in addition to the water phase. The work is applied to the channeling in unconsolidated
formations where grains are moving when adjacent to the deposition/erosion surface of the sand;
thus, formation emptiness gradually forms until development of channels. However, the work did
not honor the variations in grain size distributions and pore size distributions, as well as their
relationship that may further affect particle migration process.
Sharma et al. (1985) [60] presented an analytical solution based upon the Effective Medium
Theory to predict permeability changes during formation damaging and design the stimulation
jobs. Ohen and Civan (1991) [61] numerically presented a model the near-wellbore behavior with
37
proven damage, in order to conduct a good plan for reservoir development. Bedrikovestsky et al.
(2002, 2003, 2004) [62][63][64] proposed solutions to particle detachment based on mechanical
equilibrium, which included torque balance of drag, electrostatic, lifting, and gravity forces.
Zeinijanhromi and Bedrikovetski (2011) [3] proposed analytical solutions for the cases under
single- and two-phase flow. On a large scale, fine migration is equivalent to polymer flooding for
mobility control.
Numerical modeling of sand erosion processes, including channelization in an
unconsolidated reservoir subject to water injection, has also been presented Rostami and
Taleghani (2014) [52]. In their work, they considered both erosion and deposition processes
during fluid flow, primarily focusing on sand erosion as a result of high injection rate induced sand
failure or fracturing, which eventually leads to channelization, rather than the pressure conditions
in typical water injection. Afrough et al. (2017) [65] used the Magnetic-resonance–imaging (MRI)
method to present an analytical solution to screen reservoirs with potentials for fine migration.
Karazincir et al. (2017) [25] conducted lab experiments under various test conditions to fill out the
gap between lab and field observations. Borazjani et al. (2018) [66] proposed a semi-analytical
solution to model injected water salinity, fine migration, and permeability losses under two-phase
flow. In another study by Wennberg et al. (1995) [67], mass transfer governing equations are used
for convection and diffusion processes, to improve the calculation. Similar to Barake (2015) [59],
mobile grains (fines), clogged drains, and steady grains are considered separately, but only fines
can move through pore throats.
These solutions were proposed for problems associated with fine migration and permeability
decreases, while specifically none of them quantitatively modeled the phenomena mentioned in
the dissertation, where formation sand gets eroded due to the movement of dislodged particles
during waterflooding, and both pore size distribution and grain size distribution change. The
dynamics in the distributions are the key to changes in reservoir properties, which matter the most
38
for history matching the performance of waterflooding in unconsolidated formations. Chapter 4 to
7 will present the algorithm.
Also, the situation of more layers participating in taking injected water and adds into increased
injectivity were documented. From the profile surveys shown by Yoelin et al. (1970) ([42]), when
increasing the injection pressure from 340 psi to 570 psi, more layers toward the bottom of the
multi-layer system participated in taking the injection fluid. In between the pressure range, a
formation parting pressure was estimated, while it may not be the case. The study in regard to
the parting pressure of the unconsolidated sands will be systematically reviewed in Section 3.2.

39
Injectivity Maintenance and
Review of Step Rate Test
3.1: Injectivity monitoring and permeability variations
In waterflooding operations, close monitoring of well injectivity is required to maintain
reservoir pressure, avoid solution gas coming out of oil, keep Voidage Replacement Ratio (VRR)
as needed. Reservoir monitoring are of key importance, especially for unconsolidated reservoirs
where failure in injectivity and pressure maintenance may result in severe safety issues such as
subsidence [23] [28].
Injectivity monitoring essentially means keeping track of changes in permeability and skin.
As a considerable number of wells are under commingled flow in multi-layered sands, good
surveillance is required in several aspects. Monitoring injectivity with respect to time and individual
layer for timely control of a watered-out layer, channeling, or a thief zone being developed. Also,
it could detect injectivity loss due to formation damage created by clay swelling, accumulation of
introduced Total Dissolved Solids (TDS), water compatibility issues, or fine migration. These
aspects ensure that each layer receives sufficient water injection to maintain the pressure.
Injectivity index, defined as the ratio of injection rate to pressure change, in BPD/psi, is used
to guide injectivity of a well. Pressure change refers to the difference between bottom hole
pressure at an injection well and reservoir external pressure. In reality, one may notice injectivity
change from changes in surface injection pressure while same amount of fluid is injected. If
surface pressure drops, permeability increase may have happened, and vice versa.
In multi-layered sands, permeability variations can be mapped by running profile surveys [39]
[40]. Examples of permeability changes indicated by the profile surveys are illustrated in Section
2.1.5, from Figure 2.2 and Figure 2.4 to Figure 2.9. While profile surveys quantitatively measure
40
the fraction of total fluid accepted by each layer, temperature logging is another tool to qualitatively
display permeability variations [68] among layers, when a gap between the temperature of
injected fluid and formation temperature exists. Figure 3.1 is an example of temperature logging
as an indicator of permeability profile around an injection well. The 'Static' curve refers to original
profile of formation temperature, which increases with formation depth. When a channel or a thief
zone is being developed, permeability of the layer increases so that the layer accepts more
injected fluid, a differential in temperature profiles can reflect the situation. Usually, time-
dependence of temperature profiles is consulted to further validate the observation.

Figure 3.1: Channeling Shown From A Temperature Log In An Injection Well. ([68]. Courtesy
SPE.)

3.2: Fracturing signals from SRT vs. grain movement in multi-
layered sands
A purposeful fracturing of tight geologic formations is conducted by injecting fluid mixtures at
a pressure above the formation parting pressure. This process is usually applied on tight sands,
carbonates, shales, and diatomite, etc. The goal is to create high-conductivity pathways to recover
hydrocarbon economically from the low permeability prospects. Sometimes pumping of fluid
Flowing
Producing Zone
Channeled Zone
2-Hour Shut-in
Static
Temperature Increase
41
mixtures at or above fracturing pressure is conducted in high-permeability, unconsolidated
formation, such as in frac-pack placement [28], or in resin consolidation job to secure sand
particles in place, mainly for sand control purposes [68]. It may also happen in polymer injection
or squeezing shut-off fluids in order to control excess water production [1][10].
In waterflooding operations without the intention to create fractures in formation, finding
formation parting pressure is still important. Operators tend to maximize the injection rate for each
injection well, but the rate cannot be too high so that formation gets fractured.
Determination of formation parting pressure can be done by running the Step Rate Test (SRT)
[69] [70]. By injecting fluid at various rates and pressures, a plot of stabilized injection pressure
vs. injection rate can be generated. Figure 3.2 shows a typical example of Step Rate Test. A
sharp decrease in slope on the pressure-rate plot is interpreted as the signal of formation
breakdown. The technique is widely applied, which provides a guidance, so that during
waterflooding the injection is operated at a maximum allowable injection rate or injection pressure
gradient without fracturing formation. The breakdown pressure gradient varies with the
geomechanical properties of formations [70][71].

Figure 3.2: Illustration Of Stabilized Pressure –Rate Plot From SRT. ([68]. Courtesy SPE.)

2000
Pressure, psi
Injection Rate, B/D
Fracture
Pressure
2400
2800
3200
3600
0 500 1000 1500 2000 2500 3000
42
3.2.1: A critical review of finding formation breakdown pressure in multi-layered
reservoirs
A case using pressure – rate data is determining formation fracturing pressure in multilayered
sands is critically reviewed in this Section. Referring to the case shown in Figure 2.6, details on
surface injection pressures and rates in the polymer injection are shown in Table 3.1 [42]. The
study was conducted to provide instructions on the polymer treatment for an injection well, around
which a thief zone had been identified. However, when a multi-layered system is tested to find
its formation breakdown pressure, several facts should be realized:
1. Each layer among the stacked sands has its distinct geomechanical characteristics. As
such, running a step-rate test non-selectively to all layers may lose its validity, as one can
anticipate that not all layers break simultaneously while they carry different geomechanical
properties. Therefore, in multi-layered sands, each layer has its own parting pressure, and
conducting tests on each of the members separately is more accurate.
For the case shown in Table 3.1 and Figure 2.6, when surface pressure increased from 340
psi (Run 3) to 570 psi (Run 5), the formation parting pressure was estimated to be 470 psi. If the
interpretation was correct, the first question would be, which one(s) among the entire stacked
sands was fractured.

Table 3.1: Data From A Multi-Profile Test Before Polymer Treatment Associated With Figure
2.6. ([42]. Courtesy SPE.)
Run
No.
Wellhead
Pressure, PSIG
Total Injection
Rate, B/D
Injection Rate into
Thief zone, B/D
Injection Rate into
Other Sands, B/D
1 220 400 360 40
2 260 700 320 380
3 340 1000 380 620
4 440 1500 380 1120
470*
5 570 3300 1300 2000
Note:

*refers to formation parting pressure [42].

43
2. A slope decrease on pressure – rate plot means increase in the cumulative product of
permeability and thickness of stacked sands, ∑ 𝑘
9
ℎ
9 9
, which does not necessarily indicate a
fracture created. Here, 𝑘 and ℎ are permeability and thickness, respectively, and the subscript 𝑖
refer to the index for each among all sand layers. From Figure 2.6, it is noted that when the layers
located at 3280 – 3300 ft. and 3500 - 3620 ft. participated in accepting more injected fluid from
Run 3 to Run 5, the cumulative product ∑ 𝑘
9
ℎ
9 9
increased, which has the similar slope-decrease
effect as when a mini-fracture is created. Having more layers participating adds additional 𝑘
:
ℎ
:

into cumulative product ∑ 𝑘
9
ℎ
9 9
, while the creation of fracture(s) means an increase in effective
permeability 𝑘 of one or more layers, which also makes cumulative product ∑ 𝑘
9
ℎ
9 9
higher. Two
situations lead to the same slope decrease effect on pressure-rate plot from Step Rate Test, as
shown in Figure 3.2.
3. The study [42] also mentioned that partial collapse of the stacked unconsolidated sands
had happened, so that polymer fluid was injected to modify the profile. Therefore, the cause for
the other sand members that lost their fluid intake capability may be one or both of the following
situations:
(a) A thief zone has developed at depth of 3320 ft., due to the accumulation of formation
grain migration and sand erosion, so that the permeability of the layer becomes apparently
higher than others, and more fluid flows into the thief zone compared to other layers.
(b) Other layer members may have experienced formation damage or partial collapse as
mentioned, and the loss of permeability would result in loss of injectivity.
In a stacked system, appearance of each layer on injectivity profile depends on the fraction
of individual permeability-thickness product over cumulative permeability-thickness product of all
layers,
;
!
<
!
∑ ;
!
<
! !
. As mentioned above, either case (a) or (b) occurs, or when both of the cases hold
true, profile surveys will behave in similar way, as the injected fluid is prone to get into one among
all layers.
44
3.2.2: Particle migration
Sand particle migration may not only result in damage of formation [48] [59]; may also lead
to the development of a washed-out zone, a thief zone, or fluid channeling [42] [57]. Same
behavior of sand migration leads to seemingly opposite results for two different mechanisms. In
the case as formation damage happens, migrated fine particle partially or fully bridges pore
throats along their pathway of movement [47]. While in the case of thief-zone development,
dislodged sand particles can pass through pore throats and leave more voids in local area. It
widens the pathways that ease both fluid flow and migration of following particles. Accumulation
of continuous sand loss gradually leads to development of a thief zone or fluid channel.
Occurrence of particle migration for distance requires not only big-enough pore throats that
do not block moving particles, but also in-situ pressure surges that are high enough to dislodge
formation particles. As mentioned in Section 2.3, it has been a common operational practice that,
periodically backflow is conducted in injection wells and producing wells by injecting non-acid
fluids for the purpose of temporarily removing formation damage in the near-wellbore region to
enhance fluid flow. The non-acid fluids injected are not reactive to dissolve formation particulates.
Instead, the fluid was injected to mobilize the damaging particles in the near-wellbore area. A
necessary condition for particle mobilization is that the local pressure perturbation or apparent
fluid velocity must exceed the critical threshold of the particles.
The routine practice of backflows stands as an evidence for particle migrations, in addition to
lab observations [47]. In the case mentioned in Section 3.2.1, by increasing injection pressure at
surface, more layers participated in taking the fluid [42], which occurrence may follow the same
mechanism. In summary, in multi-layered systems, for either creations of a fracture or having
more layer's participation in taking injected fluid, a pressure-rate plot from Step Rate Test would
have same behavior as slope decrease. Without additional information, the pressure at the
'bending' point on the curve may not indicate a true fracturing signal of stacked sands. It could
mean additional layers get activated in taking fluid injection.
45
Modeling of Sand Grain
Movement

From the discussions in Section 1.3.7 and Section 2.1.3, it is anticipated that factors including
but not limited to pore size exclusion and hydraulic forces work together to dislodge formation
particles and cause their movement. This Chapter describes a new numerical model as a
representation of sand erosion process induced by particle movement.
The modeling includes four parts:
1. Modeling of fluid flow in a reservoir in terms of pressure distribution based on the
discretized diffusivity equation.
2. Modeling the movement of unbonded/suspended particles along with the fluid flow
between an injection well and a producing well, where the flowing velocities of particles are
calculated as apparent flow velocities in porous media with the injection well completed as open
hole. The assumption is that, once particle velocity exceeds the critical velocity calculated for a
particle of certain size, the particle will be mobilized and migrate along with the fluid flow as long
as the pore size allows. As particles of various sizes move, they may clog some pore throats when
particle sizes exceed the sizes of the pore throats they encounter. This way the local grain
accumulation and permeability reduction happen. On the other side, if the particles pass through
the relatively bigger pores, sand erosion and permeability increase take place. The particles keep
moving until blocked by smaller pores later or get produced. New grain size distribution, as weight
percentage of local particles of various sizes, is calculated as the net result of grain accumulation
and loss at the current location and time. The model allows particles to be produced out of the
producing well if they are able to travel from original place to the producing well.
46
3. Pore size distribution varies following the Markov Chain probability transition process.
Giving possibilities to scenarios of pore size changes is necessary. In reality, either pore
enlargement or shrinkage could happen. Suppose a small particle gets detached off the pore wall
and pass through, the pore throat may get enlarged, which further eases the migration of bigger
sized particles afterward. As a comparison, if the particle encounters a relatively smaller pore
throat, it may partially or fully clog the pore throat. The pore clogging then negatively affects the
movement of other particles. Therefore, pore size variation is modelled stochastically, considering
pore enlargement, pore size shrinkage, and full plugging.
The initial selection of transition matrices for the Markov Chain process could be done by
honoring the several cases of formation pore clogging, enlargement, both, or no change in the
pore sizes, as defined in Section 4.2. Also, multiple realizations would be generated by assigning
various probability transition matrices satisfying the requirement associated with the chosen
scenario(s). For the realizations generated, values of ultimate permeability increase can be
calculated to compare with profile surveys, if available. If they do not match, the probability
transition matrices will be re-selected until validated by profile surveys. The stochastic process
allows inverse modeling and validation of lab/field measurements.
4. Due to the net changes in grain population densities, porosity and permeability at each
location are updated for next timestep. When net grain population densities decrease at current
timestep, an increase in porosity and permeability make particle movement easier in the next step,
and vice versa.
Assumptions, referred and newly defined equations, and the Markov chain process with
assigned stochastic matrices for various scenarios, are discussed in detail, followed by the
demonstration of the entire algorithm.

47
4.1: Analytics
4.1.1: Assumptions
The set of equations is defined with the following global assumptions:
1. Laminar and linear flow in a one-dimensional system is simulated. Initial properties are
homogeneously distributed over the reservoir, including porosity, permeability, grain size
distribution, and pore size distribution. As the problem is to examine particle migrations under
typical waterflooding conditions without the intention to fracturing formation, non-Darcy flow is not
considered [25].
2. Gravity is not taken into the calculation of reservoir pressure distribution.
3. Freshwater is the only fluid injected. The temperature of injected water is the same
formation temperature, as continuous circulation of injected water makes water temperature
increase to the level of formation temperature. Depth-temperature function affects kinematic and
dynamic viscosities of water. However, at a certain depth, formation temperature holds constant.
4. Formation grains are water wet and spherical. They vary in size and follow a grain size
distribution. No inter-granular cement exists.
5. Several possible processes such as grain breakage or clay swelling, which may change
the size and geometry of individual particle, are not considered.
Other specific assumptions of equations are discussed in the following sections.

4.1.2: Pressure distribution and apparent fluid flow velocity
The diffusivity equation is honored in calculating pressure distribution 𝜕𝑃(𝑥,𝑡). Thus, the
assumptions of laminar flow and no gravity at isothermal condition are valid. As sand
accumulation and erosion happens, both local porosity and permeability evolve. Porosity 𝛷 and
permeability 𝑘 are functions of time 𝑡 and location 𝑥.
48
@1(A,C)
@C
=
D(A,C)4E
>
. (A,C)
@
?
1(A,C)
@A
?
(Eq. 4.1)
where, 𝜇 and 𝐶
C
refer to dynamic viscosity and total formation compressibility, respectively.
With a pressure distribution 𝜕𝑃(𝑥,𝑡) calculated, apparent fluid flow velocity 𝑉
@
, at each
location and time, is calculated as,
𝑉
F
(𝑥,𝑡) =
. (A,C)
4D(A,C)
∇𝑃(𝑥,𝑡) (Eq. 4.2)
It is extremely important to note that many studies have pointed out that the apparent fluid
flow velocity becomes much higher in a perforated hole than in an open hole with same wellbore
radius, due to limited flow area caused by perforations [15] [17]. A diagram in estimating apparent
flow velocity at various wellbore radii and completion conditions is shown in Figure 4.1 [15].
Chapter 5 and 6 computes a case for open-hole condition. Illustration of incorporating
perforation condition into simulation will be included in Chapter 7.

Figure 4.1: Apparent Flow Velocity At Various Completion Conditions.
([15] Courtesy SPE.)

4.1.3: Critical shear velocity for particle migration
Besides Stoke’s law, many correlations have been proposed in calculating critical shear
velocity as a threshold for grain mobilization. Paphitis’s correlation (2001) [72] is one of the most
10
-4
10
-3
10
-2
10
-1
10
0
10
1
10
2
BPD/FT
Velocity (ft/s) in Gravel Pack
100
1000
10000
1 2 3 4 5
LEGEND
Velocity Based on
Flow Rate For Cross
Sectional Area
Specified By:
1 – 18” Open Hole
2 – 12” Open Hole
3 – 6” Open Hole
4 – 4- .75” Perf
5 – 4- .5” Perf
49
plausible ones, where critical shear velocity 𝑢
∗
"
is calculated iteratively using Eq. 4.3. (Paphitis,
2001) [72].
Note that the correlated critical movability function
B
∗
.
$
, instead of acting as a single curve, is
actually a range for particle movement criteria. According to the definition, the lower limit of the
function refers to mobilization threshold of single particle, and the upper limit curve is defined as
a condition when particles are triggered to move massively. In our problem, massive movement
of particles is not needed in modeling particle migration. Therefore, the lower limit of the movability
function
B
∗
.
$
is taken into the calculation, as:
G
∗
H
"
=
).I;
JK
∗
+12𝑒
:L.;JK
∗
+0.01 𝑙𝑛𝑅𝑒
∗
+0.078 (Eq. 4.3, [72])
where 𝑅𝑒
∗
= 𝑢
∗
"
CD
E
, and 0.1 < 𝑅𝑒
∗
< 10
F
.
𝑅𝑒
∗
is defined as grain Reynolds number, 𝑢
∗
"
is critical shear velocity to dislodge a particle,
𝐺𝑆 is grain size in diameter, 𝑣 is kinematic viscosity of fluid and 𝑤
/
refers to settling velocity of a
particle. According to the definition,
𝑣 =
4
N
(Eq. 4.4)
where, 𝜇 is dynamic viscosity of fluid, and 𝜌 refers to fluid density.
As the simulations are tested for various formation depths, temperature and viscosity values
are shown in Table 4.1, assuming the surface temperature is 65℉, and the geothermal gradient
is 1.8 ℉/100 𝑓𝑡.

Table 4.1: Water Viscosity And Density.
Depth,
𝑓𝑡.
Formation
Temperature,
℉
Dynamic
Viscosity, 𝑐𝑃
Density,
𝑙𝑏/𝑓𝑡.
2

Kinematic
Viscosity, 𝑓𝑡.
G
/
𝑠
1000 83 0.8259 62.187 8.924E-6
2000 101 0.6748 61.981 7.316E-6
3000 119 0.5632 61.728 6.132E-6
4000 137 0.4779 61.432 5.227E-6
5000 155 0.4119 61.099 4.530E-6
50
In calculating settling velocity 𝑤
/
in Eq. 4.3, relationships proposed by van Rijn (1993) [73]
for quartz grains are used. The relationships (Eq. 4.5 - Eq. 4.7) are empirically classified based
on grain diameters.
For particle diameter ranging between 1 to 100 microns, the settling velocity 𝑤
/
is,
𝑤
9
=
(N
H
:N)O∗PQ
?
,* RN
(Eq. 4.5, [73])
where 𝐺𝑆 is grain diameter. 𝜌
I
is grain density. For quartz, 𝜌
I
= 165.4
JK
@L.
%
(equivalent to
2.65 ×10
2
;I
N
%
). 𝑔 is gravitational acceleration, assumed to be 32.16
@L
/
&
(equivalent to 9.801
N
/
&
).
Settling velocity 𝑤
/
, for grains falling in the size range between 100 to 1000 microns, is,
𝑤
9
=
,) R
PQ
C1+
).), (N
H
:N)O∗PQ
(
R
?
D
).;
−1 (Eq. 4.6, [73])
For grains larger than 1000 microns in diameter,
𝑤
9
= 1.1*
SN
H
:NTO∗PQ
N
,
).;
(Eq. 4.7, [73])
The settling velocity with respect to grain diameters are plotted in Figure 4.2.

Figure 4.2: Settling Velocity 𝑤
/
Calculated With Van Rijn (1993) Correlations.
51
4.1.4: Grain size distribution and grain population density
Formation grains may or may not migrate, depending on whether the apparent flow velocity
is high enough to dislodge them. Mobilizing finer particles requires lower apparent velocities. To
determine which grains can be mobilized, one may consult grain size distribution (GSD) measured.
In this study, the grain size distribution of unconsolidated formation [16] is taken.
In a grain size distribution, grain sizes in diameter are 𝐺𝑆
9OP
9OQ
'
, and weight percentages
corresponding to the grain sizes are 𝐺𝑊
9OP
9OQ
'
. Given a reservoir condition, apparent flow velocity
𝑉
@
(𝑥,𝑡) can be calculated and compared with the critical shear velocity 𝑢
∗
(𝐺𝑆). If 𝑉
@
(𝑥,𝑡) is higher
than 𝑢
∗
(𝐺𝑆), grains with diameter 𝐺𝑆 are movable. Similarly, sizes of all movable grains at each
location and time can be found, using Eq. 4.3 - Eq. 4.7.
Unconsolidated formations contain fine particles smaller than 1 micron [47], which are
movable under most of the waterflooding conditions. Also, it is not unusual that fine particle can
weigh a considerable percentage of the dry bulk [29]. Therefore, migration of fine particles and its
effect on porosity and permeability change cannot be neglected.
However, many measured GSD’s do not show good resolutions at sizes smaller than 1
micron [16]. Another issue is that data points are too sparse on GSD’s. The reasons for low data
density of GSD could be coarse measurements in sieve analysis, or unrepresentative sampling.
It has been realized that sampling from produced sands tends to show a smaller-sized distribution,
while core sampling around the wellbore after certain particle migration has happened, tends to
result in a bigger-sized distribution. Sidewall coring of formation before its initial production is
recommended. Also, one has to target the most representative location of the formation.
With coarsely measured GSD’s, simulation results of grain movement display sharp cuts on
the finer end of new GSD’s. This is against lab and field observations, where grains of various
sizes are being produced. In dealing with the situation, first preference is to try to access the
measured initial GSD with a higher resolution. If not possible, one could extrapolate the GSD
52
curve toward the finer-sized end. In addition, through curve fitting, one may interpolate the curve
to obtain a higher data density on the GSD.
When extrapolation and interpolation are needed, estimation of a proper data density is
required to balance simulation accuracy and computational cost. Inverse modeling of the process
by comparing calculated average permeability increase to measured profile survey may provide
some basic guidance on severity of grain migration and sand erosion. If sand production has
occurred, the best is to measure the sizes of produced sands to determine the data density
needed in refining the original coarse GSD. If sand production has not happened, or the data is
not available, one could start with extracting one point for every 1% to 2% of incremental weight
on the cumulative GSD curve. Our experience is that for an open-hole injection well, operated at
typical water injection pressure, at least 10 data points representing grains smaller than 50
microns are required.
For a grain size distribution that includes particles of 𝑛
I
sizes, their sizes in diameter are
𝑑
9OP
9OQ
'
. With a spherical geometry assumed, grain volumes 𝐺𝑉
9OP
9OQ
'
are,
𝐺𝑉
UV,
UV6
H
=
W
I
*𝐺𝑆
UV,
UV6
H
,
=
(Eq. 4.8)
Weight percentages, 𝐺𝑊
9OP
9OQ
'
, can be obtained from either cumulative or incremental grain
size distribution. Assume the measured dry sample has a bulk mass of 𝑚 and grain density 𝜌
/
.
Grain populations 𝐺𝑃
9OP
9OQ
'
, representing the count for particles of 𝑛
I
sizes, are
𝐺𝑃
UV,
UV6
H
=
PX
RST
RSU
H
YN
"
PZ
RST
RSU
H
=
I
W
PX
RST
RSU
H
YN
"
[PQ
RST
RSU
H
\
(
(Eq. 4.9)
Then, dimensionless grain population densities 𝐺𝑃𝐷
9OP
9OQ
'
, become
53
𝐺𝑃𝐷
UV,
UV6
H
=
P1
RST
RSU
H
∑P1
RST
RSU
H
=
VW
RST
RSU
H
XVY
RST
RSU
H
Z
(
∑
VW
RST
RSU
H
XVY
RST
RSU
H
Z
(
(Eq. 4.10)

4.1.5: Pore size variation
With the onset of sand erosion and particle deposition caused by particle movement, not only
grain size distribution change, pore size distribution is no longer at its initial condition. As fine
particles are initially suspended in pore spaces, or attached to pore walls, movement of the
particles results in pore size changes. If a moving particle partially or fully clog a pore throat of a
smaller-sized pore and deposits there, the pore throat may shrink or get reduced to zero in size.
If fluid flow erodes a fine particle off the pore wall, the pore throat size may get enlarged. In real
situations, particle deposition and erosion could happen simultaneously due to the size exclusion
nature of porous media. Various scenarios are defined in Section 4.2, for using Markov Chain
Probability Transition process to model pore size variations.
A pore size distribution has 𝑛
[
sizes initially. In representing the case of pore size
enlargement and clogging, two additional sizes are introduced. Size of zero represents fully
clogged pores. The other size introduced is greater than the largest pore on the pore size
distribution, representing particle loss induced pore enlargement. In total, the simulation carries
𝑛
[
+2 sizes in accounting for pore size variations.
At the beginning of a timestep, the vector of pore size distribution is [𝑃
LP
]
P×]Q
(
^G_
. When some
pores get bridged while the others get enlarged, net change in pore size distribution becomes
[𝑃
LG
]
P×]Q
(
^G_
. Here the subscripts 𝑡1, 𝑡2 of 𝑃 refer to current and future state, respectively. The
change in pore size distribution in reflected in probability transition matrix, [T]
]Q
(
^G_×]Q
(
^G_
. Each
54
entry vector [𝑃] is count percentage of a pore size, which has to be less or equal than 1. The
summation of all entries in vector [𝑃] is 1.
[𝑃
CL
]
1×S𝑛
𝑝
+2T
= [𝑃
C,
]
1×S𝑛
𝑝
+2T
× [T]
S6
d
8LT×S6
d
8LT
(Eq. 4.11)
One could name the number of realizations in a timestep, along the stochastic process. If
two realizations are generated in every timestep as shown in Eq. 4.12 and Eq. 4.13, for each of
the two paths, two more realizations will be generated in the next timestep. Therefore, by the end
of the simulation of 𝑛
L
timesteps, 2
Q
*
realizations are generated, as shown in Figure 4.3.
J𝑃
CL,,
K
, ×S6
d
8LT
= [𝑃
C,
]
, ×S6
d
8LT
×[T
,
]
S6
d
8LT×S6
d
8LT
(Eq. 4.12)
J𝑃
CL,L
K
, ×S6
d
8LT
= [𝑃
C,
]
, ×S6
d
8LT
×[T
,
]
S6
d
8LT×S6
d
8LT
(Eq. 4.13)
The pore size transition matrices [T
P
]
]Q
(
^G_×]Q
(
^G_
and [T
G
]
]Q
(
^G_×]Q
(
^G_
are selected to
honor the cases of pores remaining unchanged, plugged, enlarged, or a combination of those
scenarios, following the definitions listed in Section 4.2.

Figure 4.3: Number Of Realizations Generated In Stochastic Process.

Note that theoretically, more than two realizations are computable per timestep. For instance,
if 3 realizations are generated with 3 matrices, by the end of a simulation discretized for 𝑛
L
steps,
Initial Condition
!
!","
t = 0
t = Δt
t = 2Δt
t = n*Δt
!
!",$
!
!$,"
!
!$,$
!
!$,%
!
!$,&
in total 2
'
!
realizations
⋯
⋯⋯⋯⋯
55
3
Q
*
realizations will be generated. In the study, two realizations are generated per step because
that, as computation proceeds to the end, the number of realizations generated is enormous.
Among the realizations, many are quite close to others and do not increase the diversity much.
For the computational example, generating more realizations per timestep has no mathematical
significance.

4.1.6: Particle movement and migration
For a particle to move, all of the following conditions must be met:
1. Apparent velocity 𝑉
@
(𝑥,𝑡) brought by injected fluid must exceed the motion threshold 𝑢
∗
(𝑑)
of the particles. Given a grain size distribution, the critical shear velocity 𝑢
∗
is calculated for all
sizes 𝐺𝑆
9OP
9OQ
'
(𝑥,𝑡).
2. Pore sizes 𝑃𝑆
9OP
9OQ
(
^G
(𝑥,𝑡) with their distribution [𝑃
LP
]
P×]Q
(
^G_
(𝑥,𝑡) at the beginning of
current timestep, must be greater than particle sizes 𝐺𝑆
9OP
9OQ
'
(𝑥,𝑡) so that the particles can pass
through. If not, the bigger particles will stay in-situ.
3. Grain population density 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡) at current step and location is more than critical
residual fraction of its initial population density 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡 = 0). Here 𝑃𝐿𝑅 is assumed to be 0.2.
The value could be calibrated by lab measurements. The significance of defining the factor 𝑃𝐿𝑅
honors the lab observations that, there will always be some fine particles left in the porous media,
even after a long term of fluid flush.
If all of the above conditions are met, probabilities 𝑃𝑟𝑜𝑏
9OP
9OQ
'
(𝑥,𝑡) for particles of sizes 𝐺𝑆
9OP
9OQ
'

to pass through the pores are,
56
𝑃𝑟𝑜𝑏
!"#
!"$
!
(𝑥,𝑡) =
⎩
⎪
⎨
⎪
⎧
𝑃𝐿𝐹∗2[𝑃
%#
(𝑥,𝑡)]
&"#
&" $
"
, 𝑖𝑓
𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡)
𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡 =0)
> 𝑃𝐿𝑅
0, 𝑖𝑓
𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡)
𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡 =0)
£ 𝑃𝐿𝑅

(Eq. 4.14)
for all j that 𝑃𝑆
:
> 𝐺𝑆
9OP
9OQ
'
. 𝑃𝐿𝐹 is defined as particle loss factor per timestep, which needs to
be calibrated with lab measurements for each formation. Here, 𝑃𝐿𝐹 is assumed to be 0.5.
Traveling distance, 𝐷𝑖𝑠𝑡
9OP
9OQ
'
(𝑥,𝑡) of the moving particles becomes,
𝐷𝑖𝑠𝑡
UV,
UV6
H
(𝑥,𝑡) = N
Δt∗𝑉
F
(𝑥,𝑡)
Δx
, 𝑖𝑓 𝑎𝑙𝑙 3 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛𝑠 𝑎𝑟𝑒 𝑡𝑟𝑢𝑒
0, 𝑖𝑓 𝑎𝑛𝑦 𝑜𝑓 𝑡ℎ𝑒 3 𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛𝑠 𝑖𝑠 𝑓𝑎𝑙𝑠𝑒

(Eq. 4.15)
where, Δt is timestep size, and Δx is the grid block size. 𝐷𝑖𝑠𝑡
9OP
9OQ
'
(𝑥,𝑡) calculates the number
of grid blocks that the particles have travelled the farthest in current step 𝑡.

4.1.7: Update in particle population density, porosity, and permeability
As particles of sizes 𝐺𝑆
9OP
9OQ
'
(𝑥,𝑡) migrate to the neighboring and more distant locations by
end of step 𝑡, new particle distribution 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡
L^P
) prior to the next step 𝑡+1, becomes,
𝐺𝑃𝐷
!"#
!"$
!
= 𝑥+𝑥
&
,𝑡
%( #
?=𝐺𝑃𝐷
!"#
!"$
!
= 𝑥+𝑥
&
,𝑡
%
?+𝑓
(*
#
+*)
∗𝑃𝑟𝑜𝑏
!"#
!"$
!
(𝑥,𝑡)∗𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡
%
)
(Eq. 4.16)

𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡
%( #
) =𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡
%
)−𝑃𝑟𝑜𝑏
!"#
!"$
!
(𝑥,𝑡) ∗𝐺𝑃𝐷
!"#
!"$
!
(𝑥,𝑡
%
)
(Eq. 4.17)
57
For all 𝑥
:
£ 𝐷𝑖𝑠𝑡
9OP
9OQ
'
(𝑥,𝑡). 𝑓 refers to filtration factor, which is found iteratively such that
∑ 𝑓
(f
+
gf)
f
+
O&9/L
!,-
!,.
'
(f,L)
f
+
OP
= 1 to maintain material balance.
Original formation porosity 𝛷 is the void fraction in bulk,
𝛷(𝑥,𝑡
C
) = 1−
W
IZ
∑𝐺𝑃𝐷
UV,
UV6
H
(𝑥,𝑡
C
)∗_𝐺𝑆
UV,
UV6
H
(𝑥,𝑡)`
=
−𝑎 (Eq. 4.18)
where a is non-movable grain volume fraction, calculated similarly as Eq. 4.8 – Eq. 4.10.
By rearranging Eq. 4.18, movable particle volume fraction, b, is,
𝑏(𝑥,𝑡) =
W
IZ
∑𝐺𝑃𝐷
UV,
UV6
H
(𝑥,𝑡
)
)∗_𝐺𝑆
UV,
UV6
H
(𝑥,𝑡)`
=
= 1 – a - 𝛷(𝑥,𝑡
)
) (Eq. 4.19)
At new particle distribution 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡
L^P
), formation porosity 𝛷(𝑥,𝑡
L^P
) is updated to be,
𝛷(𝑥,𝑡
C8,
) = 𝛷(𝑥,𝑡
)
)+𝑏∗

j
kl
∑P1"
RST
RSU
H
(A,C
>mT
)∗[PQ
RST
RSU
H
\
(
j
kl
∑P1"
RST
RSU
H
(A,C
5
)∗[PQ
RST
RSU
H
\
(
(Eq. 4.20)
𝛷(𝑥,𝑡
C8,
) = 𝛷(𝑥,𝑡
)
)+(1−𝛷(𝑥,𝑡
)
)−𝑎 )∗(1−
∑P1"
RST
RSU
H
(A,C
>mT
)∗[PQ
RST
RSU
H
\
(
∑P1"
RST
RSU
H
(A,C
5
)∗[PQ
RST
RSU
H
\
(
)
(Eq. 4.21)
Here the subscript 0 of 𝑡 refers to initial condition.
Using the Kozeny-Carman relationship [74], permeability 𝑘(𝑥,𝑡
L^P
) is updated to be,
𝑘(𝑥,𝑡
C8,
) =
^
?
5
H
?
D(A,C
>mT
)
(
?(, :D(A,C
>mT
))
?
(Eq. 4.22)
where 𝑙
I
is mean grain size, 𝜓 is sphericity of sand particles, and A is a dimensionless
coefficient dependent upon the tortuosity of porous media.

58
4.2: Stochastic modeling using the Markov Chain Probability
Transition process
The Markov Chain process is applied in a wide variety of industries wherein future state only
depends on current state and probability transition matrix, when a series of states are during
evolution. This memory-lessness is a key Markov property. If a process satisfies the Markov
property, the Markov Chain probability transition process can be used.
[𝑃
CL
]
, ×6
= [𝑃
C,
]
, ×6
×[T]
6×6
(Eq. 4.23)
where the subscript 𝑛 refers to the total number of possibilities involved in the variations, 𝑡
P

and 𝑡
G
refer to current and future time. The probability transition matrix [T] should always be a
square matrix.
For the ease of demonstration, we set a convention for pore sizes 𝑃𝑆
:
, such that
𝑃𝑆
68,
> 𝑃𝑆
6
> ⋯ > 𝑃𝑆
L
> 𝑃𝑆
,
> 𝑃𝑆
_
= 0 (Eq. 4.24)
During the continuing particle movement process, pores could continuously shrink or enlarge
in size or stay the same within a certain timestep. Small pores could get completely bridged. Each
of the scenarios is discussed below.
Assume reservoir pore size distribution is [𝑃
LP
], when subject to waterflooding induced
particle movement, the change in pore size distribution is represented by transition matrix [𝑇],
and new pore size distribution becomes [𝑃
LG
]. As discussed in Section 4.1.5, with 2 more sizes
added, the process is shown as,
[𝑃
CL
]
, ×S6
d
8LT
= [𝑃
C,
]
, ×S6
d
8LT
×[T]
S6
d
8LT×S6
d
8LT
(Eq. 4.25)
a. If original reservoir pores are uniform in size 𝑃𝑆
P
, the vector of initial pore size distribution
[𝑃
LP
] is,
59
[P
LP
] =
⎣
⎢
⎢
⎢
⎢
⎡
P
LP,nD
/
P
LP,nD
-
P
LP,nD
&

⋯
P
LP,nD
.
P
LP,nD
.0-
⎦
⎥
⎥
⎥
⎥
⎤
=
⎣
⎢
⎢
⎢
⎢
⎡
0
1
0
⋯
0
0
⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.26)
Note that although pores may reach the diameter of 0 or 𝑃𝑆
Q^P
as particle deposition and
erosion happens, initially, the percentage of those two conditions is 0. Also, for all transition
matrices shown in the dissertation, future states refer to the columns while current states refer to
the rows.
The probability transition matrix [𝑇] is,
𝑃𝑆
o
= 0 𝑃𝑆
P
𝑃𝑆
G
⋯ 𝑃𝑆
Q
𝑃𝑆
Q^P

𝑃𝑆
o
= 0
𝑃𝑆
P
𝑃𝑆
G

⋯
𝑃𝑆
Q
𝑃𝑆
Q^P

⎣
⎢
⎢
⎢
⎢
⎡
𝑃
7,7
𝑃
P,7
𝑃
G,7
𝑃
7,P
𝑃
P,P
𝑃
G,P
⋯
𝑃
Q,7
𝑃
Q^P,7
⋯
𝑃
Q,P
𝑃
Q^P,P

𝑃
7,G
𝑃
P,G
𝑃
G,G
⋯
⋯
⋯
⋯
𝑃
Q,G
𝑃
Q^P,G
⋯
⋯
⋯

𝑃
7,Q
𝑃
P,Q
𝑃
G,Q
𝑃
7,Q^P
𝑃
P,Q^P
𝑃
G,Q^P
⋯
𝑃
Q,Q
𝑃
Q^P,Q
⋯
𝑃
Q,Q^P
𝑃
Q^P,Q^P⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.27)
where ∑ 𝑃
9,:
= 1
:OQ^P
:O7
, for all 𝑖 = 0,1,2,…,𝑛,𝑛+1.

b. If original formation pores have n sizes ranging from 𝑃𝑆
P
to 𝑃𝑆
Q
, initial pore size distribution
[𝑃
LP
] becomes,
[P
LP
] =
⎣
⎢
⎢
⎢
⎢
⎡
P
LP,nD
/
P
LP,nD
-
P
LP,nD
&

⋯
P
LP,nD
.
P
LP,nD
.0-
⎦
⎥
⎥
⎥
⎥
⎤
=
⎣
⎢
⎢
⎢
⎢
⎡
0
P
LP,P
P
LP,G

⋯
P
LP,Q
0
⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.28)
where ∑ 𝑃
LP,:
= 1
:OQ^P
:OP
.

c. If all pores do not change in size, initial pore size distribution probability [P
LP
] is the same
as in a or b, the transition matrix [𝑇], is:

Future
Current
60
𝑃𝑆
o
= 0 𝑃𝑆
P
𝑃𝑆
G
⋯ 𝑃𝑆
Q
𝑃𝑆
Q^P

𝑃𝑆
o
= 0
𝑃𝑆
P
𝑃𝑆
G

⋯
𝑃𝑆
Q
𝑃𝑆
Q^P

⎣
⎢
⎢
⎢
⎡
1
0
0
0
1
0
⋯
0
0
⋯
0
0

0
0
1
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯

0
0
0
0
0
0
⋯
1
0
⋯
0
1
⎦
⎥
⎥
⎥
⎤
(Eq. 4.29)

d. If all pores shrink in size as partially or fully clogged by particles, initial pore size distribution,
transition matrix [𝑇] is,
𝑃𝑆
o
= 0 𝑃𝑆
P
𝑃𝑆
G
⋯ 𝑃𝑆
Q
𝑃𝑆
Q^P

𝑃𝑆
o
= 0
𝑃𝑆
P
𝑃𝑆
G

⋯
𝑃𝑆
Q
𝑃𝑆
Q^P

⎣
⎢
⎢
⎢
⎢
⎡
𝑃
7,7
𝑃
P,7
𝑃
G,7
0
𝑃
P,P
𝑃
G,P
⋯
𝑃
Q,7
𝑃
Q^P,7
⋯
𝑃
Q,P
𝑃
Q^P,P

0
0
𝑃
G,G
⋯
⋯
⋯
⋯
𝑃
Q,G
𝑃
Q^P,G
⋯
⋯
⋯

0
0
0
0
0
0
⋯
𝑃
Q,Q
𝑃
Q^P,Q
⋯
0
𝑃
Q^P,Q^P
⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.30)
where ∑ 𝑃
9,:
= 1
:O9
:O7
, for 𝑖 = 0,1,2,…,𝑛,𝑛+1.

e. If all pores get enlarged in size as the particles passing through the pore throats and create
more voids, transition matrix [𝑇] is
𝑃𝑆
o
= 0 𝑃𝑆
P
𝑃𝑆
G
⋯ 𝑃𝑆
Q
𝑃𝑆
Q^P

𝑃𝑆
o
= 0
𝑃𝑆
P
𝑃𝑆
G

⋯
𝑃𝑆
Q
𝑃𝑆
Q^P

⎣
⎢
⎢
⎢
⎢
⎡
𝑃
7,7
0
0
𝑃
7,P
𝑃
P,P
0
⋯
0
0
⋯
0
0

𝑃
7,G
𝑃
P,G
𝑃
G,G
⋯
⋯
⋯
⋯
0
0
⋯
⋯
⋯

𝑃
7,Q
𝑃
P,Q
𝑃
G,Q
𝑃
7,Q^P
𝑃
P,Q^P
𝑃
G,Q^P
⋯
𝑃
Q,Q
0
⋯
𝑃
Q,Q^P
𝑃
Q^P,Q^P⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.31)
where ∑ 𝑃
9,:
= 1
:OQ^P
:O9
, for 𝑖 = 0,1,2,…,𝑛,𝑛 + 1.
f. In real cases, pore sizes could experience changes of all possibilities. In another word,
among the current pores with n sizes, some may shrink, some enlarge in size, and some get
completed plugged so their sizes become 0. Transition matrix [𝑇] is,

61
𝑃𝑆
o
= 0 𝑃𝑆
P
𝑃𝑆
G
⋯ 𝑃𝑆
Q
𝑃𝑆
Q^P

𝑃𝑆
o
= 0
𝑃𝑆
P
𝑃𝑆
G

⋯
𝑃𝑆
Q
𝑃𝑆
Q^P

⎣
⎢
⎢
⎢
⎢
⎡
𝑃
7,7
𝑃
P,7
𝑃
G,7
𝑃
7,P
𝑃
P,P
𝑃
G,P
⋯
𝑃
Q,7
𝑃
Q^P,7
⋯
𝑃
Q,P
𝑃
Q^P,P

𝑃
7,G
𝑃
P,G
𝑃
G,G
⋯
⋯
⋯
⋯
𝑃
Q,G
𝑃
Q^P,G
⋯
⋯
⋯

𝑃
7,Q
𝑃
P,Q
𝑃
G,Q
𝑃
7,Q^P
𝑃
P,Q^P
𝑃
G,Q^P
⋯
𝑃
Q,Q
𝑃
Q^P,Q
⋯
𝑃
Q,Q^P
𝑃
Q^P,Q^P⎦
⎥
⎥
⎥
⎥
⎤
(Eq. 4.32)
where ∑ 𝑃
9,:
= 1
:OQ^P
:O7
, for all 𝑖 = 0,1,2,…,𝑛,𝑛+1.

4.3: Inverse modeling
As grains dislodge and leave the original location, permeability increase 𝐾
7
(𝑥,𝑡
L^P
) is
𝐾
)
(𝑥,𝑡
C8,
) =
. (A,C
>mT
)
. (A,C
>
)
=
D(A,C
>mT
)
(
S, :D(A,C)T
?
D(A,C
>
)
(
S, :D(A,C
>mT
)T
?
(Eq. 4.33)
At the end of the simulation, average permeability 𝐾
pE-
7
(𝑥,𝑡 = 𝑛
L
) of all locations can be used
to match the increased injectivity seen from a profile survey, or history match the waterflooding
performance. 𝐾
pE-
7
(𝑥,𝑡 = 𝑛
L
) is calculated as in Eq. 4.34.
,
`
qr+
5
(A,CV6
>
)
= ∑
,
`
5
(A,CV6
>
)
AV6
s
AV,
(Eq. 4.34)
where 𝑛
f
=
t
1
uv
, 𝑛
L
=
L
*2*34
uw
.
If none of the realizations generated matches the profile survey, a refinement of 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡),
or reselection of the two Probability Transition Matrices [T
P
]
]Q
(
^G_∗]Q
(
^G_
and [T
G
]
]Q
(
^G_∗]Q
(
^G_
are
needed.
One may update 𝐺𝑆
9OP
9OQ
'
(𝑥,𝑡) by increase the value of 𝑛
I
for 1 or more. Suppose one particle
size is added per loop,
𝑛
O
(𝑙𝑜𝑜𝑝 = 𝑙+1) = 𝑛
O
(𝑙𝑜𝑜𝑝 = 𝑙)+1 (Eq. 4.35)
62
Transition matrices can be re-selected by increasing the chances of bigger pores getting
enlarged for 10%, if the average permeability increase 𝐾
pE-
7
(𝑥,𝑡 = 𝑛
L
) is smaller than the
injectivity increase shown on profile survey. Similarly, one can increase the chances of plugging
the smaller pores in the transition matrix for 10%, if the simulated sand erosion happens faster
than what is shown on profile surveys. In both situations, the sum of all possibility entries in future
states of a single pore size needs to be 1, meaning the summation of all entries in any single row
in the matrix [𝑇] is 1.
From our experience, reselection of the matrices is required for no more than several times.
Most of the computation does not need to reselect the matrices as long as the grain size
distribution has good resolution especially at the finer-sized end. With the situation, the enormous
number of realizations generated presents a good diversity. The high demand is actually on the
resolution of grain size distribution. If lab measurements on 𝑃𝐿𝐹 and 𝑃𝐿𝑅 are missing, reselection
of them can be done similarly in order to matching the profile survey or field history.

4.4: Algorithm
The algorithm is shown schematically in Figure 4.4. Step 1 calculates pressure distribution
with Eq 4.1, and apparent flow velocity using Eq. 4.2. Step 2 calculates critical shear velocity for
particles with various sizes, using Eq. 4.3 – Eq. 4.7. The Markov Chain process brings two
different transition matrices in step A following Eq. 4. 12, Eq. 4.13. The two matrices are randomly
generated for twice using Eq. 4.32. Resulted pore size distributions [𝑃
LG
]
P×]Q
(
^G_
work as inputs
in Step 3, to determine the movability of grains. Calculations in Step 3 involve Eq. 4.14 – Eq. 4.17.
Step 4 and 5 calculate new grain size distribution 𝐺𝑃𝐷
9OP
9OQ
'
(𝑥,𝑡
L^P
), and resultant updates in
porosity 𝛷(𝑥,𝑡
L^P
) and permeability 𝑘(𝑥,𝑡
L^P
) using Eq. 4.18 – Eq. 4.22. As the computation
continues to the end with all realizations generated, Step 6 calculates the harmonically averaged
permeability increase 𝐾
pE-
7
(𝑥,𝑡 = 𝑛
L
) of the eroded formation, for matching the profile survey or
63
field data. The average permeability increase is calculated with Eq. 4.33 and Eq. 4.34. If matched,
the simulation accepts the transition matrix and the matched realizations. If not, the algorithm re-
selects one among the factors 𝑃𝐿𝐹, 𝑃𝐿𝑅, and transition matrices [𝑇
P
] and [𝑇
G
]. In cases where
the coarse GSD is the source of problem, refinements of GSD are needed. With new input, the
simulation repeats from Step 1 to 6 until 𝐾
pE-
7
(𝑥,𝑡 = 𝑛
L
) calculated in the realizations match the
permeability increase shown on profile surveys.
The accepted realizations will be possible reflections of how a thief zone is being developed,
and how petrophysical properties are changing dynamically. The dynamic porosity 𝛷(𝑥,𝑡
L^P
) and
permeability 𝑘(𝑥,𝑡
L^P
) can be incorporated in history matching the performance of an
unconsolidated reservoir under waterflooding.

64

Figure 4.4: Workflow And Inverse Modeling Of The Particle Migration Algorithm.

1. Calculate pressure distribution
!(#,%) and apparent fluid flowing
velocity '
!
(#,%).
A. New pore size distributions !
"# $× &
!
'#
after
probability transition using Markov Chain Stochastic
Process !
"# $× &
!
'#
= !
"$ $× &
!
'#
×[T]
&
!
'# × &
!
'#
.
2. Calculate critical shear velocity
u*(-) using Paphitis (2001).
3. Calculate travelling distance
-./%(0
(
,#,%) from,
a. Apparent flowing velocity
compared to critical shear
velocity;
b. PSD probability transition
and size exclusion;
c. Residual particle fraction
compared to PLR assumed.
4. new particle distribution
1!2
)*$
)*&
"
#,%
"'$
as a net result
of sand erosion and accumulation.
5. Update of Porosity 4 (#,%)
and Permeability 5 #,% .
6. Check Injectivity Increase
7
+,-
.
(#,% = 0
"
) with Profile
Surveys. Re-electing PLF,
PLR, transition matrices and
refinement of GSD might be
needed.
65
Validation approach

This chapter summarizes the parameterizations applied to validation of the algorithm.
Table 5.1 includes the general reservoir properties used in the simulation. The properties are
selected to represent a high porosity, high permeability unconsolidated sand body. In the base
case, a slab reservoir is located at the depth of 2000 ft., with a hydrostatic condition assumed for
initial reservoir pressure. An injection well and a producing well are located at the two ends of the
reservoir. The schematic is shown as in Figure 5.1. Both of the wellbores are vertical and of 6-
inch diameter. The injected fluid is fresh water at formation temperature. Table 5.2 includes
parameter variations tested in sensitivity studies. The values in parenthesis are used in the base
case. The algorithm was tested for formation depths ranging from 1000 to 5000 ft. At each depth,
water injection was tested at various pressure gradients, from 0.6 to 0.85 psi/ft. Formation
thicknesses were assigned to be 30 and 50 ft., besides the 10 ft. thickness assumed in the base
case. Injector-producer distances are assumed to range from 50 to 300 ft.
Figure 5.2 presents the cumulative grain size distribution of a typical California
unconsolidated sand [16], which is used in the validation. During the computation process, it is
realized that the GSD measured above 90 weight percentiles is not accurate enough which
resulted in trivial results. As the smaller-sized particles migrate the most under typical injection
pressure in waterflooding, extrapolation of the curve on the finer end and interpolation to add
more data density were done on the GSD, as shown in Figure 5.2. The refined incremental and
cumulative grain size distributions are shown in Error! Reference source not found. and Figure
5.3, respectively. Table 5.3 includes the exact values of extrapolated and interpolated data points
on the GSD. The assumed pore size distribution is shown in Figure 5.4. At initial condition, both
GSD and PSD are homogeneously applied to the entire reservoir. Figure 5.5 shows the computed
66
critical shear velocities to mobilize formation particles, using the lower limit curve of Paphitis (2001)
[72] movability function.

Table 5.1: Reservoir Properties For Base Case.
Reservoir Property Value Unit
Reservoir Size, 𝐿
f
∗ 𝐿
x
∗ 𝐿
y
50 * 10 * 10 ft. * ft. * ft.
Grid Size, 𝛥𝑥 ∗ 𝛥𝑦 ∗ 𝛥𝑧 5 * 10 * 10 ft. * ft. * ft.
Simulation Time, 𝑡
LoLpJ
300 days
Timestep, 𝛥𝑥 30 days
Initial Porosity, 𝛷
9Q9
0.3 dimensionless
Initial Permeability, 𝑘
9Q9

2000 mD
Total Compressibility, 𝐶
L

2.5 E-5 psi
-1

Initial Reservoir Pressure, 𝑃
9Q9
1299 psi
Injected Fluid Viscosity, 𝜇
0.6748 cP
Reservoir Thickness, ℎ = 𝛥𝑧 10 ft.
Formation Depth, 𝐷 3000 ft.
Wellbore radius, 𝑟
.
0.5 ft.
Injection gradient, 𝐼𝑛𝑗
Izp3
0.75 psi/ft.
Wellbore Completion Open-Hole -

Table 5.2: Sensitivity Parameterizations.
Reservoir Property Value Unit
Initial Permeability,
𝑘
9Q9

800, 1500, (2000,) 3000 mD
Injection Gradient,
𝐼𝑛𝑗
Izp3

0.6, 0.65, 0.7, (0.75,)
0.8, 0.85
psi/ft.
Reservoir Thickness,
ℎ = 𝛥𝑧
(10,) 30, 50 ft.
Formation Depth,
𝐷
1000, 2000, (3000,)
4000, 5000
ft.
Injector-Producer Distance, 𝐿
f
(50,) 100, 200, 300 ft.

Figure 5.1: Schematic Of The Slab Reservoir Assumed In Simulation.
Injector – Producer Distance, !
!
Producer
"
"#,%&'(
= 14.7 )*+
Injector
Thickness,ℎ
67

Figure 5.2: Cumulative Grain Size Distribution Of Sampled California Unconsolidated Sand.
([16] Courtesy SPE.)

Figure 5.3: Incremental Grain Size Distribution Of Sampled California Unconsolidated Sand:
Data Extracted From Figure 5.2 by Schwartz (1969). ([16] Courtesy SPE.)

0
20
40
60
80
100%
1.0
0.5
0.1
0.05
0.01
0.005
Cumulative Percentage By Weight
Grain Diameter, inches
Extrapolated
Data Points
0.001
Formation Sand
Profile
Gravel
Profile
C = 1.42
68

Figure 5.3: Initial Grain Size Distribution Extrapolated upon Figure 5.2.

Figure 5.4: The Assumed Initial Pore Size Distribution.

69
Table 5.3: Hand Extracted And Extrapolated Grain Size Distribution Values From Figure 5.2.
([16] Courtesy SPE.)
Diameter,
in.
Diameter,
micron
Individual
Weight
Percentage,
%
Cumulative
Weight
Percentage,
%
0.00071 18 0.1 100
0.00075 19 0.25 99.8
0.00083 21 0.4 99.5
0.00094 24 0.75 99
0.00110 28 1 98
0.00130 33 1 97
0.00142 36 1 96
0.00157 40 1 95
0.00169 43 1 94
0.00181 46 1 93
0.00193 49 1.5 92
0.00217 55 2.5 90
0.00248 63 2.75 87
0.00299 76 6 84.5
0.00402 102 9.75 75
0.00500 127 9.75 65
0.00598 152 9 55.5
0.00701 178 8.25 47
0.00799 203 8.25 39
0.00902 229 6.25 30.5
0.01000 254 7.35 26.5
0.01299 330 7.45 15.8
0.01500 381 5.05 11.6
0.02000 508 4.3 5.7
0.02500 635 4.35 3

70

Figure 5.5: Critical Shear Velocity At Various Depths And Formation Temperatures.

71
Results

This chapter presents the simulation results based on general properties and sensitivity test
conditions defined in Chapter 5.
In a simulation with 𝑡 timesteps, allowing 2 pathways for each step, the number of realizations
will be 2
L
by end of last step, as shown in Figure 4.3. In our simulations, ten timesteps are
computed, to balance the satisfactory resolution of dynamic presence of properties and
computational cost. As the particle movement proceeds over time, grain size distributions
generated in first three steps, among the substantial number of realizations (1024 realizations by
the end of the simulation), are presented, see Figure 6.1. It represents the base case, which
shows a 10 ft.-thick unconsolidated sand body with an initial permeability 2 Darcy, located at 3000
ft. depth, subject to a water injection gradient of 0.75 psi/ft. The injector-producer distance is 50
ft.
Starting from Figure 6.1(a), the initial condition of grain size distribution, two realizations are
generated with two distinct transition matrices [𝑇
P
] and [𝑇
G
] by end of the first timestep, as shown
in Figure 6.1(b). As the simulation continues, four realizations are generated by the second
timestep, as shown in Figure 6.1(c). As the simulation continues to the third timestep, in total 8
realizations are generated, as shown in Figure 6.1(d1) and Figure 6.1(d2). Comparing the varying
GSD’s with the initial state shown in Figure 6.1(a), the grains toward the finer (blue) end are
changing more drastically than the coarser grain (yellow) end. One may distinguish the fine
particles from the bar charts to observe the dynamics of GSD along the waterflooding process,
from Figure 6.1(a) to Figure 6.1(d1) and Figure 6.1(d2).

72

Figure 6.1(a): Initial GSD In Base Case.

Figure 6.1(b): Two Realizations Of GSD By The End Of The First Timestep In Base Case.

73

Figure 6.1(c): Four Realizations Of GSD By The End Of The Second Timestep In Base
Case.

Figure 6.1(d1): Four Realizations Of GSD By The End Of Third Timestep In Base Case.
74

Figure 6.1(d2): Another Four Realizations GSD By The End Of Third Timestep In Base
Case.
Figure 6.1: Initial Status, As Well As All Realizations For The First, Second, And The Third
Steps In Base Case Described in Table 5.1.

Figure 6.2 shows the evolution of porosity profiles with time. Starting from initial
homogeneous porosity distribution to the 5
th
timestep at 150
th
day, formation sand near the
wellbore of the injector gets rapidly eroded as a result of the high local pressure perturbation. In
the timesteps afterward, sand erosion substantially slows down, and the porosity change is getting
to the equilibrium. It supports the observations that, the amount of movable sand particles is finite
within a formation. Even when pore sizes are enlarged significantly, the bigger-sized particles
stay immobile as the critical shear velocity needed to trigger their movements cannot be reached
under typical injection pressures during waterflooding.
Figure 6.3 shows the changes in permeability profiles for the same case. Similar as porosity,
permeability increase is faster at the beginning of particle migration process. Then the increase
75
slows down and gradually reaches to the equilibrium. Curves in Figure 6.2 and Figure 6.3 are
plotted with averaged values of all 1024 realizations generated by the end of the simulation.
Figure 6.4 shows averaged permeability increase vs. time. For linear systems, permeability
is averaged harmonically using Eq. 4.34.

Figure 6.2: Averaged Porosity Evolution For The Base Case Described In Table 5.1.

Figure 6.3: Averaged Permeability Evolution For The Base Case Described In
Table 5.1.
76

Figure 6.4: Harmonically Averaged Permeability Increase Vs. Injection Duration For The Base
Case Described In Table 5.1.

As shown in the base case, particle loss and sand erosion override sand particle
accumulation, so that the overall permeability of the formation increases over time. Cases testing
the sensitivities of layer thickness, injection pressure gradient, formation depth and initial reservoir
permeability are computed. As mentioned in Chapter 4, profile surveys can validate the
realizations and minimize the number of possible cases. Here, to demonstrate the significance of
various factors, permeability increases of all 1024 realizations at the end of the simulation are
averaged and summarized in Table 6.1 – Table 6.12. The results are based on the slab reservoir
geometry shown in Figure 5.1. Both the injector and producer are completed as open holes, with
a distance of 50 ft.
Again, the sand migration and erosion issue are much more severe in perforated injection
wells, as the flowing area is significantly limited by perforations. It also holds true in the injection
wells around which the sand face is partially damaged, where effective thickness is apparently
smaller than the thickness of entire hydrocarbon-bearing interval. This will be discussed in
Chapter 7 in detail.
77
In Table 6.1 – Table 6.12, all the values are ratios of increased permeability over the original
permeability of the sand body. Value 1 means the permeability remains the same during
waterflooding. From the results, we made the following observations.
1. As formation depth increases, particle migration and sand erosion are more severe. This
is because that, as a product of formation depth and injection pressure gradient, injection pressure
becomes higher for a deeper formation. For the same reason, as injection gradient increases,
injection pressure is also higher. When outlet pressure at the producing well is kept as 14.7 psi,
higher injection pressure delivers a higher injection rate, which mobilizes the bigger-sized sand
particles more easily. When more sand particles travel far and eventually leave the formation,
more void space are left behind in the reservoir. Accordingly, permeability increases significantly
with the depth and pressure gradient, as noticed in each of the tables.
2. High initial permeability sand is more prone to particle migration and sand erosion. For
instance, comparing Table 6.1 and Table 6.4, keeping all other conditions the same, a 10 – ft.
thick formation of 3 Darcy initial permeability, located at 5000 ft. depth getting fluid injection under
a pressure gradient of 0.75 psi/ft, shows an average permeability increase being 4.814 times of
its original value (3 Darcy), while the permeability increase for a formation with an 800 mD initial
permeability is 2.6292 times of its initial value (800 mD). This is because high initial permeability
offers more effective pathways for both fluid and sand particles. The behavior is observed from
comparing within the group of Table 6.1 – Table 6.4, the group of Table 6.5 – Table 6.8, the group
of Table 6.9 – Table 6.12.
3. Assuming the sand face around an open-hole injection well is clean, meaning the entire
interval contributes to fluid flow, thickness does not affect sand erosion much. In real cases,
thickness of the interval seems to play an important role, which will be elaborated in Chapter 7.
4. It should be noted that as the data density on grain size distribution (GSD) and pore size
distribution (PSD) increases, the percentages of trivial results significantly decrease. This is
78
especially true for the cases where more data points are extracted toward the finer end of GSD
and PSD.
5. The values listed in tables are averaged among all 1024 realizations for each
computational case to reveal the effect of each factor on particle migration. However, when
examine certain individual realizations, the resulted permeability increase can be substantially
higher. Those extreme realizations are generated when the randomly assigned transition matrix
represents erosion more than sand accumulation. Therefore, selection of more distinct matrices
leads to more diversified realizations.

Table 6.1: Permeability Increase Ratio For 10 Ft. Thick Sand With 800 Md Initial Permeability.
Injection gradient,
psi/ft. 0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1
3000 1 1 1 1 1 1
4000 1 1.0242 1.053 1.1367 1.5921 2.4382
5000 3.1022 3.2337 3.2384 3.2384 3.2384 3.2384

Table 6.2: Permeability Increase Ratio For 10 Ft. Thick Sand With 1500 Md Initial Permeability.
Injection gradient,
psi/ft. 0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1.0070
3000 1.0862 1.1722 1.7504 2.4963 3.0662 3.1088
4000 3.1137 3.1137 3.1137 3.1137 3.1137 3.1137
5000 3.1137 3.1137 3.2060 3.3152 3.4025 3.5004

79
Table 6.3: Permeability Increase Ratio For 10 Ft. Thick Sand With 2000 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1.0071 1.0241 1.0528 1.1066 1.1465
3000 3.0107 3.1980 3.2014 3.2014 3.2014 3.2014
4000 3.2014 3.2014 3.2014 3.2014 3.3830 3.4758
5000 3.5086 3.7966 4.0320 4.2906 4.6245 5.0872

Table 6.4: Permeability Increase Ratio For 10 Ft. Thick Sand With 3000 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1.2362 1.9277 2.6103 3.1014 3.2252 3.2320
3000 3.2320 3.2320 3.2320 3.2320 3.2320 3.3618
4000 3.5767 3.8089 4.0904 4.4520 4.6718 5.1948
5000 5.6402 6.5352 7.1875 8.0571 9.8843 10.3858

Table 6.5: Permeability Increase Ratio For 30 Ft. Thick Sand With 800 Md Initial Permeability.
injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1
3000 1 1 1 1 1 1
4000 1 1.0239 1.0524 1.1350 1.5737 2.3809
5000 3.0370 3.1700 3.1743 3.1743 3.1743 3.1743

Table 6.6: Permeability Increase Ratio For 30 Ft. Thick Sand With 1500 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1.0071
3000 1.0874 1.1755 1.7916 2.5976 3.2115 3.2591
4000 3.2659 3.2659 3.2659 3.2659 3.2659 3.2659
5000 3.2659 3.2659 3.4420 3.5653 3.6477 3.8588

80
Table 6.7: Permeability Increase Ratio For 30 Ft. Thick Sand With 2000 Md Initial Permeability.
injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1.0071 1.0243 1.0532 1.1075 1.148
3000 3.0691 3.2592 3.2634 3.2634 3.2634 3.2634
4000 3.2634 3.2634 3.2634 3.2893 3.4542 3.5788
5000 3.6501 3.8691 4.1678 4.5715 4.8078 5.3617

Table 6.8: Permeability Increase Ratio For 30 Ft. Thick Sand With 3000 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1.2313 1.8769 2.5243 2.9922 3.1107 3.1176
3000 3.1176 3.1176 3.1176 3.1176 3.1176 3.2436
4000 3.3998 3.5686 3.8538 4.1115 4.4744 4.8493
5000 5.2560 6.0983 6.8722 7.7966 9.5164 9.9120

Table 6.9: Permeability Increase Ratio For 50 Ft. Thick Sand With 800 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1
3000 1 1 1 1 1 1
4000 1 1.0238 1.0522 1.1338 1.5679 2.3675
5000 3.0296 3.1609 3.1644 3.1644 3.1644 3.1644

Table 6.10: Permeability Increase Ratio For 50 Ft. Thick Sand With 1500 Md Initial Permeability.
injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1 1 1 1 1.0071
3000 1.0881 1.1766 1.7905 2.5793 3.1786 3.2256
4000 3.2318 3.2318 3.2318 3.2318 3.2318 3.2318
5000 3.2318 3.2318 3.3852 3.5155 3.5843 3.7509

81
Table 6.11: Permeability Increase Ratio For 50 Ft. Thick Sand With 2000 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1 1.007 1.024 1.0525 1.1058 1.1456
3000 3.0463 3.2358 3.2404 3.2404 3.2404 3.2404
4000 3.2404 3.2404 3.2404 3.2431 3.4282 3.5487
5000 3.5948 3.8394 4.1187 4.4466 4.7576 5.2770

Table 6.12: Permeability Increase Ratio For 50 Ft. Thick Sand With 3000 Md Initial Permeability.
Injection
gradient,
psi/ft.
0.6 0.65 0.7 0.75 0.8 0.85
Depth, ft.
1000 1 1 1 1 1 1
2000 1.236 1.9246 2.5989 3.0996 3.2123 3.2182
3000 3.2182 3.2182 3.2182 3.2182 3.2182 3.3492
4000 3.5261 3.8121 4.0691 4.4777 4.6573 5.2272
5000 5.698 6.6 7.4211 8.2767 10.4211 10.9843

Besides the sensitivity factors tested, it is noticed that pressure drawdown per lateral distance
un
uf
affects particle migration strongly. For instance, if increase the distance between the injector
and producer, the formation of base case may not show sand erosion, keeping the outlet pressure
at the producing well as 14.7 psi. In the base case, where injector-producer distance is 50 ft.,
obvious sand erosion and permeability increase are seen. If changing the distance to be 100 ft.,
porosity changes from its original value of 0.3 to only 0.3006, as shown in Figure 6.5(a). The
permeability change is also trivial, see Figure 6.6(a).
The sensitivity of injector-producer distance is also tested for 200 and 300 ft. Porosity profiles
are shown in Figure 6.5(b) and (c), where the porosity is calculated using Eq. 4.18 – Eq. 4.21.
The associated permeability profiles are shown in Figure 6.6(b) and (c), and the permeability
calculations follow Eq. 4.22 and Eq. 4.33. Figure 6.5 and Figure 6.6 show that, as the distance
between the producer and the open-hole injector increases further, sand erosion is not seen, as
82
indicated by the constant porosity and permeability profiles. The reason is that pressure change
per linear distance, instead of pressure drawdown, determines flow rate and velocity. The injector-
producer distance ∆𝑥 inversely affects pressure drawdown per distance
un
uf
. Given a large injector-
producer distance, the same injection gradient and pressure become insufficient to deliver the
same apparent flow velocity to trigger grain mobilization. Particles can be mobilized only when
apparent flow velocity exceeds their critical shear velocities.

Figure 6.5(a): Porosity Profiles When Injector-Producer Distance Is 100 Ft.

83

Figure 6.5(b): Porosity Profiles When Injector-Producer Distance Is 200 Ft.

Figure 6.5(c): Porosity Profiles When Injector-Producer Distance Is 300 Ft.
Figure 6.5: Porosity Profile Evolutions For The Base Case Mentioned In Table 5.1 Except For A
Change In Injector-Producer Distance.

84

Figure 6.6(a): Permeability Profiles When Injector-Producer Distance Is 100 Ft.

Figure 6.6(b): Permeability Profiles When Injector-Producer Distance Is 200 Ft.

85

Figure 6.6(c): Permeability Profiles When Injector-Producer Distance is 300 Ft.
Figure 6.6: Permeability Profile Evolutions For The Base Case Mentioned In Table 5.1 Except
For A Change In Injector-Producer Distance

Critical flow velocity must be met to mobilize the local grain particles, in addition to the physics
of size exclusion in porous media. Assuming the smallest grain particle is 18 microns in diameter,
following the grain size distribution in Figure 5.2, 'safe' injection gradients are calculated for sand
layers with various depths, initial permeabilities, and injector-producer distances. By ‘safe’
injection gradient, it means below which even the smallest particles of 18-micron diameter will not
get dislodged, meaning no possibility to create a thief zone. Open-hole injector with clean sand
in the near-wellbore area is assumed. The bottomhole pressure is assumed to be 14.7 psi at the
producing well. Results are shown in Table 6.13.
For example, for a 3000 ft.- deep unconsolidated sand, if its initial permeability is 2000 mD,
and the injector-producer distance is 100 ft., the 'safe' injection gradient is calculated as 0.73 psi/ft.
If the smallest formation grain is more than 18 microns in diameter, sand migration and erosion
will not happen when water injection is operated below a pressure gradient of 0.73 psi/ft. Although
sand thickness does not affect the sand erosion much in clean formations, note that the critical
injection gradients calculated in Table 6.13 are based on a 10 ft.-thick layer. Similarly, the table
86
include computations for cases where the distances between the injector and producer are 50,
100, and 150 ft, respectively.

Table 6.13: Critical Injection Gradient To Incur Sand Migration And Erosion In Open-Hole
Injection Wells
Initial
Permeability,
mD
800 1500 2000 3000
Injector-
Producer
Distance, ft.
50 100 150 50 100 150 50 100 150 50 100 150
Depth,
ft.
1000 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9
2000 0.9 0.9 0.9 0.84 0.9 0.9 0.63 0.9 0.9 0.42 0.84 0.9
3000 0.9 0.9 0.9 0.49 0.9 0.9 0.36 0.73 0.9 0.2 0.49 0.73
4000 0.6 0.9 0.9 0.26 0.65 0.9 0.2 0.48 0.73 0.13 0.26 0.48
5000 0.35 0.86 0.9 0.19 0.43 0.69 0.14 0.28 0.5 0.09 0.19 0.29

87
Conclusions and Discussions

This chapter includes key conclusions of the base case and sensitivity tests done in Chapter
5 and 6. Moreover, effect of well completion type and formation damage on grain movement and
erosion are explained with examples.

7.1: Conclusions
For unconsolidated formations under waterflooding, particle migration may result in sand
erosion and permeability increase. The corresponding profile damage highly depends on injection
gradient, formation depth, pore size distribution, grain size distribution, and formation temperature.
Based on the model developed in analyzing the grain migration process, we have concluded the
following:
1. The most important factors affecting the computation resolutions are grain size distribution
(GSD) and pore size distribution (PSD). Accurate measurements on GSD and PSD curves with
high data density toward the finer size end are preferred. From the computations, 10+ data points
on GSD for grains smaller than 50 microns, are required to generate satisfying diversity in
realizations. The higher data density on GSD and PSD, the less choppy results are received on
porosity and permeability profiles.
2. Another controlling factor is pressure drawdown per distance. Given a 5000-ft. deep
formation, assuming the smallest grain is 18 microns in diameter, if the producer is distant from
the open-hole injector for 500 ft. or more, particle migration is not happening even at an injection
gradient of 0.85 psi/ft. As a comparison, with the same setting, if the injector-producer distance is
100 ft. or less, the sand erosion can be easily noticed. Therefore, it is not the absolute value of
88
injection pressure or injection gradient that decides whether or not sand migration happens, but
the pressure change per distance that matters.
Note that in radial systems, pressure drawdown per distance in the near-wellbore region is
higher than that for a linear system, which triggers particle migration more easily.
3. Formation depth affects the computation of sand migration from two aspects. On one side,
with the same injection gradient, the deeper formation can take a higher injection pressure, which
can mobilize more bigger-sized particles. On the other side, in the deeper layer where the
formation temperature is higher, not only the water density is lower, both dynamic and kinematic
viscosities of water become lower. Decreases in viscosities make the particle migration harder.
In summary, depth-pressure function and depth-temperature function act against each other, but
former overrides the latter one. The net effect in the deeper formation is to cause more severe
particle migration and sand erosion, if the two formations at different depths have same formation
strength.
4. Initial permeability is an important factor that influences particle migration. Given two cases
where only the initial permeabilities are different, the one with higher initial permeability is more
prone to particle migration and result in sand erosion. It is analogous to the process, as the
formation with a high initial permeability gets damaged more easily and severely by the end of
the drilling operation. The reason is that the formation with a higher permeability allows more
particle invasion which plug up the pores. It is in alignment with Barkman and Davidson (1972)
[12]. For the same principle, in modeling of particle migration induced sand erosion, the inherent
formation particles get moved more easily in a high-permeability porous media.
5. As mentioned above, reservoir temperature affects the viscosities of injected fluid, which
further affect the critical shear velocity to dislodge formation particles, as shown in Table 4.1 and
Figure 5.5. As waterflooding continues, the temperature of injected fluid gradually gets closer to
reservoir temperature. When fluid temperature increases, it becomes harder to mobilize the
89
particle of same size. Therefore, instead of assuming a constant temperature, incorporating the
temperature – viscosity function generates more reasonable results.
6. Assuming open-hole completion and clean formation, layer thickness affects the
computation less significant as other factors. Comparing Table 6.1 to Table 6.5 and Table 6.9, no
obvious trend is seen to conclude whether sand migration prefers to happen in thin or thick layers.
From comparisons among other tables, such as the group composed of Table 6.2, Table 6.6, and
Table 6.10, the group of Table 6.3, Table 6.7, and Table 6.11, as well as the group of Table 6.4,
Table 6.8, and Table 6.12, similar conclusions are reached.
Field observations may suggest a higher possibility of developing a thief zone in thin layers.
It should be noted that, in those cases the layers may also have high initial permeabilities which
are prone to get eroded, as discussed above. Another cause would be the damage developed in
partial sand face around the injection wells due to drilling and completion, which limited the
effective thickness and increased the apparent flow velocity to cause more sand erosion. It will
be discussed in detail in the next section.
Based on our computation, if assuming the entire sand face of the hydrocarbon-bearing
interval is clean, the unconsolidated sand with higher initial permeability subject to a higher
injection rate is more prone to developing a thief zone, which could be thin or thick.

7.2: Discussions
7.2.1: Completion method and formation damage affecting sand migration in
unconsolidated formations
As discussed in Section 1.3 and Section 4.1.2, in designing gravel-packed screens, it has
been realized that the apparent fluid velocity is much higher in a perforated well than in an open-
hole of same wellbore radius, due to the limited flow area caused by perforations [15] [17].
Similarly, if two wellbores are completed in same way, the apparent fluid velocity in a narrower
90
wellbore is higher. In addition to the Figure 4.1 which shows the estimations of apparent fluid
velocity at various completion conditions and wellbore radii, Figure 7.1 shows that the smaller
perforations further increase the apparent flow velocity [15].

Figure 7.1: Fluid Drag Forces Vs. Completion Conditions. ([15] Courtesy SPE.)

The well completion condition, perforation size and density, etc., were not incorporated into
the sensitivity tests presented in Chapter 5 and 6, for the following reasons:
1. Rigorous mathematical relationships in calculating or correlating apparent flow velocity
was not available.
2. Initial permeability of formation is a key factor affecting particle migration, but in either
Figure 4.1 or Figure 7.1, the formation permeability value is missing, so the curves could not be
incorporated into the simulation directly.
However, a general comparison between the cases shown in Figure 4.1 and the validation
example described in Chapter 5 can be done. In the case 4 shown in Figure 4.1, the test was
done at the rate of 100 BPD/ft. in a perforated well, where the perforation density was 4 perfs./ft.
and the perforation diameter was 0.75 in. The apparent flow velocity was tested to be higher than
0.01
0.1
1.0
10
0.01
Legend
(1) 100 BPD/ft 8’ OH
(2) 10 BPD/ft 4 ½ SPF
(3) 20 BPD/ft 4 ¾ SPF
0.1 1.0
ΔP Per Unit Particle Length, psi
N
RE
/(1-Φ)
Cohesive Strength
= 0.35 psi
(Ref: Fig 2-A)
FLUID DRAG FORCE
(Calculated Pressure
Drop Based on 30° API
Oil Flow Through Clean
Formation Sand Pack
Oil Viscosity = 10 cp)
0.003 0.01 0.1
Apparent Velocity, ft/s
(1)
(3)
(2)
91
1 ft./s. Consult the critical shear velocity calculated in Figure 5.5, such a high apparent flow
velocity is able to mobilize grain particles of 0.025 inch in diameter, equivalent to 635 microns.
Note that the diameter of 635 microns corresponds to the largest grains measured from the
unconsolidated sand sample, as shown in Figure 5.2 [16]. However, such big grain particles do
not move at all when the injector was completed as open hole, according to computed sequential
grain size distributions shown in Figure 6.1.
Partial formation damage acts in a similar way as perforations in terms of increasing apparent
flow velocity in the near-wellbore region around the injection well. For a thick or thin layer, if a
portion of the sand face gets damaged so that the effective thickness and flow area are limited,
apparent flow velocity also becomes higher. As a result, the water injection at the same rate can
only flow into the undamaged zones. On profile surveys, it may behave as that the thief zone is
being developed in a thin section of the entire interval.
In summary, high apparent flow velocity could be a result of either perforations or partial
formation damage around injection wells, which may affect sand movement much more
significantly than the key factors mentioned in Chapter 6, including pressure drawdown per lateral
distance, initial formation permeability, injection gradient, formation depth and thickness.
Therefore, investigations on effective flow area and apparent flow velocity should be conducted
before simulating particle migration and sand erosion processes.
Figure 4.1 shows that the apparent flow velocity in a perforated hole is estimated to be 10
2

to 10
3
times higher than in an open hole. Assuming a ten-time increase in apparent flow velocity
and keep all other conditions the same as in Table 5.1, GSD evolutions of the first three timesteps
are generated and shown in Figure 7.2(a) to (d2). In the case, a 10-ft. thick formation with initial
permeability of 2000 mD, located at the depth of 3000 ft., is subject to an injection gradient of 0.75
psi/ft and the apparent flow velocity is increased for 10 times.
Start with the initial condition shown in Figure 7.2(a), particles adjacent to the injection well
migrate substantially just by the end of the first timestep, as shown in Figure 7.2(b), see the
92
several grid blocks at the left end of the slab reservoir. Comparing to the open-hole case shown
in Figure 6.1(b), particles of relatively bigger sizes which do not migrate in the open-hole case
start to travel under the increased apparent flow velocity. As it proceeds to the second and third
timestep, grains from more distant locations join the migration, as shown in Figure 7.2(c), (d1),
and (d2).

Figure 7.2(a): Initial GSD.

Figure 7.2(b): Two Realizations Of GSD By The End Of The First Timestep.

93

Figure 7.2(c): Four Realizations Of GSD By The End Of The Second Timestep.

Figure 7.2(d1): Four Realizations Of GSD By The End Of Third Timestep.

94

Figure 7.2(d2): Another Four Realizations Of GSD By The End Of Third Timestep.
Figure 7.2: Evolution Of Grain Size Distribution On Initial Status, As Well As All Realizations For
The First, Second, And The Third Steps.

Figure 7.3 and Figure 7.4 show the dynamic changes in porosity and permeability profiles for
the same case of Figure 7.2. Sub-figures in Figure 7.3 and Figure 7.4 are generated based on
different injector-producer distance. From Figure 6.5 and Figure 6.6, with the assumption of open-
hole injector, particle migration is trivial when the injector-producer distance is more than 100 ft.
However, in the case where the wellbore is perforated, assuming the apparent flow velocity
increases for 10 times, sand migration resultant porosity and permeability increases are very
obvious, given the injector-producer distances to be 50, 100, 200, and 300 ft., respectively. Figure
7.3(a) to (d) and Figure 7.4(a) to (d) demonstrate the significant effect of perforations on particle
migrations. With less perforation density and diameter, the sand erosion will be more severe.

95

Figure 7.3(a): Porosity Variations When The Injector-Producer Distance Is 50 Ft.

Figure 7.3(b): Porosity Variations When The Injector-Producer Distance Is 100 Ft.

96

Figure 7.3(c): Porosity Variations When The Injector-Producer Distance Is 200 Ft.

Figure 7.3(d): Porosity Variations When The Injector-Producer Distance Is 300 Ft.
Figure 7.3: Porosity Profile Evolutions For The Base Case Mentioned In Table 5.1 Except For A
Change In Injector-Producer Distance And A Change From Open-Hole To Perforated Well.

97

Figure 7.4(a): Permeability Variations When The Injector-Producer Distance Is 50 Ft.

Figure 7.4(b): Permeability Variations When The Injector-Producer Distance Is 100 Ft.

98

Figure 7.4(c): Permeability Variations When The Injector-Producer Distance Is 200 Ft.

Figure 7.4(d): Permeability Variations When The Injector-Producer Distance Is 300 Ft.
Figure 7.4: Permeability Profile Evolutions For The Base Case Mentioned In Table 5.1 Except
For A Change In Injector-Producer Distance And A Change From Open-Hole To Perforated
Well.

Figure 7.5 compares average permeability changes over time among the four cases with
various injector-producer distances. Average permeability increases for all cases are significant
99
as shown in Figure 7.5. Permeability is harmonically averaged using Eq. 4.34. Even when the
injection-producer distance is 300 ft., sand erosion induced average permeability increase for 3.2
times of its initial value is still very severe which can lead to uneven injectivity profiles.
This section highlights the importance of including wellbore completion type in estimating
apparent flow velocity for simulating particle migration and sand erosion. Similar principle applies
to partially damaged sand face around injection wells, which also creates the 'jetting' effect, just
as what perforations do.

Figure 7.5: A Comparison Among The Cases With Various Injector-Producer Distances In A
Perforated Injection Well

7.2.2: Low pore volume injected in developing a thief zone
The dissertation is focused on introducing an algorithm into history matching the
waterflooding performance in unconsolidated formations where sand migration and erosion take
place. Associated changes of petrophysical properties such as porosity and permeability are
drastic so they much be modelled and incorporated into reservoir simulation. The algorithm was
tested for the base case with a 10 ft.-thick layer and cases for 30- or 50 ft.-thick layers. However,
100
It should be realized that the actual thief zones formed under the continuous waterflooding could
be very thin, to the level of inches or less. For that matter, developing a thief zone on field scale
and scale of numerical models may not be on the same order of magnitude in terms of evaluating
pore volume injected (PVI) involved to create a thief zone.
Review the real cases in Figure 2.7 to Figure 2.9, it is noticed that the actual thief zones are
very thin compared to the thicknesses of the entire formations. The values are summarized in
Table 7.1.

Table 7.1: Thickness of Thief Zone Compared to the Formation Thickness
Case
Thief Zone
Thickness, ft.
Formation
Thickness, ft.
Figure 2.7 5 200
Figure 2.8 5 650
Figure 2.9 10 560

For example, considering the thief zone thickness and formation thickness shown in Figure
2.8, if assuming the well placement follows a 5-spot pattern. Calculations indicate that the pore
volume of the entire formation is 2041 times more than the volume of the thief zone, see Table
7.2. Thus, the actual PVI involved on a field scale will be 2041 times smaller to result into a thief
zone. As such, to compare any variables plotted against Pore Volume Injected (PVI) documented
at a field scale with a scale used in lab experiments or numerical simulations can be very
misleading. Attentions on comparison between thief zone thickness and thickness of the entire
layer are required.

Table 7.2: Assumed Properties For Computing The Case Shown In Figure 2.9.
Parameter Value Unit
Injector-Producer Distance 200 ft.
Well Placement 5-spot pattern -
Thief Zone Dimension 200 * 10 * 5 = 1E4 ft.
3

Pore Volume Geometry ¼ of a cylinder -
Pore Volume ¼ * 3.14 * 200
2
* 650 = 2.041E7 ft.
3
101
7.3: Highlights
The modeling of unconsolidated sand erosion is conducted under typical waterflooding
conditions, meaning no intention to create fractures. It features the stochastic process of honoring
possibilities involved in sand migration process. As particles migrate, they may plug up certain
pore throats and accumulate locally which reduces local permeability. On the other hand, particles
may pass through certain pore throats which enhances local permeability. Pore size variations
are realized using the Markov Chain process. Grain size distributions and pore size distributions
are generated any each location and time, which enable the updates of porosity and permeability.
Pressure drawdown per distance, injection gradient, initial permeability, and formation depth
are key factors affecting the severity of sand erosion. Temperature should be considered in order
to deliver more accurate results. Most importantly, as many of the unconsolidated formations are
cased and perforated, also, their high permeability allows considerable invasion during drilling
and completion, evaluations of apparent fluid velocity must be done prior to modeling particle
migration. Perforations and partial formation damage could affect particle migration much more
significantly than the factors mentioned above.
In modeling particle migration and sand erosion, the algorithm developed use easily
accessible data, including pore size distribution, grain size distribution, general reservoir
properties such as initial porosity, permeability, thickness, fluid viscosity, etc. If available, profile
surveys help inversely validate the realizations generated. The algorithm is coupled with a typical
reservoir simulator so that is can be used for history matching the unconsolidated formations
under waterflooding.

102
Suggested Future Work
8.1: Extensions of modeling in multi-phase flow with different
grain wettability, and consideration of formation strength
Since the early discussion [11] of wettability of formation grains, additional experimental work
indicated that wettability is one of the factors affecting invasion distance as a function of residual
oil saturation [56]. However, certain properties to be adopted for refinement of our algorithm need
further investigations. Behaviors of grain particle movements have been classified based on
various grain wettability conditions [29], as shown in Figure 8.1. It is also noticed from lab
measurements that the differences in particle migration rates are substantial when mobile fluid
saturation varies.
Few mathematical relationships representing those conditions are currently available for
modeling of sand migration process. From the previous discussions, grain size is a key factor
affecting particle migration. Also, movement of same-sized grains should be separately
considered when they are water-wet, oil-wet, or of mixed-wettability. As the saturation of wetting
phase changes, grain movement velocity and pathway also change.
From Figure 8.1(a), if movable particles are wet by the flowing phase in a single-phase
system, they continue to flow along with the fluid unless trapped in pore throats. Our algorithm is
developed based on the single-phase condition. In a two-phase environment, if the flowing phase
is at its connate saturation, movable particles will stay in-situ as long as they are wet by flowing
phase, as shown in Figure 8.1(b). It posts an important message that the residual saturation of
flowing phase that wets the grains should be considered in simulation. Figure 8.1(c) demonstrates
that if both wetting phase and non-wetting phase are flowing, fine particles move along with the
wetting phase.
103
If particles are of mixed wettability, as long as flow rate is high enough to trigger their
movement, particles travel along the interface between two fluid phases, see Figure 8.1(d). The
significance of wettability is further shown in Figure 8.1(e), as when mutual solvent flushes
through porous media, movable particles massively migrate with the flowing solvent.
Moreover, in multi-phase flow, saturation of the flowing phase strongly affects formation
strength, see Figure 8.2 [15]. As described in Section 2.2.1, cohesive strength, as a major factor,
differentiates unconsolidated formations from consolidated sands. Cohesive strength also
dominates the evaluation of formation strength for types of unconsolidated sands. When formation
strength varies as a result of changes in saturation, grains with little to no cement will be triggered
to move to very different extents. In our algorithm, one assumption is that there is no cement
bonding the movable fine particles. If to expand the algorithm for multi-phase environment, where
the cement level stands differently due to the changes in fluid saturation, modeling of particle
movement needs to include the functions. To sum up, saturation of flowing phase not only
influences the movability and flowing pathway of particles, but also affects formation strength.
Therefore, in the case of history matching waterflooding performance in an unconsolidated
reservoir, quantitative investigations of grain wettability and changes in fluid saturations are
needed, to refine the algorithm.

Figure 8.1(a): When A Single-Fluid Phase Presents, Fines (Wet By The Single-Phase)
Move With The Flowing Fluid Unless Bridged At Pore Throats.

Mobile Fines
Fluid Flow Direction
Fines Bridged at
Pore Restriction
(a)
104

Figure 8.1(b): Water-Wet Fines Are Immobile When The Water Phase Is Immobile.

Figure 8.1(c): Water-Wet Fines Not Bridged At Pore Restrictions Are Mobile When Both
Water And Oil And Flowing.

Figure 8.1(d): Fines Of Mixed Wettability Constrained To Move Only Along The Oil/Water
Interface.

Immobile
Water-Wet
Fines
Connate
water
(Immobile)
Oil
(b)
Mobile Water
Mobile Water
(c)
Oil
Mixed
Wettability
Fines
Connate
water
Oil
(d)
105

Figure 8.1(e): Mutual Solvents Release Fines Held By Wetting And Interfacial Forces,
Causing Them To Migrate At High Concentrations.
Figure 8.1: Observations On Wettability And Multi-Phase Affecting Particle Movement.
([47]. Courtesy SPE.)

Figure 8.2(a): Cohesive Test Run On-Air/Water System.

Water
Mutual Solvent
Oil
Oil
Oil
Oil
(e)
0
0
Air/Water Capillary Pressure, psi
2
4
6
Saturation of Wetting Phase, %
20 40 60 80
100
Pore Entry
Pressure
(a)
106

Figure 8.2(b): Sand Size Distribution Of The Sample Tested.

Figure 8.2(c): Cohesive Strength Generated With Measurements Shown In
Figure 8.2(a).
Figure 8.2: Formation Strength Affected By Saturation In Multi-Phase Condition. ([15]. Courtesy
SPE.)

8.2: Other aspects of work to be clarified related to the sand
particle movement process
Generally, sand movement processes have been studied for the sake of understanding
formation damage, sand erosion, and sand production, etc. However, many observations from
0
20
40
60
80
100
Cumulative Weight Percentage, %
Screen Opening, U.S. Mesh
60 80 100 140 200 270 325
(b)
0 20 40 60 80 100
Water Saturation, %
Cohesive Strength, psi
0.2
0
0.4
0.6
0.8
1
Oil=water System
(Calculated from Air-Brine
Capillary Pressure Curve)
(c)
107
labs and fields need to be further clarified, which not only help sand control and water control
procedures, but also benefit the modeling of particle migration.
As related to formation damage, several parts of the problem are hard to quantify, such as
the initial skin of formation created by drilling and completion. It should be realized that the initial
skin may not exist or stands differently for each well [29]. More importantly, damaging particles
should be separately treated from suspending fine particles in porous media. These particles, that
cause formation damage, have not gone through enough diagenesis or cementation when
considering fine movement. As such, real situation complicates the modeling process.
Some field practices are opposite from the previous experimental studies. For instance, water
filtering in the field often uses sand filters with backwashes when seeing a significant drop in
filtration speed. Also, it is known that many oil field operators periodically backflush producing
wells in recovering the productivity of those wells. The same has been applied on injectors to
recover injectivities as well. However, one study [13] has tested with cores and concluded that,
once the samples have been plugged, backwashing does not improve the permeability to a
considerable extent. While reasons for the differentiations remain unclear, it reminds researchers
that calibration of lab experiments with field observations should be done before simulating
particle migration process.
From the literature review, it is also noted that pore pressure, as a factor affecting water
injectivity, needs further investigations. The studies [45] [46] [47] presented analytical or semi-
analytical solutions to model pressure transient behaviors in formations embedded with pancake-
shaped high conductivity zones. However, the relationship between pore pressure profile and
elastic changes in pore space and geometry have not been modeled, while pore space and
geometry are important for particle migration due to the size exclusion nature of porous media.
Also, lab tests on TDS of injected water and compatibility between injected water and formation
water are needed in modeling the capability of particle movements under such various
environments.
108
In the work done by Davies and Davies (2001) [32], three samples of different rock types are
tested, where each has its own packing feature, grains angularity, and cement level. In reality,
more tests are needed for real sands with embedded clays of different types, and their impact on
water injectivity. Except for the analytical solution to deep filtration process [58], most of the factors
considered in literature have barely been upscaled, before being incorporated into numerical
simulations.
Grain breakage has been discussed but not extensively correlated in laboratories, which lead
to changes in grain size and geometry. Unsteady state water injection with suspending particles
and oil residuals has been examined in lab tests but not mathematically correlated or numerically
simulated. However, it is not negligible, as in unsteady state, particle mobilizations are much more
severe due to the higher pressure perturbations.
In order to simulate particle movement processes more comprehensively and accurately for
history matching the performance of unconsolidated formations, the above aspects are of interest
and need to be considered.
109
Nomenclature
𝑎 a fraction of bulk volume taken by non-movable grains, dimensionless
𝐴 coefficient describing tortuosity of porous media, dimensionless
𝑏 a fraction of bulk volume taken by movable grains, dimensionless
𝐶 coefficient describing porous media sand sorting, dimensionless
𝐶
L
total compressibility, 1/𝑝𝑠𝑖
𝜓 grain sphericity, dimensionless
𝐷 formation depth, 𝑓𝑡.
𝐷𝑖𝑠𝑡 traveling distance of grains, 𝑓𝑡.
𝐷
17
,𝐷
87
the diameter of grain at 40 and 90 weight percentiles, 𝑖𝑛.
Δt simulation timestep size, 𝑑𝑎𝑦
Δx,Δy,Δz
𝑓
simulation grid size in x, y, and z-direction, 𝑓𝑡.
filtration factor, dimensionless
𝛷,𝛷
9Q9
porosity and initial porosity, dimensionless
𝑔 gravitational acceleration, 𝑓𝑡./𝑠
G

𝐺𝑃 grain population, count
𝐺𝑃𝐷 grain population density, dimensionless
𝐺𝑆 grain size, 𝑖𝑛., or 𝜇𝑚
𝐺𝑉 individual grain volume, 𝑖𝑛
2
or 𝜇𝑚
2

𝐺𝑊 incremental grain weight percentage, dimensionless
ℎ formation thickness, 𝑓𝑡.
𝐼𝑛𝑗
Izp3
injection gradient, 𝑝𝑠𝑖/𝑓𝑡.
𝑘,𝑘
9Q9
permeability and initial permeability, 𝑚𝑑
110
𝐾
7
permeability increases, dimensionless
𝐾
pE-
7
average permeability increase, dimensionless
𝑙 number of loops, count
𝑙
I
mean grain size, 𝑖𝑛. or 𝜇𝑚
𝐿
f
,𝐿
x
,𝐿
y
simulation reservoir dimension in x, y and z directions, 𝑓𝑡.
𝑚 porous media bulk mass, 𝑙𝑏
2

𝜇 fluid dynamic viscosity, 𝑐𝑃
𝑣 kinematic viscosity, 𝑓𝑡.
G
/𝑠
𝑛
I
number of grain sizes, count
𝑛
[
number of pore sizes, count
𝑛
L
number of timesteps in simulation, count
𝑛
f
number of grids along x-direction in simulation, count
𝑃,𝑃
9Q9

𝑃𝐿𝐹
𝑃𝐿𝑅
reservoir pressure and initial reservoir pressure, 𝑝𝑠𝑖
the fraction of particle loss per timestep, dimensionless
critical residual particle fraction, dimensionless
𝑃𝑟𝑜𝑏 probability of one-grain particle passing through the larger pores
𝑃𝑆 pore size, 𝜇𝑚
𝑃
LP
,𝑃
LG
pore size distribution at current and future states, dimensionless
𝑅𝑒
∗
grain Reynolds number, dimensionless
𝜌 fluid density, 𝑙𝑏/𝑓𝑡.
2

𝜌
I
grain density, 𝑙𝑏/𝑓𝑡.
2

𝑡 time, 𝑑𝑎𝑦
𝑡
LoLpJ
simulation total time duration, 𝑑𝑎𝑦
𝑢
∗
critical shear velocity, 𝑓𝑡./𝑠
111
𝑉
@
fluid velocity, 𝑓𝑡./𝑠
𝑉
I
grain velocity, 𝑓𝑡./𝑠
𝑤
/
settling velocity, 𝑓𝑡./𝑠
𝑥 location in the x-direction, 𝑓𝑡.

112
References
[1] B. Bailey et al., “Water Control,” Oilfield Review, p. 22, 2000.
[2] F. L. Hill, S. Monroe, and R. Mohanan, “Water Management-An Increasing Trend in the Oil
and Gas Industry,” in SPE/EAGE European Unconventional Resources Conference and
Exhibition, Vienna, Austria, 2012, doi: 10.2118/154720-MS.
[3] A. Zeinijahromi and P. Bedrikovetski, “Controlling Excessive Water Production Using
Induced Formation Damage,” in SPE European Formation Damage Conference and
Exhibition, Budapest, Hungary, 2015, doi: 10.2118/174229-MS.
[4] K. Lee and J. Neff, Eds., Produced Water. New York, NY: Springer New York, 2011.
[5] T. Bakke, J. Klungsøyr, and S. Sanni, “Environmental impacts of produced water and drilling
waste discharges from the Norwegian offshore petroleum industry,” Marine Environmental
Research, vol. 92, pp. 154–169, Dec. 2013, doi: 10.1016/j.marenvres.2013.09.012.
[6] M. T. Stephenson, “A survey of produced water studies,” Technological/environmental
issues and solutions., pp. 1–11, 1992.
[7] P. McCormack, P. Jones, M. J. Hetheridge, and S. J. Rowland, “Analysis of oilfield produced
water and production chemicals by electrospray ionization multistage mass spectrometry
[ESI-MSn],” Water Res., vol. 35, no. 15, pp. 3567–3578, 2001.
[8] M. T. Stephenson, “Components of produced water: a compilation of results from several
industry studies,” presented at the Proceedings international conference on health, safety
and environment, Hague, Netherlands, Nov. 1991, pp. 25–38.
[9] K. S. Chan, “Water Control Diagnostic Plots,” p. 9.
[10] I. Kozyrev, P. Zyatikov, and V. Deeva, “Solving Excess Water Production Problems In
Productive Formation,” presented at the MATEC Web of Conferences, 2016, vol. 72, doi:
10.1051/matecconf/20167201049.
113
[11] H. W. Milton, T. R. Laverty, G. R. Bystrom, and E. T. Guerrero, “Effect Of Grain Size, Grain
Shape And Clay Particles On Oil Recovery By Water Flooding.,” Jan. 1961.
[12] J. H. Barkman and D. H. Davidson, “Measuring Water Quality and Predicting Well
Impairment,” Journal of Petroleum Technology, vol. 24, no. 07, pp. 865–873, Jul. 1972, doi:
10.2118/3543-PA.
[13] D. H. Davidson, “Invasion And Impairment Of Formations By Particulates,” presented at the
54th ATCE of the Society of Petroleum Engineers of AIME, Jan. 1979, doi: 10.2118/8210-
MS.
[14] T. Hofsaess and W. Kleinitz, “30 Years of Predicting Injectivity after Barkman & Davidson:
Where are we today?,” presented at the SPE European Formation Damage Conference,
Jan. 2003, p. 12, doi: 10.2118/82231-MS.
[15] J. L. Durrett, W. T. Golbin, J. W. Murray, and R. E. Tighe, “Seeking a Solution to Sand
Control,” Journal of Petroleum Technology, vol. 29, no. 12, pp. 1664–1672, Dec. 1977, doi:
10.2118/6210-PA.
[16] D. H. Schwartz, “Successful Sand Control Design for High Rate Oil and Water Wells,”
Journal of Petroleum Technology, vol. 21, no. 09, pp. 1193–1198, Sep. 1969, doi:
10.2118/2330-PA.
[17] L. D. Lott, R. Viney, and S. J. Abbott, “Survey of Sand Control (Southern District Study
Commite on Sand Control),” Jan. 1954.
[18] A. R. Jennings, “Laboratory Studies of Fines Movement in Gravel Packs,” SPE Drilling &
Completion, vol. 12, no. 04, pp. 275–281, Dec. 1997, doi: 10.2118/36420-PA.
[19] M. S. Bruno, A. Dorfmann, K. Lao, and C. Honeger, “Coupled particle and fluid flow modeling
of fracture and slurry injection in weakly consolidated granular media,” Jan. 2001.
[20] J. Hwang and M. M. Sharma, “Filtration in Frac Packs and its Impact on Injector
Performance,” SPE Production & Operations, vol. 29, no. 03, pp. 204–215, Aug. 2014, doi:
10.2118/168135-PA.
114
[21] D. A. Shumbera, D. W. Ritter, R. C. Ellis, R. C. Jannise, and G. Lochte, “Improved Water
Injector Performance in a Gulf of Mexico Deepwater Development Using an Openhole Frac
Pack Completion and Downhole Filter System: Case History,” p. 11.
[22] N. Stein and D. W. Hilchie, “Estimating the Maximum Production Rate Possible from Friable
Sandstones Without Using Sand Control,” Journal of Petroleum Technology, vol. 24, no. 09,
pp. 1157–1160, Sep. 1972, doi: 10.2118/3499-PA.
[23] J. A. Robertson, J. A. Blesener, and S. S. Hoo, “Subzone Redevelopment of the Long Beach
Unit, Wilmington Oil Field: A Case Study,” Journal of Petroleum Technology, vol. 39, no. 10,
pp. 1229–1236, Oct. 1987, doi: 10.2118/15101-PA.
[24] S. J. Jeu et al., “Case Histories of Unique Multi-Zone Intelligent Deepwater Sand-Control
Completions - How Successes Were Achieved and the Associated Lessons Learned.
Offshore Technology Conference,” Jan. 2008, doi: doi:10.4043/19622-MS.
[25] O. Karazincir, W. Williams, and P. Rijken, “Prediction of Fines Migration through Core
Testing,” in SPE Annual Technical Conference and Exhibition, San Antonio, Texas, USA,
2017, doi: 10.2118/187157-MS.
[26] B. H. Adams, “Stress-Sensitive Permeability in a High-Permeability Sandstone Reservoir-
The Kuparuk Field,” in SPE California Regional Meeting, Ventura, California, 1983, doi:
10.2118/11718-MS.
[27] Z. Yang, I. Ershaghi, J. J. Mondragon, and S. Hara, “Method for Handling the Complexities
Associated with History Matching the Performance of a Highly Stress-Sensitive Formation,”
p. 15.
[28] M. N. Mayuga, “Geology and Development of California Giant Wilmington Oil Field,” AAPG
Mern., vol. 14, pp. 158–184, 1970.
[29] T. W. Muecke, “(done) Formation Fines and Factors Controlling Their Movement in Porous
Media,” Journal of Petroleum Technology, vol. 31, no. 02, pp. 144–150, Feb. 1979, doi:
10.2118/7007-PA.
115
[30] F. L. Monus, F. W. Broussard, and D. Schlumberger, “Fracturing Unconsolidated Sand
Formations Offshore Gulf of Mexico,” p. 15.
[31] D. Norman, “Overview: Frac Pack,” Journal of Petroleum Technology, vol. 56, no. 09, p. 70,
Sep. 2004, doi: doi:10.2118/0904-0070-JPT.
[32] J. P. Davies and D. K. Davies, “Stress-Dependent Permeability: Characterization and
Modeling,” SPE Journal, p. 12, 2001.
[33] I. N. Sneddon, “The distribution of stress in the neighbourhood of a crack in an elastic solid,”
Proceedings of Royal Society, vol. 187, no. 1009, Oct. 1946, doi:
https://doi.org/10.1098/rspa.1946.0077.
[34] M. A. Biot, L. Masse, and W. L. Medlin, “A Two-Dimensional Theory of Fracture Propagation,”
SPE Production Engineering, vol. 1, no. 01, pp. 17–30, Jan. 1986, doi: 10.2118/11067-PA.
[35] Y. Abousleiman, Y. G. Yuan, and J. C. Roegiers, “Propagation/recession of a penny-shaped
fracture in a permeable diffusive formation,” Jan. 1994, doi: 10.2118/28086-MS.
[36] X. Zhang, E. Detournay, and R. Jeffrey, “Propagation of a penny-shaped hydraulic fracture
parallel to the free-surface of an elastic half-space,” International Journal of Fracture, p.
115:125, 2002, doi: org/10.1023/A:1016345906315.
[37] D. R. Holbert and R. O. Perry, “Selective Plugging of Water-injection Wells with
Polymerizable Materials,” 1958, doi: API-58-034.
[38] J. O. Robertson and F. H. Oefelein, “Plugging Thief Zones in Water Injection Wells,” Journal
of Petroleum Technology, vol. 19, no. 08, pp. 999–1004, Aug. 1967, doi: 10.2118/1524-PA.
[39] H. N. Hall, “How to Analyze Waterflood Injection Well Performance,” World Oil, Oct. 1963.
[40] Q. Qi, S. Pepin, A. Al-Jazzaf, and I. Ershaghi, “Errors Associated with Waterflood Monitoring
Using the Hall Plot for Stacked Reservoirs in the Absence of Profile Surveys,” Apr. 2017,
doi: 10.2118/185694-MS.
116
[41] S. D. Yoelin, “Water Injection in the Lower Jones Sands, Huntington Beach Offshore Field,”
Journal of Petroleum Technology, vol. 20, no. 02, pp. 127–134, Feb. 1968, doi:
10.2118/1989-PA.
[42] S. D. Yoelin, C. D. Howald, and D. R. Holbert, “Production Logging as Used To Solve Water
Injection Problems in the Huntington Beach Offshore Field,” Journal of Petroleum
Technology, vol. 22, no. 09, pp. 1083–1088, Sep. 1970, doi: 10.2118/2684-PA.
[43] J. C. Jr. Calhoun, Fundamentals of Reservoir Engineering. University of Oklahoma Press,
Norman, OK, 1953.
[44] W. R. Watkins, “How To Diagnose a Thief Zone,” Society of Petroleum Engineers, Jul. 1973,
doi: 10.2118/4416-PA.
[45] A. C. Gringarten and H. J. Ramey, “Unsteady-State Pressure Distributions Created by a Well
With a Single Horizontal Fracture, Partial Penetration, or Restricted Entry,” Society of
Petroleum Engineers Journal, vol. 14, no. 04, pp. 413–426, Aug. 1974, doi: 10.2118/3819-
PA.
[46] P. Valko and M. J. Economides, “Transient Behavior of Finite Conductivity Horizontal
Fractures,” SPE Journal, vol. 2, no. 02, pp. 213–222, Jun. 1997, doi: 10.2118/38436-PA.
[47] Q. Qi and I. Ershaghi, “Fall-off Diagnostic of Eroded Zones during Waterflooding of
Unconsolidated Formations,” presented at the SPE Western Regional Meeting, Apr. 2019,
doi: 10.2118/195376-MS.
[48] A. Satman, “An Analytical Study of Transient Flow in Stratified Systems With Fluid Banks,”
presented at the 56th Annual Fall Conference and Exhibition of the Society of Petroleum
Engineers of AIME, 1981, p. 12.
[49] T. W. Muecke, “Formation Fines and Factors Controlling Their Movement in Porous Media,”
Journal of Petroleum Technology, vol. 31, no. 02, pp. 144–150, Feb. 1979, doi:
10.2118/7007-PA.
117
[50] F. Civan, Reservoir Formation Damage - Foundamentals, Modeling, Assessment and
Mitigation. 2017.
[51] G. O. Suman, R. C. Ellis, and R. E. Snyder, Sand Control Handbook: Prevent Production
Losses and Avoid Well Damage with These Latest Field-proven Technique. Gulf Publishing
Company, 1983.
[52] A. Rostami and D. Taleghani, “Modeling Particle Mobilization in Unconsolidated Formations
Due to Fluid Injection,” p. 10, 2014.
[53] I. Ershaghi, R. Hashemi, S. C. Caothien, and D. Abdassah, “Injectivity Losses Under Particle
Cake Buildup and Particle Invasion,” Jan. 1986, doi: doi:10.2118/15073-MS.
[54] A. C. Payatakes, C. Tien, and R. M. Turian, “A new model for granular porous media: Part I.
Model formulation.,” AIChE Journal, vol. 19, no. (1), pp. 58–76, 1973, doi: DOI:
10.1002/aic.690190110.
[55] A. C. Payatakes, R. Rajagopalan, and C. Tien, “Application of Porous Medium Models to the
Study of Deep Bed Filtration,” The Canadian Journal of Chemical Engineering, vol. 52, no.
6, pp. 722–731, 1974.
[56] M. A. J. Ali, P. K. Currie, and M. J. Salman, “The Effect of Residual Oil on Deep-Bed Filtration
of Particles in Injection Water,” p. 8, Feb. 2009, doi: doi:10.2118/107619-PA.
[57] R. G. Guedes, F. A. H. Al-Abduwani, P. Bedrikovetsky, and P. K. Currie, “Deep-Bed Filtration
Under Multiple Particle-Capture Mechanisms,” SPE Journal, vol. 14, no. 03, pp. 477–487,
Sep. 2009, doi: 10.2118/98623-PA.
[58] A. S. L. Vaz, P. G. Bedrikovetsky, C. J. A. Furtado, and A. L. Serra de Souza, “Upscaling of
Deep Bed Filtration From Core to Near-Wellbore Region,” in SPE International Symposium
and Exhibition on Formation Damage Control, Lafayette, Louisiana, USA, 2010, doi:
10.2118/128372-MS.
[59] J. F. B. Barake, “Modelling of Channelization in Unconsolidated Formations Due to Liquid
Injection,” 2015.
118
[60] M. M. Sharma, Y. C. Yortsos, and L. L. Handy, “Deposition and Release of Clays in
Sandstones,” presented at the SPE Symposium on Oilfield and Geothermal Chemistry, Apr.
1985, doi: 10.2118/13562-MS.
[61] H. A. Ohen and F. Civan, “Predicting Skin Effects Due to Formation Damage by Fines
Migration,” 1991, p. 12.
[62] P. Bedrikovetsky et al., “Porous Media Deposition Damage from Injection of Water with
Particles,” presented at the ECMOR VIII - 8th European Conference on the Mathematics of
Oil Recovery, Freiberg, Germany, Sep. 2002, doi: 10.3997/2214-4609.201405931.
[63] P. Bedrikovetsky and P. Tran, “Damage Characterization of Deep Bed Filtration from
Pressure Measurements,” p. 11.
[64] P. Bedrikovetsky et al., “Inverse Problems for Treatment of Laboratory Data on Injectivity
Impairment,” p. 7.
[65] A. Afrough, M. S. Zamiri, L. Romero-Zerón, and B. J. Balcom, “Magnetic-Resonance
Imaging of Fines Migration in Berea Sandstone,” SPE Journal, vol. 22, no. 05, pp. 1385–
1392, Oct. 2017, doi: 10.2118/186089-PA.
[66] S. Borazjani, L. Chequer, T. Russell, and P. Bedrikovetsky, “Injectivity Decline During
Waterflooding and PWRI due to Fines Migration,” in SPE International Conference and
Exhibition on Formation Damage Control, Lafayette, Louisiana, USA, 2018, doi:
10.2118/189521-MS.
[67] K. E. Wennberg and A. Hansen, “Modelling Fines Mobilization, Migration and Clogging,” p.
9.
[68] G. J. Plisga, “Chapter 31: Temperature in Wells,” Society of Petroleum Engineers, 1987, p.
7.
[69] L. L. Lacy and H. G. Hudson, “New Step-Rate Test Analysis for Fracture Evaluation,” p. 12.
[70] A. M. Garon, C. Y. Lin, and V. A. Dunayevsky, “Simulation of Thermally Induced Waterflood
Fracturing in Prudhoe Bay,” p. 12.
119
[71] D. B. Williams and D. W. Sherrard, “Impact of Inducing Fractures at Prudhoe Bay,” JPT, Oct.
1989, doi: 10.2118/16358-PA.
[72] D. Paphitis, “Sediment movement under unidirectional flows: an assessment of empirical
threshold curves,” Coastal Engineering, vol. 43, no. 3–4, pp. 227–245, Aug. 2001, doi:
10.1016/S0378-3839(01)00015-1.
[73] L. C. van Rijn, Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas. Aqua
publications, 1993.
[74] P. C. Carmen, “Fluid flow through granular beds,” Transactions, Institution of Chemical
Engineers, vol. 15, pp. 150–166, 1937.

Abstract (if available)

Abstract Recovery of hydrocarbons from unconsolidated sands under waterflooding frequently results in excess water and sand production. Such operational practices may cause reservoir subsidence and damage to production facilities with associated environmental and safety hazards. During the lifespan of a reservoir, these two issues can severely and negatively affect operations by overloading subsurface and surface facilities and pipelines, thereby decreasing hydrocarbon productivity and substantially increasing operational cost. These often happen during water injection into unconsolidated formations, where unbonded grains in some layers can migrate along with fluid flow under injection pressure. Left behind are unswept zones and layers, as more in-situ grains are dislodged and carried away from the near-wellbore area, causing erosion and development of high permeability zones. Associated changes in reservoir properties, such as porosity and permeability, are not negligible and must be modeled. Modeling of sand migration, as a predictive tool for identification of thief zones in a timely manner, should be incorporated as a critical part of reservoir simulation studies. ❧ In this study, a new algorithm modeling sand particle migration is constructed assuming single-phase flow in unconsolidated sands composed of unbonded, spherical grains following typical grain size distributions. The formation pore size distribution assumed is also in the range reported for these types of formations. The algorithm developed is capable of providing statuses of grain size distribution at any stage of waterflooding. As an indication of a washed-out zone or a thief zone being developed, the computed sequential grain size distributions model the migration of various-sized particles, which are further used in computing the dynamics of reservoir permeability. ❧ Given adequate pore sizes allowing grain migration, our method calculates the movement of grains of various sizes when local apparent fluid velocity exceeds critical shear velocity of grains. In this way, the physics of grain mobilization is better honored as one can anticipate that grains are subject to gravitational forces and viscous forces in a porous system. ❧ As possibilities exist for both pore throat bridging and pore throat enlargement due to the local particle movement, variations in pore size distribution are modeled stochastically by incorporating the Markov Chain probability transition process. During the computation, it is noted that detailed measurements of grain size distribution, especially for the smaller-sized end, are needed. Monitoring of formation temperature, which affects fluid viscosity and particle mobilization, is also important. The algorithm developed requires minimal measurements of reservoir characteristics, where only general reservoir properties, pore size distribution, and grain size distribution are required. It also allows validation and re-selections of particle loss factor and transition matrices through comparison to periodically measured profile surveys. ❧ The study includes a computational example using a typical grain size distribution for unconsolidated sands presented in the literature with the assumption of other reservoir properties. The results show that at each location and time, several realizations of grain size distribution are mapped. Permeability is calculated based on counting the particle population density for each realization of grain size distribution. This way, with a profile survey, realizations that match the real permeability increase shown on the survey can be found, as the possible reflections of reservoir dynamics. In this study, the computed permeability increase is tabulated for scenarios considering formation depth, initial permeability, injection pressure gradient, and layer thickness. Moreover, factors including pressure drawdown per unit length, well completion type, and partial formation damage are emphasized. ❧ With the improvement in accuracy of representing formation strength under varying water saturation and grain wettability conditions, the sand migration model presented here can guide the prediction and timely remediation measures from an early stage of reservoir development for recovery under water injection. More importantly, it must be included into reservoir simulations to history match the waterflooding data in unconsolidated sands.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Modeling and simulation of complex recovery processes

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

A method for characterizing reservoir heterogeneity using continuous measurement of tracer data

PDF

Subsurface model calibration for complex facies models

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

The study of CO₂ mass transfer in brine and in brine-saturated Mt. Simon sandstone and the CO₂/brine induced evolution of its transport and mechanical properties

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Reactivation of multiple faults in oilfields with injection and production

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Flood front tracking and continuous recording of time lag in immiscible displacement

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Application of data-driven modeling in basin-wide analysis of unconventional resources, including domain expertise

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

Asset Metadata

Creator Qi, Qianru (author)

Core Title Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 12/12/2020

Defense Date 04/16/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag fine particle migration,OAI-PMH Harvest,sand erosion,thief zone,unconsolidated sands,waterflooding

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Ershaghi, Iraj (committee chair), Chang, Jincai (committee member), Paul, Donald (committee member), Roberts, Richard (committee member), Shing, Katherine (committee member)

Creator Email qianru.qi7@gmail.com,qianruqi@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-412824

Unique identifier UC11667439

Identifier etd-QiQianru-9101.pdf (filename),usctheses-c89-412824 (legacy record id)

Legacy Identifier etd-QiQianru-9101.pdf

Dmrecord 412824

Document Type Dissertation

Rights Qi, Qianru

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA