Controlling Viscous Fingering
Tim H. Beeson-Jones
Supervisor: Prof. Andrew W. Woods
BP Institute
University of Cambridge
This dissertation is submitted for the degree of
Doctor of Philosophy
Clare College April 2018

To Mum, Dad and Jess.

Declaration
I hereby declare that except where specific reference is made to the work of others,
the contents of this dissertation are original and have not been submitted in whole or
in part for consideration for any other degree or qualification in this, or any other uni-
versity. This dissertation is my own work and contains nothing which is the outcome
of work done in collaboration with others, except as specified in the text and Acknowl-
edgements. This dissertation contains fewer than 40,000 words including appendices,
bibliography, footnotes, tables and equations and has fewer than 50 figures.
Tim Beeson-Jones
July 2017
“Climate change presents a unique challenge for
economics: it is the greatest and widest-ranging
market failure ever seen.”
Sir Howard Stern (2006), Stern Review: The Economics of Climate Change

Acknowledgements
First and foremost, I thank my supervisor Andy Woods. I am very grateful for the
patience, support and guidance he has given me during the course of my studies at the
BPI. I hope I have learned something of his approaches to posing problems, solving
problems and communicating the results, approaches which are both unique and world
class. I also hope his humour and enthusiasm for science will remain with me during
my career.
The people at the BPI have created a wonderful community for which I shall soon
feel nostalgic. I am very appreciative of Lotty Gladstone’s efforts towards managing
and keeping safe the brilliant flow labs and also her guidance with experiments. An-
drew Pluck’s ability to make (fix) the things I’ve needed (broken) has been crucial.
Chris Richardson and Patrick Welche have been incredibly helpful in their provision
of a reliable computer system and computing assistance. Finally, I am very grateful
to Catherine Pearson who keeps the BPI a well-oiled machine!
As for my colleagues, I thank Alex Evans, Peter Dudfield, Marcus Horsley, Ed
Hinton and Neeraja Bhamidipati for their crucial mathematical assistance, Finn Box
for sharing his image processing expertise, Caren Lacy and Sefano Rocco for their help
in the lab and Martin Lippert, Jeffrey Poon and Richard Alloway for their camaraderie.
As for the wider university, I am extremely grateful to all those at Clare College
that have made it my home close to my home for eight short years. It is a beautiful
college. In particular, I would like to thank Phil and all the porters for looking out
for us students and Lee and Bonna for keeping us well fed. Also, I am very grateful to
viii
Steve & Sarah Morgan for their tireless efforts at the UPC and all the guys at CUSSC
for showing me how to ski.
There has not been a dull moment amongst my marvellous friends. I am profoundly
grateful to my long-suffering housemates Andy Wedlake, Will Wynell-Mayow, AP
Stafford, Alex Fokas, Ben Stretch and Fred Dudbridge but most of all to the longest
sufferer Lois Salem, who has endured 5 years.
Finally, I thank my family for their immense love and support.
Abstract
Viscous fingering occurs when one fluid displaces another fluid of a greater viscosity
in a porous medium or a Hele-Shaw cell. Linear stability analysis is used to predict
methods of suppressing instability. Then, experiments in which nonlinear growth
dominates pattern formation are analysed to explore the nonlinear impact of strategies
of suppressing finger growth.
Often, chemical treatment fluid is injected into oil reservoirs in order to prevent
sand production. This treatment fluid is usually followed by water injection to clean
up the well. We explore the potential for viscous instability of the interface between
the treatment fluid and the water, and also the treatment fluid and the oil, as a
function of the volume of treatment fluid and the injection rate and viscosity ratios of
the different fluids. For a given volume of treatment fluid and a given injection rate,
we find the optimal viscosity of the treatment fluid to minimise the viscous instability.
In the case of axisymmetric injection, the stabilisation associated with the azimuthal
stretching of modes leads to a further constraint on the optimisation of the viscosity.
In the case of oil production, polymers may be added to the displacing water
in order to reduce adverse viscosity gradients. We also explore the case in which
these polymers have a time-dependent viscosity, for example through the slow release
from encapsulant. We calculate the injection flow rate profile that minimises the
final amplitude of instability in both rectilinear and axisymmetric geometries. In a
development of the model, we repeat the calculation for a shear-thinning rheology.
Finally, experiments are analysed in which the nonlinear growth of viscous fingers
develops to test the influence of different injection profiles on the development of insta-
xbility. Diffusion Limited Aggregation (DLA) simulations are performed for comparison.
In all cases, the evolving pattern has a saturation distribution, with an inner zone in
which the fingers are static and an outer zone in which the fingers advance and grow.
In the very centre of the viscous fingering patterns, there is a small fully-saturated
region. In the experiments, the mass distribution in the inner zone varies with radius
as a power law which relates to the fractal dimension for the analogue DLA simula-
tions. In the outer region the saturation decreases linearly with radius. The radius
of the inner frozen zone is approximately 2/3 of the outer radius in the cases of DLA
and – after a period of evolution – the viscous fingering experiments. This allows the
radial extents of the inner and outer zones to be predicted. The ratio of each radius
to the extent of the fully-saturated region is independent of the injection profile and
corresponds to values for DLA.
Publications
This thesis contains two published studies as follows:
• Chapter 2: Beeson-Jones, T. H. and Woods, A. W. (2015). On the selection
of viscosity to suppress the Saffman Taylor instability in a radially spreading
annulus. Journal of Fluid Mechanics, 782:127–143.
• Chapter 3: Beeson-Jones, T. H. and Woods, A. W. (2017). Control of viscous
instability by variation of injection rate in a fluid with time-dependent rheology.
Journal of Fluid Mechanics, 829:214–235.
“Knowledge is Porridge.”
Stewart Pearson, The Thick of It

Table of contents
List of figures xvii
List of tables xxvii
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Linear Stability Analysis: Predicting Interfacial Control . . . . . . . . . 4
1.3 Nonlinear Growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Annular Saffman-Taylor Suppression 13
2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Single-Interface Stability: Treatment Injection Phase . . . . . . . . . . 17
2.3.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.2 Absolute Stability of mode n . . . . . . . . . . . . . . . . . . . . 18
2.3.3 Dynamic Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.4 Single-Interface Control . . . . . . . . . . . . . . . . . . . . . . 20
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase . . . . . 22
2.4.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.2 Absolute Stability of Mode n . . . . . . . . . . . . . . . . . . . 25
2.4.3 Dynamic Stability . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5 Control Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
xiv Table of contents
3 Viscous instability with time-dependent rheology 39
3.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.3 Newtonian Flow in a Rectilinear Geometry . . . . . . . . . . . . . . . . 44
3.3.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3.2 Effect of a Gradual Increase in Viscosity . . . . . . . . . . . . . 46
3.4 Injection from a point source . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4.1 The Constant Viscosity Regime . . . . . . . . . . . . . . . . . . 52
3.4.2 Effect of a Gradual Increase in Viscosity . . . . . . . . . . . . . 60
3.5 Injection of a Shear-Thinning Fluid . . . . . . . . . . . . . . . . . . . . 63
3.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4 Dynamics of Fractal Growth 71
4.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2.1 Viscous Fingering and the Confinement of Growth . . . . . . . . 72
4.2.2 Fractal Objects, Diffusion Limited Aggregation and the Deter-
mination of the Fractal Dimension . . . . . . . . . . . . . . . . 73
4.2.3 The Fractal Dimension of Viscous Fingering Patterns . . . . . . 75
4.2.4 Chapter Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3.1 Apparatus, Experimental Conditions and Protocol . . . . . . . . 79
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.4.1 Measured Flow Rates . . . . . . . . . . . . . . . . . . . . . . . . 80
4.4.2 VF Pattern Formation . . . . . . . . . . . . . . . . . . . . . . . 82
4.4.3 An Empirical Saturation Model . . . . . . . . . . . . . . . . . . 83
4.4.4 Comparison to Diffusion Limited Aggregation . . . . . . . . . . 93
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.7 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Table of contents xv
5 Conclusions 105
5.1 Main contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.2 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
References 111

List of figures
1.1 The displacement of stable displacement of oil by water and the unstable
displacement of water by oil in a Hele-Shaw cell five-spot pattern. From
Mungan (1971). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 The displacement of glycerol by air in a rectilinear Hele-Shaw cell at the
onset of instability (a) and after the development of a dominant finger
(b). From Saffman and Taylor (1958). . . . . . . . . . . . . . . . . . . . 10
2.1 Dual-interface variable definitions. The reservoir fluid (fluid 3) and
post-treatment fluid (fluid 1) are separated by an annulus of treatment
fluid (fluid 2). The leading and trailing interfacial base states are at
radii R2 and R1 respectively, with perturbation amplitudes Bn(t) and
An(t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Variation of the maximum value of P at which the labelled mode, n =
2, ..., 6, becomes absolutely stable as a function of V . For a given mode
n, then for values of P smaller than this maximum, all modes n′ such
that n′ < n are absolutely stable. . . . . . . . . . . . . . . . . . . . . . 19
xviii List of figures
2.3 Stability diagram for a radial displacement flow with V = 0.1, showing
the highest Γ(n) = Q∗n(R
∗
2)R
∗
2 for which the flow is stable to pertur-
bations of mode n (dashed lines) plotted as a function of the ratio of
the viscosity difference between the reservoir fluid and the treatment
fluid, to the viscosity difference between the post-treatment fluid and
the reservoir fluid. Solid line represents the locus of stability for all
modes Γ(nmin) as a function of P . If a more viscous displacing fluid is
used, higher flow rates are permitted. . . . . . . . . . . . . . . . . . . 21
2.4 The non-dimensional difference in viscosity across the leading interface,
P , at which the trailing (dashed) and leading (solid) interfaces become
absolutely stable to the given mode as a function of the bounding vis-
cosity ratio, V for the limit Rˆ → 0. In (a) we show which interface
is absolutely stable for mode 2. In regions (i) and (iii) the leading or
trailing interface alone is respectively absolutely stable, whereas in re-
gion (ii) both interfaces are absolutely stable and in region (iv) both
interfaces impose a flow rate limitation for stability. In (b) we show
how this extends to higher modes. . . . . . . . . . . . . . . . . . . . . 27
2.5 The critical values of the non-dimensional leading viscosity jump, P , for
absolute stability as a function of Rˆ for V = 0.3. In regions (ii) and (iii)
modes 2 and 3 are absolutely stable respectively. As time and hence Rˆ
increase, the annulus thins and these lines converge. As such, absolute
stability is transient for modes 2 and 3 for this value of V . . . . . . . . 28
List of figures xix
2.6 The critical values of the non-dimensional difference in viscosity across
the leading interface, P , at which the dual-interface system is absolutely
stable as a function of V for modes 2 (solid lines) and 3 (broken lines)
for three values of Rˆ. The Rˆ = 0 curves can be recognised from figure
2.4 with cusps at Vlower(n). As Rˆ increases, the annulus thins and the
stability tends to that of a single interface in a system in which the post-
treatment fluid displaces the reservoir fluid, with a vertical boundary
at Vupper(n). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
xx List of figures
2.7 Variation of the maximum flow rate, Γ = QR2, for stability of the sys-
tem to all modes (thick solid line) as a function of the non-dimensional
leading viscosity jump, P , for V = 0.1 and for two values of Rˆ, Rˆ =
0.1, 0.5. Also shown in the figures are a series of thin lines on which
Det(M) = 0 for each of modes 2,...,6. In (a), for n > 2, there are two
dashed lines for each mode: one of these lines has relatively small values
of Γ for small P and Γ then increases monotonically with P , whilst the
other has relatively small values of Γ for large P and Γ then decreases
monotonically with P . The point of intersection of these two solu-
tion branches corresponds to the point Det(M) = 0 and Tr(M) = 0.
However, for n = 2, there are two non-intersecting lines for which
Det(M) = 0 and these are shown as solid thin lines. The lower branch
corresponds to the maximum value of Γ for which all perturbations of
mode 2 are stable, while the upper branch corresponds to the minimum
value of Γ for which all perturbations of mode 2 are unstable. In (b), we
show solutions for Det(M ) = 0 for modes n = 2, ..., 5 as thin solid lines,
and for each such mode, there is a lower branch on which Tr(M ) < 0
and and upper branch on which Tr(M ) > 0, and these branches do not
intersect. Note, however, that owing to the finite vertical scale, the up-
per branch of mode 2 does not appear in the graph. Also shown are two
dashed thin lines, corresponding to Det(M) = 0 for mode 6. The value
of Γ on these two dashed lines increase and decrease monotonically with
P , respectively, and their point of intersection again corresponds to the
point Det(M) = 0 and Tr(M) = 0. . . . . . . . . . . . . . . . . . . . 33
2.8 The mode that requires the lowest flow rate to remain stable (nmin) as a
function of P for V = 0.1 and Rˆ = 0.1, 0.3, 0.5. As Rˆ increases coupling
causes lower modes to impose the strictest limit on flow rate. There are
some small spikes in nmin when the respective range of values of P is
narrow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
List of figures xxi
2.9 (a) The maximum stable Γ = QR2 as a function of P for a series of
increasing R2 where V = 0.1. The post-treatment fluid is injected at
R2 = 5. As such, for 1 ≤ R2 ≤ 5 the stability of the system is bounded
by the single-interface analysis (broken line). Subsequently, the stability
of the trailing interface must be considered and the stability is bounded
by the dual-interface analysis (solid lines). (b) The maximum Q for
stability varying with R2 for the cases A, B and C in (a). . . . . . . . 36
2.10 Time to sweep to R2 = 10 scaled with the time taken for a single
interface in which post-treatment fluid displaces the original reservoir
fluid for a series of increasing S as a function of the ratio of the leading
viscosity jump to the global viscosity jump. The optimum viscosity
choice balances the adverse viscosity ratios at the leading and trailing
interfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1 Left column: θ = 1 (slow gelling). Right column: θ = 10 (rapid gelling).
Top row: viscosity ratio profiles with parameters θ, V0 = 0.01 and
Vf = 0.9 (primary axis) and the optimal velocity profiles (secondary
axis). Middle row: the evolution of
√
τkmax for the optimal (dashed
lines) and constant (solid lines) flow rate profiles. Bottom row: the
evolution of
√
τσ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.2 (a) The final amplitude of perturbations of wavenumber k for the viscos-
ity invariant case (θ = 0, grey) and a gelling case (θ = 1) for a constant
flow rate (solid line) and the optimal flow rate (broken line). (b) The
final maximum amplitudes varying as a function of θ for Vf = 0.9. (c)
The ratio of the these two final amplitudes as a function of θ for a series
of Vf . (d) The effect of increasing V0 for a series of Vf and θ = 2. (e)
A contour plot of the final amplitude ratio varying with τ and θ for
V0 = 0.01 and Vf = 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . 50
xxii List of figures
3.3 The flow rate, Q∗(t), that minimises
∫ 1
0 σmaxdt for Λ = Λ
† (dashed line),
Λ = 3Λ† and Λ = 10Λ† (solid lines) and Λ = 100Λ† (dashed-dotted line)
while R0 = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.4 (a) The evolution of the upper bound of growth rate, σ∗max, for V = 0
(broken lines), V = 0.1 (solid lines), R0 = 0.2 and a series of values of
Λ
Λ†
that increases from the bottom of the figure upwards. (b) A regime
diagram for the conditions σ∗max < 0 ∀ t (broken lines; region S) and
σ∗max > 0 ∀ t (solid lines, region U) for a series of values of V increasing
from left to right (as labelled). . . . . . . . . . . . . . . . . . . . . . . 55
3.5 (a) The evolution of the most unstable mode nmax for R0 = 0.1 and a
series of Λ
Λ†
that increases from the bottom of the page upward. (b-g)
The growth rate of discrete modes n for a constant injection rate (2nd
row) and the optimal injection rate (3rd row) for V = 0, R0 = 0.1 and
Λ = 3Λ† (left column), Λ = 10Λ† (centre column) and Λ = 100Λ† (right
column). The mode number increases from the LHS to the RHS of
the figures. (h-j) The amplitude of perturbations following the optimal
injection rate. Initially, the lowest modes have the largest amplitude.
Mode 2 (11) is highlighted with thick line in orange (purple). The
growth-rate bound derived from treating n as a continuous variable is
shown as the thick black line. . . . . . . . . . . . . . . . . . . . . . . . 59
3.6 The final amplitude varying with Λ
Λ†
for V = 0, R0 = 0.1 for constant
injection (solid lines) and optimal injection (broken lines) computed
with discrete modes (black) and a continuous series of modes (grey). . 60
3.7 The critical value of Λ0, above which the system is always unstable for
the gelling case in which V0 = 0.3 and Vf = 0.5. A series of values of
gelling rate, θ, is shown. . . . . . . . . . . . . . . . . . . . . . . . . . . 61
List of figures xxiii
3.8 The optimal flow rate, Q∗(t), for a series of different gelling rates θ
for R0 = 0.35, V0 = 0, Λf = 1250 and Vf = 0.9 found using (3.27)
(solid lines) and the invariant viscosity case (broken line). In (a), the
timescale of viscosity change compared to the flow is small, whereas in
(b) it is large. In (a) the numerical solution to (3.19) is shown (crosses). 62
3.9 (a) The final amplitude of each mode n for the viscosity invariant case
(θ = 0, grey lines), for a constant flow rate (solid lines) and the equiva-
lent optimal flow rate (broken lines) where V0 = 0, Vf = 0.5, Λf = 1250
and R0 = 0.35. A gelling example is also shown (θ = 1, black). (b)
The variation of both maximum amplitudes with θ. (c) The ratio of
the constant to optimal injection final amplitude plotted against θ for
two values of Λf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.10 Evolution of amplitudes for a shear-thinning fluid in which the apparent
viscosity ratio varies according to V0 = 0.01, Vf = 0.9 and θ = 10
for constant injection (U = 1; solid lines) and a series of values of m
that decreases from the Newtonian case (m = 1; black) at bottom of
the figure upwards: m = 0.99 (blue), m = 0.95 (red) and m = 0.78
(yellow). Also included are the evolution of amplitudes in the case of
optimal injection (cf. figure 3.11(a) broken lines). . . . . . . . . . . . . 67
3.11 Injection with Vf = 0.9. (a) Optimal injection flow rate profiles. (b)
Variation of the final amplitude with gelling speed parameter θ for con-
stant injection (solid lines) and optimal injection (broken lines). The
values of m in the cases of shear thinning correspond to those in figure
3.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.1 The evolution of the pattern area, A. . . . . . . . . . . . . . . . . . . . 82
4.2 VF pattern development time series for the experiment keys and time
intervals, δt, as shown in the captions. Scales in [cm]. . . . . . . . . . . 84
xxiv List of figures
4.3 The evolution of saturation for three different examples of inlet con-
dition with the model (4.5) fitted (black lines). The time series are
evenly-spaced in time and are shown starting at an initial time ti at
intervals ∆t, with these values as shown in the captions as (ti, ∆t). . . 85
4.4 Schematic of the saturation profile. . . . . . . . . . . . . . . . . . . . . 85
4.5 The final saturation profile across experiments. In all figures, the broken
line is (4.4) with D = 1.70. . . . . . . . . . . . . . . . . . . . . . . . . 87
4.6 Variation of the fitting parameter Rb with C¯. The experiment key
number increases from left to right. . . . . . . . . . . . . . . . . . . . . 88
4.7 The evolution of the ratio Rtip/Rfr varying with Rtip/b for the injection
profiles as labelled (c.f. figure 4.5 for legend keys). . . . . . . . . . . . . 89
4.8 The evolution of the radii Rtip (left column) and Rfr (right column).
Top row: constant injection, middle row: constant pressure, bottom
row: ramped injection (c.f. figure 4.5 for legend keys; the experiment
key number increases from bottom left to top right). . . . . . . . . . . 92
4.9 (a) A DLA cluster of 105 particles. (b) The radial saturation profiles
of DLA clusters of size ranging from 105 to 106 particles in ten equal
increments. The simulation to produce each cluster size was repeated
seven times, and that shown is the mean saturation of all repeats. . . . 94
4.10 Evolution of (a) the ratio Rtip/Rfr, (b) Rtip and (c) Rfr as a function
of the number of particles released. . . . . . . . . . . . . . . . . . . . . 95
4.11 Comparison of the radii Rin (a) and Rtip (b) for the viscous fingering
experiments (points) and DLA (broken lines). . . . . . . . . . . . . . . 97
4.12 The evolution of Utip [cm/s] (left axis) and Ca (right axis). (a) Praud
and Swinney (2005) experiments (bottom to top: P1 to P3). (b) Con-
stant flow rate experiments C1-C3. . . . . . . . . . . . . . . . . . . . . 99
4.13 Saturation profiles for the final patterns of May and Maher (1989) (+),
Thomé et al. (1989) (o), and Lajeunesse and Couder (2000) (their fig.
17(a) - (boxes)) and (fig 17.(b) - (diamonds)). . . . . . . . . . . . . . . 100
List of figures xxv
4.14 The box counting algorithm applied to different regions of two final
VF patterns. From the top series to the bottom: Ov = Overall, Sc =
Screened (frozen) and Ac = Active (growth). . . . . . . . . . . . . . . . 101

List of tables
4.1 Experimental setup. Experiments C1-C3 and R1-R3 were performed by
the authors of this chapter, whereas P1-P3 were performed by Praud
and Swinney (2005) and C4 and R4 were performed by Bischofberger
et al. (2015). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.2 Measured area and flow rate. . . . . . . . . . . . . . . . . . . . . . . . . 81
4.3 The error in the frozen region, RMSEfr, and shape parameters of the
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.4 The image processing MATLAB functions used and a brief description
thereof. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Chapter 1
Introduction
The theme of this thesis is the control of viscous fingering, the need for which is outlined
in §1.1. In §1.2 we introduce the classic theoretical treatment of viscous fingering:
the linear stability analysis. We discuss previous work that has used linear stability
analysis as a theoretical tool to develop methods of controlling viscous fingering and
in doing so we outline chapters two and three of this thesis, which propose novel
control methods. Then, in §1.3, we describe new experiments to explore the nonlinear
development of viscous fingering and outline chapter four, in which we develop a
new empirical model to describe the patterns that form during this nonlinear growth.
This provides insight into the efficacy of the control methods developed in chapters 2
and 3 for the fully nonlinear growth.
1.1 Motivation
When one fluid is displaced by another in a porous medium the interface between
the fluids may be stable or unstable depending on their respective mobilities in the
medium. In a homogeneous medium, if a less viscous fluid displaces a more viscous
fluid then the interface is subject to the viscous fingering instability, whereby a series
of fingers of the less viscous fluid penetrates into the fluid being displaced.
2 Introduction
Viscous fingering can be problematic. In the oil and gas industry the process of
waterflooding involves the injection of water into the oil reservoir in order to drive the
oil towards the production well and to maintain the reservoir pressure. Upon water
reaching the production well, an event known as ‘breakthrough’, further injected water
may follow the path of least resistance along a channel of water, rather than continue
to displace the oil towards the production well. Viscous fingering may cause this
channel to form more readily, reducing the time until breakthrough and causing larger
regions of uncaptured oil to be left in the reservoir (Lake, 1989). This is illustrated
in figure 1.1, which shows the immiscible displacement of one fluid by another at the
time that the displacing fluid breaks through to the production well. The experiment
was carried out in a thin gap between two parallel plates, effectively behaving as
a two-dimensional flow. The flow geometry is the ‘five spot’ pattern and there is
no flow across boundaries. Figure 1.1(a) illustrates the displacement of water by a
viscous oil from the injection well to the production well. Since the oil is more viscous,
the displacement is not subject to viscous fingering and the majority of the water is
displaced out of the cell at the time of breakthrough. However, some water does remain
since the fastest flow is down the centreline. In contrast, figure 1.1(b) illustrates the
opposite case: the interface is unstable and fingers of water penetrate through the
oil zone. A large finger in the centreline of the flow extends to the production well,
causing breakthrough before the majority of the oil can be swept.
Other difficulties can arise from viscous fingering in the oil and gas industry. If the
rock surrounding the injection well, where the fluid pressure may be large, is poorly
consolidated then sand can be dislodged from the rock at this site. This sand may
then travel with the fluids to the production well and erode the well and associated
machinery. This process is called ‘sand production’ and can be very destructive. To
avoid sand production, a chemical treatment fluid can be deployed to consolidate the
rock close to the injection well. One drawback of chemical treatment is that it might
block pore throats in the immediate vicinity of the injection well (Paraskeva et al.,
2000). Thus, following the chemical treatment a second fluid, that cleans up the
1.1 Motivation 3
Injection
Production
oil
water
(a) Oil displacing water
Injection
Production
water
oil
(b) Water displacing oil
Fig. 1.1 The displacement of stable displacement of oil by water and the unstable
displacement of water by oil in a Hele-Shaw cell five-spot pattern. From Mungan
(1971).
well and displaces the chemical treatment away from the well, may be injected. The
chemical treatment then forms an annulus that travels through the reservoir. In the
presence of adverse viscosity ratios the annulus of fluid can break up due to viscous
fingering (Cardoso and Woods, 1995). Talaghat et al. (2009) describe such fingering
as a problem during sand consolidation that should be mitigated in the minimum
amount of time since the well is out of production during chemical treatment.
Research interest is not limited to the oil and gas industry. Other examples of
where viscous fingering is undesirable are catalyst regeneration in fixed bed reactors,
filtration processes and sugar refining (Hill, 1952; Homsy, 1987). In other contexts,
the instability may be exploited, for example to enhance fluid-fluid mixing (Jha et al.,
2011). Turning to natural phenomena, the degree of mixing of fluids of different
viscosity may be a factor that controls the composition of magma during its ascent in
the conduit of a volcano (McBirney, 1984). In biology, Lubkin and Murray (1995) used
a model from the field of viscous fingering to describe the early development of lung
tissue. More recently, Callan-Jones et al. (2008) noted the similarity of the growth of
a particular type of biological cell fragment and viscous fingering. Viscous fingering
is also used as an archetype of interfacial instability that leads to pattern formation
4 Introduction
(Casademunt, 2004). The subject is not short of applications and has proven to be a
rich area of scientific research over the past half century.
1.2 Linear Stability Analysis: Predicting Interfa-
cial Control
As a preface to introducing the scientific works upon which we have built, we will
first introduce a simplifying analogue to porous media flow. The flow of fluid in a
porous medium is governed by Darcy’s law, which states that the velocity field is
proportional to pressure gradients. This is also true of flow in the narrow gap between
two parallel plates. This experimental device is called a Hele-Shaw cell (Hele-Shaw,
1898). Since the gap is narrow, motion perpendicular to the plates may be neglected
and gap-averaged quantities are used. The equivalence is imperfect for two-phase flow
since the Hele-Shaw model neglects pore-scale dispersive capillary forces that affect
the propagation of the interface as well as other two-phase flow effects (Perkins and
Johnston, 1969). In a porous medium, surface tension only acts on the scale of the
pores. This force is not modelled in the Hele-Shaw geometry (Homsy, 1987). Therefore,
there is a limitation on the validity of the relation of our analysis to a porous medium
because the planform curvature we describe, which acts to select the fastest growing
wavelength of instability, would not be present. The conclusions reached in later
chapters are strictly only valid for Hele-Shaw flows, however, an effective capillary
pressure, which depends on the proportion of each fluid phase in pore, does exist
and acts to cut-off small wavelengths of instability in porous media (Woods, 2015).
Nonetheless, the Hele-Shaw cell is a useful simplifying tool that provides insight into
viscous fingering phenomena and is interesting in its own right (Homsy, 1987). We
will now introduce our research in the context of the literature, which is mostly based
on Hele-Shaw flow.
Whilst the criterion for stability was first derived by Hill (1952), the viscous fin-
gering instability is often referred to by the study that first derived the dispersion
1.2 Linear Stability Analysis: Predicting Interfacial Control 5
relation. Saffman and Taylor (1958) performed a linear stability analysis (LSA) of the
interface between two immiscible fluids of differing viscosity in a rectilinear geometry.
If the displaced fluid is more viscous then, according to Hele-Shaw flow, the pressure
gradient in the longitudinal direction is steeper. Thus, a forward perturbation of the
less viscous fluid is surrounded by lower pressure fluid and growth is encouraged. It
was shown by Saffman and Taylor (1958) that a fastest growing wavelength exists
because long wavelength perturbations are more stable and the large curvature of
short wavelength perturbations stabilizes those perturbations through the action of
interfacial tension, which was included in their study using a formulation proposed by
Chouke et al. (1959).
Following this work, Bataille (1968) extended the linear stability theory to the more
complex case of the flow of a less viscous fluid injected into a more viscous fluid from a
point source. The complexity arises from the evolution of the base state. Furthermore,
the wavenumber of the circular modes is quantised. At early stages in the flow the cur-
vature may be sufficiently large that interfacial tension stabilizes all azimuthal modes.
In the experiments presented by Paterson (1981) the fingers featured the azimuthal
wavelength predicted by linear stability analysis at the onset of instability. However,
Maxworthy (1989) found that this prediction of wavelength is only valid for small
capillary numbers, and for larger capillary numbers the wavelength that is selected is
larger than that predicted by Paterson (1981). Paterson (1981) adopted the simplest
formulation for the interfacial pressure jump, which was introduced by Chouke et al.
(1959) and Saffman and Taylor (1958). This condition only accounts for interfacial
tension as a cause for the pressure jump at the interface. In fact, Dias and Miranda
(2013b) have shown that an accurate prediction of the selected wavelength can be made
for all capillary numbers if the pressure jump condition also includes: i) the effect of a
thin wetting film that trails behind the interface (Park and Homsy, 1984) and ii) the
viscous stresses at the interface (Kim et al., 2009). These developments demonstrate
that linear stability theory can accurately predict the wavelength of instability at the
onset of instability. In order to simplify the analysis, we neglect the contribution to
6 Introduction
the pressure jump that originates with normal stresses and ignore the viscous film
that trails behind the interface. Strictly, the pressure jump condition we use is only
valid under conditions where the wavelength of instability is much greater than the
gap width of the Hele-Shaw cell and the displacing fluid is perfectly wetting. We will
now discuss the approaches to controlling the instability that have been developed,
which have involved: modifying the geometry of the system; modifying properties of
the fluids; and varying the injection flow rate.
It is possible to alter the stability by modifying properties of the Hele-Shaw cell
such that additional terms are introduced to the growth rate equation. Al-Housseiny
and Stone (2013) introduced a small tapering angle to the cell gap and demonstrated
that it is possible to stabilise an otherwise unstable flow with a convergent angle.
There was good agreement between their new linear stability analysis prediction of
the stability threshold and experimental realisations of the flow. In another approach,
replacing the top plate with an elastic film has been shown to delay the onset of
instability (Pihler-Puzović et al., 2012) since the perturbations in pressure ahead of
the advancing front are reduced by the bending of the film. Rather than introducing
new terms to the growth rate equation, a variety of approaches involve altering the
magnitude of viscous destabilisation relative to interfacial tension, and the gradual
stretching of the interface in the case of radial spreading, in order to select the fastest
growing mode for the duration of the flow.
Cardoso and Woods (1995) found the injection profile Q(t) = αct
−1/3 leads to the
total stabilisation of an axisymmetric displacement, where Q(t) is the time-evolving
injection flow rate and αc is the constant of proportionality such that mode 2 is selected
to be the fastest growing mode and is dynamically stabilised throughout injection by
the continuous matching of destabilisation and stabilisation terms. With the same time
dependency, Brenner et al. (1990) suggested that the injection profile Q(t) = αt−1/3
would give self-similar solutions for the nonlinear growth of a finger in a 90◦ sector
cell. Li et al. (2009) numerically and experimentally confirmed that for any particular
α > αc, the linear stability theory prediction of the fastest growing mode is selected
1.2 Linear Stability Analysis: Predicting Interfacial Control 7
by nonlinear interactions as that which becomes manifest after a brief adjustment
period. These nonlinear interactions were theoretically analysed by Dias and Miranda
(2010) who performed a mode-coupling analysis to derive the criterion that controls
the selection of the fastest growing mode. This work was extended by Zheng et al.
(2015) who showed that the fastest growing mode could be selected by introducing
time dependency to the gap width rather than the injection flow rate. For a constant
injection flow rate, they found that uniformly varying the plate separation such that
the gap width, b(t), followed the power law b ∝ t1/7 led to the fastest growing mode
predicted by linear stability analysis being fixed over time in the experiments. The
approach of interfacial control through mode selection has also been used in other
contexts, which we will now discuss.
The potential of a three fluid system in reducing adverse viscosity gradients has
also been explored (Mungan, 1971). The linear stability analysis for a three fluid
system was first introduced by Nayfeh (1972) in a rectilinear geometry. This was
later developed by Cardoso and Woods (1995) in an axisymmetric geometry. In the
latter case, the annular intermediate fluid was thin in radial extent and also of very
large viscosity. Experimentally, this resulted in the break up of the viscous annulus
into droplets, the number of which was correctly predicted by linear stability analysis.
Guided by the variation in injection rate required to stabilise a single interface, Cardoso
and Woods (1995) found the equivalent variation in injection rate needed to stabilise
the thin viscous annulus. In chapter two we generalise that analysis and develop an
equation governing the variation of growth rate of the instability as a function of time
for an annulus of any thickness or viscosity. We then go on to find the evolving critical
injection flow rate such that the system is dynamically stable and then show how to
select the viscosity of the annulus such that the overall injection time is minimised.
This addresses the industrial need of rapidly and uniformly treating an oil reservoir
with an annulus of chemical treatment fluid, as described in §1.1.
Following treatment of the reservoir to prevent sand production, oil recovery can
begin. In most situations, constraints on time may not permit injection with a flow rate
8 Introduction
that is sufficiently slow for stability and there will be growth of the instability. Often,
polymers are added to the displacing fluid to reduce adverse viscosity gradients and
improve stability (Sorbie, 1991). Mungan (1971) investigated the case of a three fluid
system in which the intermediate fluid was a polymer solution. He demonstrated that
programming the polymer concentration in a spatially varying manner leads to the very
effective suppression of instability. With the viscosity of the solution approximately
equal to the viscosity of oil at the leading interface and gradually reducing to that of
the water at the trailing interface, he showed a significant improvement in the recovery
of the oil phase compared to the recovery when the same mass of polymer is deployed
at a uniform concentration. This problem was then theoretically analysed by Gorell
and Homsy (1983), who arrived at the same conclusion through the variational calculus
framework they developed. However, this strategy may be operationally challenging to
deploy in the field. In chapter three we propose an alternative process: we consider
the case of a single interface displacement in which the displacing fluid uniformly
undergoes a gradual increase in viscosity, as might occur if a polymer is gradually
released from a soluble encapsulant (Gun and Routh, 2013), if the polymer undergoes
a change in viscosity that is thermally activated (Tran-viet et al., 2014) or if some
other reaction, such as cross-linking or polymerisation, gradually takes place. To
optimise this system, we build on the intriguing framework first proposed by Dias et al.
(2012). Rather than seeking an injection rate which stabilises a specific mode for the
duration of the injection, Dias et al. (2012) sought to minimise the final amplitude
of perturbations following the injection of a finite volume of fluid in a finite time. In
the case where the displacement started at a considerable radius from the injection
point and was unstable for the duration of the flow, Dias et al. (2012) found that the
optimal injection profile, Q(t), is linear, Q ∝ t, and leads to a significant suppression of
instability in both experimental and computational realisations of the flow. This has
been further confirmed in Hele-Shaw flow simulations by Huang and Chen (2015), who
found that in contrast to the immiscible case if the fluids are miscible then the fingering
is more severe with the linear injection scheme. Drawing on these results, we find that
1.3 Nonlinear Growth 9
in the case that the viscosity of the displacing fluid gradually increases, the optimal
flow rate profile involves the injection of more fluid later in the flow as compared to
the constant viscosity case in both rectilinear and axisymmetric displacements. Many
polymer solutions feature a reduction in the effective viscosity at larger shear rates on
account of the alignment of the polymer chains (Kamal et al., 2015). This property
is known as shear-thinning. We extend our analysis of a Newtonian displacing fluid
to a shear-thinning fluid and find that, in this case, the optimal flow rate involves
increasing the injection flow rate to a lesser extent as compared to the Newtonian
case.
These control methods are based on the linear stability analysis of the interface,
which is valid when the amplitudes of instability are modest. As the instability devel-
ops fully, nonlinear calculations are needed to describe the evolution of the interface.
In the most unstable and developed cases a heavily branched pattern forms for which
bulk properties, such as the fractal dimension, are often invoked to describe the pat-
tern.
1.3 Nonlinear Growth
To illustrate the development of viscous fingering beyond the linear regime, Saffman
and Taylor (1958) used air to displace glycerine in a rectilinear Hele-Shaw cell, as
reproduced in figure 1.2. Figure 1.2(a) illustrates the onset of instability and there
are seven fingers with similar shapes whereas figure 1.2(b) depicts a later stage of the
flow, when one finger has become dominant.
In nonlinear growth, such as that illustrated in figure 1.2(b), the injected fluid
is diverted from neighbouring fingers to supply the dominant fingers. We shall refer
to this interactive mechanism as screening (Praud and Swinney, 2005). Lajeunesse
and Couder (2000) classify the subsequent behaviour of the dominant finger as either
stable or unstable depending on its confinement, where the confinement is the ratio
of the most unstable wavelength according to linear stability analysis (Chouke et al.,
10 Introduction
air
glycerol
(a) Onset
air
glycerol
(b) Dominant finger
Fig. 1.2 The displacement of glycerol by air in a rectilinear Hele-Shaw cell at the onset
of instability (a) and after the development of a dominant finger (b). From Saffman
and Taylor (1958).
1959; Saffman and Taylor, 1958) to the width of the channel. If this wavelength
and the channel width are of a similar scale then confinement is strong and a stable
finger propagates through the channel. In contrast, if the wavelength is much smaller
than the channel width then the confinement is weak and the finger undergoes the
tip-splitting instability whereby the fingers continually branch.
To generalise the behaviour in rectilinear cells to those with non-parallel sides,
Thomé et al. (1989) investigated the case of fingers propagating through sector-shaped
cells and found self-similar fingers in the case of strong confinement. As in the case of
rectilinear geometry, the role of interfacial tension in selecting the angular width of the
self-similar fingers has been understood in sector geometry (Combescot and Ben Amar,
1991). Lajeunesse and Couder (2000) explored the growth of weakly confined fingers
from the rectilinear geometry to the full circular geometry. It was found that above a
90◦ sector angle two long-lived structures are seen to coexist, whereas for sector angles
smaller than this critical angle non-dominant branches are eventually screened from
growth. In the full circular geometry, Thomé et al. (1989) proposed that the fingers
behave as if they are divided into virtual sector cells with the virtual walls defining the
effective confinement, and therefore stability, of the fingers. Lajeunesse and Couder
1.3 Nonlinear Growth 11
(2000) showed that if the confinement is weak there will be many generations of the
tip-splitting instability. In this case, complex patterns form and bulk quantities, such
as the fractal dimension, are used to describe the fingering. We now briefly introduce a
number of key studies that investigated this weakly confined regime of viscous fingering,
which becomes the focus of our final chapter.
Experiments performed by Praud and Swinney (2005) involved air injected at
constant pressure, which led to the flow rate gradually increasing with time. Highly
branched, complex fingering patterns were formed and the fractal dimension was found
to be 1.70. In contrast to Praud and Swinney (2005), May and Maher (1989) performed
constant flow rate experiments. In some of their experiments the fractal dimension
was found to be 1.70, whereas in others was not. Since there is evidence to suggest that
modifying the injection flow rate can suppress instability in the absence of tip splitting
(Dias et al., 2012; Huang and Chen, 2015), it is of interest to understand whether the
nonlinear pattern is different depending on the injection rate as a function of time.
In chapter four we re-analyse the Praud and Swinney (2005) and May and Maher
(1989) data, which present tip-splitting patterns. We assess the fractal dimension by
analysing the radial distribution of air saturation (azimuthally averaged area fraction).
There is no data available in the literature of the pattern of finger growth with a linearly
increasing flow rate, so we performed a series of experiments for this situation. We
discover a number of universal features in the variation of saturation with radius. In
all cases, it is seen that inside of a particular radius the pattern is frozen. Growth
of the pattern is concentrated at the tips as has previously been observed (Daccord
et al., 1986; Praud and Swinney, 2005). However, the detailed delineation between
the inner frozen pattern and the outer growing pattern is new. In the frozen region,
the variation of saturation decreases with radius according to a power law, whereas in
the outer part of the pattern the saturation decreases linearly to zero at the front of
the pattern. Surprisingly, the overall radius of the pattern is found to be independent
of the time dependence of the source flux for a fixed volume of injected fluid. The
saturation distributions obtained in the experiments are also very similar to those
12 Introduction
obtained from the 2D numerical diffusion limited aggregation models that we perform
following Witten and Sander (1981) and Kuijpers et al. (2014). We conclude with a
description of the possible applications of this work.
Chapter 2
On the selection of viscosity to
suppress the Saffman-Taylor
instability in a radially spreading
annulus
2.1 Summary
We examine the stability of a system with two radially spreading fronts in a Hele-Shaw
cell in which the viscosity increases monotonically from the innermost to the outermost
fluid. The critical parameters are identified as the viscosity ratio of the inner and
outer fluids and the viscosity difference between the intermediate and outer fluids as
a fraction of the viscosity difference between the inner and outer fluids. There is a
minimum viscosity ratio of the inner and outer fluids above which, for each azimuthal
mode, the system is stable to perturbations of that mode at any flow rate. This
condition is directly analogous to the result for a single interface. Below this minimum
ratio, the system may be stable at any flow rate early in the flow. However, once the
inner radius reaches a critical fraction of the outer radius, this absolute stability ceases
to apply owing to the coupling of the inner and outer interfaces. We determine the
14 Annular Saffman-Taylor Suppression
maximum flow rate, as a function of time, in order that all modes remain stable due to
the effects of interfacial tension. These criteria for stability are then used to select the
viscosity of the intermediate fluid so that a fixed volume of the intermediate and then
inner fluid can be added to the system in the minimum time with the system remaining
stable throughout. The optimal viscosity for this intermediate fluid depends on the
relative volume of the inner and intermediate fluid and also on the overall viscosity
ratio of the innermost fluid and the original fluid in the cell, with the balance being to
suppress the early time instability of the outer interface and the late time instability of
the inner interface. We discuss application of this approach to a problem of injection
of treatment fluid in an oil well.
2.2 Introduction
In processes which involve production of fluids from a subsurface porous layer through
a well, interventions occur whereby from time to time two fluids are pumped in se-
quence into the porous rock from the well. One example of this occurs when reservoir
engineers inject a chemical treatment to consolidate the oil-bearing sands (Paraskeva
et al., 2000) or introduce an agent that inhibits the precipitation of salts in the event
of any mixing between geologic and injected water (Woods, 2015). After the chemical
treatment is injected into the rock from the well, post-treatment fluid is typically used
to clean up the well and as a result the finite volume of treatment fluid is displaced
by the post-treatment fluid, and forms an annulus some distance beyond the well. Ta-
laghat et al. (2009) describe how it is difficult to obtain a uniform front of the injected
chemical owing to the time constraint of injecting the treatment as quickly as possible,
since the well is unproductive during such interventions. Any break up or fingering
of the treatment fluid might lead to regions in which the treatment fluid is unable to
form a competent bond designed to prevent sand production.
The challenge for such injection arises where adverse viscosity ratios are present. As
well as two-phase flow effects, the interfaces are subject to the possible development
2.2 Introduction 15
of viscous fingering. In this study, in order to build fundamental understanding of
such coupled interface flows, we analyse the stability of an idealised flow involving two
unstable interfaces spreading radially from a central source in a Hele-Shaw geometry.
We use the analysis to select an optimum viscosity for the treatment fluid so as to
prevent any instability during the injection process while maximising the injection
rate.
Paterson (1981) pioneered the study of the single-interface radial viscous insta-
bility in a Hele-Shaw geometry, determining the linear growth rates of the different
azimuthal modes. This analysis has been developed in a number of directions. For
example, Miranda and Widom (1998) investigated weakly nonlinear tip-splitting phe-
nomena, whilst Parisio et al. (2001) explored viscous fingering on the surface of a
sphere. There has been interest in the control of such phenomena. Cardoso and
Woods (1995) explained how fingering could be prevented by keeping the interface
linearly stable at low flow rates. Li et al. (2009) developed a nonlinear analysis and
a technique to control the final shape of the interface. In the context of Hele-Shaw
flow, the possibility of control through varying the gap width has been explored (Al-
Housseiny and Stone, 2013; Dias and Miranda, 2013a).
The two-interface problem was first introduced by Nayfeh (1972) and later devel-
oped by Cardoso and Woods (1995) who performed a theoretical and experimental
dual-interface analysis in the special case in which the inner fluid is highly viscous.
With a stable trailing interface and an unstable leading interface, they modelled the
formation of drops from the annulus of intermediate fluid. However, in the present
problem of well treatment, both interfaces are likely to be unstable. Recently, Gin and
Daripa (2015) examined the growth rates of instabilities in a dual-interface system as
a function of the viscosity of the three fluids in the system. In that analysis less at-
tention was placed on the conditions for overall stability; however, this is relevant for
the injection of treatment fluid followed by a volume of post-treatment clean-up fluid,
and forms the focus of the present work. Our ultimate aim is to select the viscosity of
the intermediate (treatment) fluid so as to minimise the overall injection time for both
16 Annular Saffman-Taylor Suppression
the treatment fluid and the subsequent post-treatment fluid, whilst ensuring overall
stability of the system. We assume that the original reservoir (outer) fluid is more
viscous than the post-treatment (inner) fluid, and that the treatment fluid is of in-
termediate viscosity. For simplicity, we assume that the treatment fluid is immiscible
with both the reservoir fluid and the post-treatment fluid, so that there is interfacial
tension at both interfaces. In this regard, we note that if two fluids are only weakly
soluble and hence only partially miscible, this can also lead to an effective interfacial
stress (Korteweg, 1901; Pojman et al., 2006).
First, in §2.3 we analyse the linear stability of a single interface as a reference
calculation and identify that there is a maximum flow rate below which all modes are
stable. We relate this maximum flow rate to the viscosity of the displacing fluid. In
§2.4 we explore the dual-interface system and determine conditions for both absolute
stability of the lowest modes of the system, and also dynamic stability of all modes
through control of the flow rate. In §2.5 we use these results to find the optimal choice
of the viscosity of the treatment fluid so as to minimise the injection time.
2.3 Single-Interface Stability: Treatment Injection Phase 17
Fig. 2.1 Dual-interface variable definitions. The reservoir fluid (fluid 3) and post-
treatment fluid (fluid 1) are separated by an annulus of treatment fluid (fluid 2). The
leading and trailing interfacial base states are at radii R2 and R1 respectively, with
perturbation amplitudes Bn(t) and An(t).
2.3 Single-Interface Stability: Treatment Injection
Phase
2.3.1 Formulation
The flow configuration is shown in figure 2.1. In this section, we consider the stability
of the leading interface alone. Paterson (1981) showed that the linear growth of
perturbations of the form
bn = Bn(t)e
inθ (2.1)
during radial displacement of a single interface is given by
1
Bn
∂Bn
∂t
=
Qc23n
2πR22
− Q
2πR22
− κTn(n
2 − 1)
R32(µ2 + µ3)
. (2.2)
In (2.2), Bn is the amplitude of azimuthal mode n, Q is the volumetric flow rate per
unit depth, R2 is the radial position of the interface, µi is the viscosity of species
i, κ is the permeability of the porous medium (equivalent to h2/12 where h is the
Hele-Shaw cell gap spacing), T is the interfacial tension, the viscosity contrast cij =
(µj − µi)/(µj + µi), the treatment fluid is labelled 2 and the displaced reservoir fluid
18 Annular Saffman-Taylor Suppression
is labelled 3. The first term on the right-hand side represents viscous destabilisation
and the second and third terms represent the stabilising effects of stretching of the
interface and interfacial tension respectively.
2.3.2 Absolute Stability of mode n
For each mode n, equating the viscous destabilisation term with the stretching term
gives the maximum viscosity contrast for which the interface is stable for any flow
rate, at all radii
c23 ≤ 1/n. (2.3)
We refer to this as absolute stability. We envisage that the treatment fluid is of
intermediate viscosity between the reservoir fluid and the post-treatment fluid and so
introduce a parameter P , the ratio of the viscosity difference between the treatment
fluid and the reservoir fluid, to the viscosity difference between the post-treatment
fluid (labelled 1) and the reservoir fluid. We also introduce a viscosity ratio between
the post-treatment fluid and reservoir fluid V . They are respectively given by
P =
µ3 − µ2
µ3 − µ1 , V =
µ1
µ3
. (2.4)
Hence, we can rewrite (2.3) in terms of P and V :
P ≤ 2
(1− V )(n+ 1) . (2.5)
The values of P for absolute stability of mode n are shown in figure 2.2 for n = 2, ..., 6.
As V increases or P decreases higher modes become absolutely stable. Both regions
are labelled for mode 2, otherwise each region is labelled with the highest mode that
is absolutely stable.
2.3 Single-Interface Stability: Treatment Injection Phase 19
P
V
5
4
2 st’
2 unst’
6
3
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
Fig. 2.2 Variation of the maximum value of P at which the labelled mode, n = 2, ..., 6,
becomes absolutely stable as a function of V . For a given mode n, then for values of
P smaller than this maximum, all modes n′ such that n′ < n are absolutely stable.
2.3.3 Dynamic Stability
If the viscosity contrast is higher than that given by (2.3), the viscous destabilisation
term must be equated with the interfacial tension and stretching terms in order to
control mode n. This yields a maximum flow rate Qn(R2) for which mode n is stable
at a finite radius, R2,
Qn(R2) =
1
R2
2πκTn(n2 − 1)
(n(µ3 − µ2)− (µ3 + µ2)) . (2.6)
It is convenient to scale Qn(R2) by Qref , the flow rate at which mode 2 would go
unstable at the well radius (Rw) if the reservoir fluid were displaced by an inviscid
20 Annular Saffman-Taylor Suppression
species. Hence
Qref =
12piκT
Rwµ3
, (2.7)
R∗2 =
R2
Rw
, (2.8)
t∗ = tQref
piR2w
, (2.9)
Q∗n(R2) =
Qn(R2)
Qref
, (2.10)
Γ(n) = Q∗n(R
∗
2)R
∗
2 =
n(n2−1)
6((n+1)P (1−V )−2) . (2.11)
Henceforth, the stars are dropped for convenience. The minimum value of Γ(n),
Γ(nmin) say, identifies the highest dimensionless flow rate such that the interface is
stable to all such azimuthal modes (2.1). Higher modes are stabilised by interfacial
tension, and lower modes by stretching. As an approximation for large nmin, we can
find the solution n = nmin by finding a solution of
dΓ
dn = 0. This is given by the real,
positive root of the equation
n3(2P (1− V ))− n2(6− 3P (1− V ))− (P (V − 1) + 2) = 0. (2.12)
This provides a good approximation for large n, but is inexact since n is discrete. Note
that there is only one real positive root of (2.12) if P (1 − V ) < 2, which is always
true if the treatment fluid is of intermediate viscosity. In figure 2.3 we plot Γ(n) for
modes n = 2− 8 (dotted lines) and Γ(nmin) (solid line) as a function of P . The figure
illustrates how for smaller viscosity jumps, P , the interface remains stable at larger
flow rates Γ (2.11).
2.3.4 Single-Interface Control
The time t taken for the interface to advance to radius R2(t) starting with radius
R2(0) = 1 when injected at the highest dimensionless stable stable flow rate
Γ(nmin)
R2(t)
is
t =
∫ R2
1
2R′22
Γ(nmin)
dR′2 =
2
3
(R32 − 1)
Γ(nmin)
(2.13)
2.3 Single-Interface Stability: Treatment Injection Phase 21
 
 
Γ
P
Γ(2)Γ(8)
Γ(nmin)
Stable
Unstable
0 0.2 0.4 0.6 0.8 1
100
101
102
103
104
Fig. 2.3 Stability diagram for a radial displacement flow with V = 0.1, showing the
highest Γ(n) = Q∗n(R
∗
2)R
∗
2 for which the flow is stable to perturbations of mode n
(dashed lines) plotted as a function of the ratio of the viscosity difference between the
reservoir fluid and the treatment fluid, to the viscosity difference between the post-
treatment fluid and the reservoir fluid. Solid line represents the locus of stability for
all modes Γ(nmin) as a function of P . If a more viscous displacing fluid is used, higher
flow rates are permitted.
22 Annular Saffman-Taylor Suppression
and this implies that the dynamically stable flow rate, Q(R2), is
Q(R2) = Qn=nmin(R2) =
Γ(nmin)(
1 + 3
2
Γ(nmin)t
) 1
3
. (2.14)
In the limit t≫ 1
Q(R2) = Γ(nmin)
2
3
(
3
2
t
)− 1
3
. (2.15)
This generalises the result of Cardoso and Woods (1995) (their equation (5.2)) who
analysed the limiting case of a large viscosity contrast in which only mode 2 needs to
be considered.
This calculation identifies the maximum flow rate to ensure stability during the
period in which there is one interface. During this phase it is desired to make Γ(nmin)
as large as possible which can be achieved by making the injected fluid more viscous
(such that c23 → 0). However, by doing so the subsequent clean-up phase (in which
post-treatment fluid displaces the chemical treatment) tends to develop an unstable
interface. The destabilisation of the trailing interface in an annular system could lead
to the break up of the annulus of treatment fluid. Hence, we seek an optimum viscosity
to minimise the time the well is out of production for given volumes of treatment
chemical and post-treatment fluid. In order to inform this selection of viscosity we
model the stability of a two-interface system in §2.4.
2.4 Dual-Interface Stability: Post-Treatment Fluid
Injection Phase
2.4.1 Formulation
In this section we calculate the maximum flow rate for stability for a three fluid system,
in which there are two interfaces at positions R1(t)+a(θ, t) and R2(t)+b(θ, t), where a
and b are small perturbations to the shape of the interface. In general we can express
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase 23
the perturbations of the interfaces in the form of a power series:
a(θ, t) =
∞∑
n=1
An(t)e
inθ
b(θ, t) =
∞∑
n=1
Bn(t)e
inθ.
(2.16)
A schematic of the variables used in the calculation is shown in figure 2.1. The velocity
potential φ is related to velocity v and the pressure p by Darcy’s law for flow in porous
media
v = −κ
µ
∇p = −∇φ (2.17)
and the flow is incompressible hence the velocity potentials obey Laplace’s equation,
∇2φ = 0. The solution is assumed to comprise a steady solution φ0j = − Q2pi lnr + cj
and a first-order perturbation φ1j =
∑∞
n=1 φ
1
j,n(r, t)e
inθ where for mode n, the linearised
perturbation has radial structure φ1j,n given by Laplace’s equation (for j = 1, 2, 3)
φ11,n = αn(t)
(
r
R1
)n
,
φ12,n = βn(t)
(
r
R1
)−n
+ γn(t)
(
r
R2
)n
,
φ13,n = ǫn(t)
(
r
R2
)−n
.
(2.18)
Owing to the orthogonality of the modes, it follows that for each mode n at the trailing
interface, R1, the continuity of velocity and the jump in pressure owing to interfacial
tension are given by (Cardoso and Woods, 1995)
∂v01
∂r
An + v
1
1,n =
∂v02
∂r
An + v
1
2,n =
dAn
dt
; (2.19)
∂
∂r
(
φ01µ1
κ
)
An +
φ11,nµ1
κ
=
∂
∂r
(
φ02µ2
κ
)
An +
φ12,nµ2
κ
− T (1− n
2)
R21
An (2.20)
24 Annular Saffman-Taylor Suppression
where v1j,n = − ∂∂rφ1j,n. Similar conditions may be written for the leading interface at
R2. Substituting the velocity potentials (2.18) into the boundary conditions (2.19,
2.20) and eliminating αn(t), βn(t), γn(t) and ǫn(t) yields the following set of coupled
ordinary differential equations (ODEs) for each mode n:
dAn
dt
= f1
(
Q(n− f−11 )
2πR21
− Tκn(n
2 − 1)
R31(µ2 − µ1)
)
An︸ ︷︷ ︸
(i)
+ f2
(
Qn
2πR22
− Tκn(n
2 − 1)
R32(µ3 − µ2)
)
Bn︸ ︷︷ ︸
(ii)
dBn
dt
= f3
(
Qn
2πR21
− Tκn(n
2 − 1)
R31(µ2 − µ1)
)
An︸ ︷︷ ︸
(iii)
+ f4
(
Q(n− f−14 )
2πR22
− Tκn(n
2 − 1)
R32(µ3 − µ2)
)
Bn︸ ︷︷ ︸
(iv)
where
f1 =
c12(1− c23Rˆ2n)
1 + c12c23Rˆ2n
, f2 =
c23(1 + c12)Rˆ
n+1
1 + c12c23Rˆ2n
,
f3 =
c12(1− c23)Rˆn−1
1 + c12c23Rˆ2n
, f4 =
c23(1 + c12Rˆ
2n)
1 + c12c23Rˆ2n
,
Rˆ =
R1
R2
.
(2.21)
The terms labelled (i) and (iv) on the right-hand sides of (2.21) describe local be-
haviour in which we can recognise terms for viscous destabilisation, stretching of the
interface and interfacial tension (cf. (2.2)). The terms labelled (ii) and (iii) couple the
two interfaces and in the limit of a thick annulus (Rˆ→ 0) they become weak, leaving
two separate single interfaces. The coupled ODEs can be written in matrix form
d
dt
(
An
Bn
)
= M
(
An
Bn
)
. (2.22)
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase 25
The growth rates of the system (λ+, λ−) are the eigenvalues of M . In the limit of a
thin annulus (Rˆ→ 1; R1 = R2) the growth rates have the asymptotic form
λ+ =
Qc13n
2πR22
− Q
2πR22
− 2κTn(n
2 − 1)
R32
1
µ1 + µ3
,
λ− = − Q
2πR22
.
(2.23)
In equation (2.23) the largest growth rate is that which would occur with species 1
displacing species 3 (no annulus; cf. (2.2)), but with twice the individual interfacial
tension. In the limit of a stable inner interface (µ1 ≫ µ2, µ3) the result matches that
given by Cardoso and Woods (1995). These asymptotic limits were also derived by
Gin and Daripa (2015).
To understand the controls on the onset of instability, and hence determine the crit-
ical flow rate to ensure stability, it is of interest to study the nature of the instabilities
as the three controlling parameters Rˆ, V and P vary.
2.4.2 Absolute Stability of Mode n
The single-interface stability problem reveals the fascinating result that there is a
viscosity contrast below which a mode is stable at all flow rates and radii (absolute
stability), based on the balance between the stretching of the interface owing to the
radial spreading and the viscous destabilisation. We can draw from that result to
provide insight into the two-interface problem.
For small values of Rˆ, we expect the two interfaces to be largely decoupled. We
expect that as P increases the outer interface becomes progressively more unstable.
To assess this transition in more detail, it is useful to consider how the condition for
absolute stability depends on the viscosity change between the treatment and reservoir
fluids, P , and also the ratio between the post-treatment and reservoir fluid viscosity
V . We can then generalize these results for Rˆ = O(1) where coupling between the
interfaces becomes significant.
26 Annular Saffman-Taylor Suppression
Absolute stability results from balancing viscous destabilisation with stretching,
and in the dual-interface system this is done by setting Det(M ) = 0 with T = 0. This
yields
c12c23(1− Rˆ2n)n2 − (c12 + c23)n+ 1 + c12c23Rˆ2n = 0. (2.24)
In the limit Rˆ → 0, (2.24) gives the condition for absolute stability of mode n from
§2.3.2 for two separate single interfaces,
c12 ≤ 1
n
,
c23 ≤ 1
n
.
(2.25)
These conditions are combined by writing them in terms of P , giving the condition
for absolute stability for mode n
1− n+1
n−1V
1− V ≤ P ≤
2
(1− V )(n+ 1) . (2.26)
Below the lower bound the trailing interface is locally not absolutely stable, and like-
wise above the upper bound the leading interface is not. Equating these conditions
gives a critical V below which neither mode can be absolutely stable
Vlower(n) =
(
n− 1
n+ 1
)2
. (2.27)
In figure 2.4(a) we show the upper P bound of the absolutely stable region as a solid
line (familiar from §2.3.2) and the lower bound P as a dashed line for mode 2 as a
function of V . The system is only absolutely stable in region (ii). Elsewhere, viscous
destabilisation on one or both interfaces is possible and can only be stabilised by
controlling flow rate. In figure 2.4(b) we extend this to higher modes by indicating
which modes (n′ ≤ n) are absolutely stable in each domain. This figure is important
since it illustrates how, at the onset of instability, the lowest modes in the system may
not be unstable at all, depending on the parameters P and V . To determine which
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase 27
P
V
(iv)
Neither
Trailing
Both
Leading
(iii)
(ii)
(i)
0 0.1 0.2 0.3
0.4
0.6
0.8
1
(a) mode 2
P
V
5 6
3
2
4
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
(b) modes 2-6
Fig. 2.4 The non-dimensional difference in viscosity across the leading interface, P , at
which the trailing (dashed) and leading (solid) interfaces become absolutely stable to
the given mode as a function of the bounding viscosity ratio, V for the limit Rˆ→ 0. In
(a) we show which interface is absolutely stable for mode 2. In regions (i) and (iii) the
leading or trailing interface alone is respectively absolutely stable, whereas in region
(ii) both interfaces are absolutely stable and in region (iv) both interfaces impose a
flow rate limitation for stability. In (b) we show how this extends to higher modes.
mode will be unstable first, the effect of interfacial tension must be considered (see
§2.4.3).
As fluid continues to be supplied to the system, the inner radius grows and Rˆ
increases. This leads to a progressively increasing level of coupling between the leading
and trailing interfaces. In the case V > Vlower(n = 2), we identified values of P for
which mode 2 becomes absolutely stable with Rˆ ≪ 1; as Rˆ increases, we expect that
the coupling will cause these bounding values of P to converge, and eventually that
there is a value of Rˆ for which there is no longer a region of absolute stability for
mode 2. The bounding values of P can be found by solving (2.24) at chosen values of
V, Rˆ and n. We illustrate this trend in figure 2.5 in which we show how the bounding
values of P vary with Rˆ for mode 2 where V = 0.3. The boundary for mode 3 is also
plotted and is contained within that of mode 2. In this case, there are no regions of
absolute stability for mode 4 since 0.3 < Vlower(n = 4).
28 Annular Saffman-Taylor Suppression
P
Rˆ
(ii)
(iii)
0 0.2 0.4 0.6 0.8 1
0.1
0.3
0.5
0.7
0.9
1
Fig. 2.5 The critical values of the non-dimensional leading viscosity jump, P , for
absolute stability as a function of Rˆ for V = 0.3. In regions (ii) and (iii) modes 2 and
3 are absolutely stable respectively. As time and hence Rˆ increase, the annulus thins
and these lines converge. As such, absolute stability is transient for modes 2 and 3 for
this value of V .
In the limit of a thin annulus, Rˆ→ 1, (2.24) has the asymptotic form
c13 ≤ 1
n
. (2.28)
This condition can be obtained by considering a single interface with with post-
treatment fluid displacing reservoir fluid directly (see §2.3.2). It can also be written
as a function of V
Vupper(n) =
n− 1
n+ 1
. (2.29)
Equation (2.29) shows that in the case V > 1/3, the instability that would result
from post-treatment fluid directly displacing the reservoir fluid is only unstable to
higher modes, n ≥ 3 , and so the two-interface problem is also only unstable to these
higher modes. In this case, depending on the value of P we expect that the most
unstable mode could correspond to much higher values of n. Figure 2.6 shows how
the absolutely stable region of figure 2.4 evolves as Rˆ increases for mode 2 (solid lines)
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase 29
 
 
Rˆ = 1
Rˆ = 0.7
Rˆ = 0
P
V
Vlower(2)
Vupper(2)
0 0.1 0.2 0.3 0.4 0.5 0.6
0
0.2
0.4
0.6
0.8
1
Fig. 2.6 The critical values of the non-dimensional difference in viscosity across the
leading interface, P , at which the dual-interface system is absolutely stable as a func-
tion of V for modes 2 (solid lines) and 3 (broken lines) for three values of Rˆ. The
Rˆ = 0 curves can be recognised from figure 2.4 with cusps at Vlower(n). As Rˆ increases,
the annulus thins and the stability tends to that of a single interface in a system in
which the post-treatment fluid displaces the reservoir fluid, with a vertical boundary
at Vupper(n).
and mode 3 (broken lines). As Rˆ increases, the bounds tend towards a vertical line
at Vupper(n). The system is only absolutely stable if an interface between the post-
treatment and reservoir fluid is absolutely stable. The cusp at Vlower(n) for mode 3 is
larger than the cusp for mode 2, as expected from figure 2.5.
We have found that if Vlower(n) < V < Vupper(n) there are solutions of (2.24) that
are transient absolutely stable bounds of P , i.e. initially mode n can be absolutely
stable but later in the flow interfacial tension must be considered for dynamic stability.
Below this range there is no absolute stability and above it mode n remains absolutely
stable at all radii. If mode n is absolutely stable, we must consider the dynamic
stability of a higher mode to find the maximum flow rate for the system’s stability.
30 Annular Saffman-Taylor Suppression
2.4.3 Dynamic Stability
In general, the condition for stability is found by balancing the viscous destabilisation
with the stretching and interfacial tension terms to ensure λ+ = 0, since if either
growth rate were positive the perturbations would grow. The eigenvalues are given by
λ± =
Tr(M)
2
±
√√√√(Tr(M)2
4
−Det(M )
)
. (2.30)
It can be seen from (2.30) that if Tr(M ) ≤ 0 and Det(M) = 0, then λ+ = 0 and
λ− < 0. If we solve the relation Det(M ) = 0 we find two solutions for Qn(R), but
only one of these is in a region where Tr(M) ≤ 0 and therefore represents stability.
To illustrate the evolution of the stability of the different modes as the fluids invade
the pore space and the trailing front catches up with the leading front, in figure 2.7 we
present the maximum value of Γ = QR2 (thick line) below which the system is stable
for the cases Rˆ = 0.1 and Rˆ = 0.5 where V = 0.1.
The left hand panel corresponds to Rˆ = 0.1, which occurs at an early time in the
injection process when the two interfaces are far apart. In this limit then to leading
order the stability of each interface may be approximated by the local stability of
that interface (2.21). In this case, the thick line showing the maximum flow rate for
stability is approximately given by the mode 2 stability bound of the trailing interface
for small P , but as P increases, this intersects the mode 6 stability bound of the
leading interface for this value of Rˆ. As P continues to increase, the critical flow rate
for stability is given by successively lower modes, until coinciding, approximately, with
the mode 2 stability of the leading interface.
In this limit of small Rˆ, since the interfaces are far apart, the stability curves
Det(M ) = 0 for the higher modes are well approximated by the stability curves for that
mode on each of the two interfaces individually, and the effects of the coupling between
the interfaces are small. Indeed, in figure 2.7(a) for all the higher modes, n > 2, with
Rˆ = 0.1 there is a common point at which Det(M ) = 0 and Tr(M ) = 0 which
connects the solutions Det(M ) = 0 with Tr(M) < 0 and the solutions Det(M) = 0
2.4 Dual-Interface Stability: Post-Treatment Fluid Injection Phase 31
with Tr(M ) > 0. For a given mode, n > 2, as we move across the point Det(M ) = 0
and Tr(M) = 0, on the line Det(M) = 0 we can interpret this point as corresponding
to that at which the stability of the mode changes from being dominated by one
interface to the other. The main exception to this for Rˆ = 0.1 is for the lowest mode,
mode 2 (solid thin line). Here, it is seen that for mode 2 the solution of the relation
Det(M) = 0 does not intersect the solution Tr(M ) = 0; instead for mode 2 two
non-intersecting branches of the solution Det(M ) = 0 may be seen. Tr(M ) < 0
on the lower branch and Tr(M) > 0 on the upper branch. The coupling between
the interfaces leads to a loss of the special solution, Det(M) = 0 and Tr(M ) = 0,
for the lowest mode, mode 2. Instead there is a smooth adjustment along the line
Det(M) = 0 from the solution dominated by the leading interface when P ≈ 1 to the
trailing interface when P ≈ 0.
However, as more fluid is injected and the trailing interface catches up with the
leading interface, Rˆ = 0.5 (right hand panel) then the distance between the interfaces
becomes smaller and so some of the higher modes, notably modes 3-5, in the case
Rˆ = 0.5, also become highly coupled. As a result, for each of these modes there is
no longer a solution Det(M ) = 0 with Tr(M) = 0. Again, two independent solution
branches for Det(M) = 0 emerge for each of these modes. On the lower branch in
(P,Γ) space, Tr(M) < 0, and on a separate higher solution branch Det(M) = 0 we
have Tr(M ) > 0. Now the lower solution branch determines the maximum flow rate,
Γ, for which that mode is stable to all perturbations.
In figure 2.7(b), the dark line again shows the overall stability boundary for all
modes in terms of the maximum values of Γ for which the modes are stable, as a
function of P . For small P , the stability threshold is determined by mode 2. However,
for this larger value of Rˆ, Rˆ = 0.5, figure 2.7(b) shows that, as P increases, this first
intersects the stability boundary for the higher modes with the mode 4 stability curve.
With further increases in the value of P , this overall stability threshold is given by
mode 3 and then back to the mode 2 at very large P .
32 Annular Saffman-Taylor Suppression
It is also seen from the figure that higher modes only become unstable for large
values of Γ and so are not rate-limiting in determining the maximum flow rate for
stability. Indeed, in figure 2.8 we illustrate the most unstable mode for three repre-
sentative values of Rˆ as a function of P . In accord with the form of figure 2.7, it is
seen that as the value of P increases from close to 0, the most unstable mode jumps
from mode 2 to mode 6 when P has value of order 0.45, but then as P continues to
increase this most unstable mode gradually falls again to mode 2.
This figure, in conjunction with figure 2.7, demonstrates how lower modes impose
stricter limitations on the maximum stable flow rate as Rˆ increases.
Our analysis has identified that, for small Rˆ, there are values of P and V for which
the system is unstable only to higher modes. If V > Vlower (2.27), all modes n
′ < n
can be absolutely stable for a range of P . As Rˆ increases, these lower modes also
become progressively unstable. If V > Vupper(n) (2.29), all modes n
′ < n are stable
for all P and Rˆ. As fluid continues to be injected, the effect of interfacial tension,
which is the mechanism ensuring stability, becomes progressively weaker. As a result,
in order to ensure the stability of the system, the injection rate should progressively
decrease. Depending on the values of P and V , the mode imposing the strictest limit
on flow rate to maintain stability may become lower with increasing Rˆ. Owing to
the transition in this strictest-limit mode with the injection of a finite volume of fluid,
the minimum time required for injection whilst maintaining stability depends on the
choice of P .
2.5 Control Strategy
The aim of this section is to identify the optimal value of viscosity of chemical treat-
ment to minimise injection time. This will depend on the viscosity ratio of the post-
treatment fluid and reservoir fluid, V , but also on the volumes of chemical treatment
fluid and post-treatment fluid to be deployed.
2.5 Control Strategy 33
Γ
P
6
Stable
5
4
3
2
2
6
0 0.2 0.4 0.6 0.8 1
100
101
102
103
(a) Rˆ = 0.1
Γ
P
4
3
22
Tr(M) > 0
6
0 0.2 0.4 0.6 0.8 1
100
101
102
103
(b) Rˆ = 0.5
Fig. 2.7 Variation of the maximum flow rate, Γ = QR2, for stability of the system
to all modes (thick solid line) as a function of the non-dimensional leading viscosity
jump, P , for V = 0.1 and for two values of Rˆ, Rˆ = 0.1, 0.5. Also shown in the figures
are a series of thin lines on which Det(M ) = 0 for each of modes 2,...,6. In (a), for
n > 2, there are two dashed lines for each mode: one of these lines has relatively small
values of Γ for small P and Γ then increases monotonically with P , whilst the other
has relatively small values of Γ for large P and Γ then decreases monotonically with
P . The point of intersection of these two solution branches corresponds to the point
Det(M) = 0 and Tr(M) = 0. However, for n = 2, there are two non-intersecting
lines for which Det(M ) = 0 and these are shown as solid thin lines. The lower branch
corresponds to the maximum value of Γ for which all perturbations of mode 2 are
stable, while the upper branch corresponds to the minimum value of Γ for which all
perturbations of mode 2 are unstable. In (b), we show solutions for Det(M) = 0
for modes n = 2, ..., 5 as thin solid lines, and for each such mode, there is a lower
branch on which Tr(M) < 0 and and upper branch on which Tr(M ) > 0, and these
branches do not intersect. Note, however, that owing to the finite vertical scale, the
upper branch of mode 2 does not appear in the graph. Also shown are two dashed
thin lines, corresponding to Det(M ) = 0 for mode 6. The value of Γ on these two
dashed lines increase and decrease monotonically with P , respectively, and their point
of intersection again corresponds to the point Det(M ) = 0 and Tr(M) = 0.
34 Annular Saffman-Taylor Suppression
 
 
Rˆ = 0.5
Rˆ = 0.3
Rˆ = 0.1
P
n
m
in
0 0.2 0.4 0.6 0.8 1
0
2
4
6
8
Fig. 2.8 The mode that requires the lowest flow rate to remain stable (nmin) as a
function of P for V = 0.1 and Rˆ = 0.1, 0.3, 0.5. As Rˆ increases coupling causes lower
modes to impose the strictest limit on flow rate. There are some small spikes in nmin
when the respective range of values of P is narrow.
During treatment, a finite volume of fluid is injected with behaviour characterised
in §2.3, followed by a finite volume of post-treatment fluid (see §2.4). Using the
analysis of §2.3 and §2.4 we can identify the maximum possible flow rate for each
value of P , as a function of time during the injection. By integrating this flow rate
with time, we can then assess the value of P which minimises the total injection time.
As an example of such a calculation, in figure 2.9(a) the maximum Γ = QR2 for
stability is shown as a function of P for different values of the outer radius of the
treatment fluid R2. In these calculations, the post-treatment fluid is added to the
system when R2 = 5, and the whole process is completed when R2 = 10. During the
treatment phase, 1 < R2 < 5, the maximum flow rate to ensure stability of the system
may be found using the single-interface analysis of §2.3 (dotted line). However, once
the post-treatment fluid starts to be injected, for R2 ≥ 5, the stability of the trailing
interface should also be taken into account in determining the maximum flow rate for
stability of the system. We show how this maximum flow rate changes as a function
of R2 in figure 2.9(a), illustrating the bounds for the values R2 = 5.2, 5.5, 6.5 and 10
(solid lines) as obtained from the analysis of §2.4. For a specific choice of P we see that
Γ needs to be systematically reduced once R2 > 5 and both interfaces are migrating
2.5 Control Strategy 35
through the system, and this has been illustrated schematically with the lines A-A’ ,
B-B’ and C-C’.
For clarity, in figure 2.9(b) we illustrate how Q = Γ/R2 varies as R2 increases,
following each of the lines A-A’, B-B’ and C-C’ from figure 2.9(a). It is seen that
for larger values of P (e.g. curve C-C’), the maximum stable flow rate during the
treatment phase, R2 < 5, is smaller but as R2 increases further during the clean-up
phase, and the trailing interface becomes rate-limiting, the maximum flow rate may
be smaller for smaller values of P . The trade-off between these early and late time
effects can lead to a choice of P that minimises the overall time taken to inject the
treatment fluid.
In order to illustrate this trade-off, in figure 2.10 we show the total injection time
required for the leading interface to reach the value R2 = 10 as a function of the
viscosity of the treatment fluid, as parametrised by P . Curves have been shown for
different values of the ratio of volume of treatment fluid to the total volume of fluid
injected, as parametrised by S = 1 − Rˆ2. The overall time taken to sweep is found
by integrating Q(R2). For each case, there is an overall optimum viscosity for the
treatment fluid. Note that in the figure, time is scaled with the time it would take
to inject the same total volume of fluid with a single interface between the reservoir
fluid and the post-treatment fluid (cf. (2.13)). As S becomes larger, the optimal value
of P decreases; this is because the effect of the instability at the trailing interface
becomes smaller since less post-treatment fluid is added, and the interfacial tension
on this trailing interface is more effective for smaller radii. Therefore, there is benefit
in preferentially increasing the viscosity of the treatment fluid to maximise the average
injection rate.
As an illustration of the potential benefits of this analysis in terms of injection
time, we note that for example, in the case S = 0.25, by increasing the viscosity of
the treatment fluid to the optimal viscosity, the minimum injection time to maintain
overall stability is about 30% of the injection time when the treatment fluid has the
same viscosity as the post-treatment fluid.
36 Annular Saffman-Taylor Suppression
 
 
Γ
P
C
B
A
5.5
R = 5.2
10
1 < R < 5
6.5
C’B’A’
0 0.2 0.4 0.6 0.8 1
100
101
102
(a)
 
 
C - C’
B - B’
A - A’
A
B
C
Q
R2
2 4 6 8 10
10−1
100
101
102
(b)
Fig. 2.9 (a) The maximum stable Γ = QR2 as a function of P for a series of increasing
R2 where V = 0.1. The post-treatment fluid is injected at R2 = 5. As such, for
1 ≤ R2 ≤ 5 the stability of the system is bounded by the single-interface analysis
(broken line). Subsequently, the stability of the trailing interface must be considered
and the stability is bounded by the dual-interface analysis (solid lines). (b) The
maximum Q for stability varying with R2 for the cases A, B and C in (a).
2.6 Conclusions 37
 
 
S = 0.5
S = 0.25
S = 0.175
τ
P
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
Fig. 2.10 Time to sweep to R2 = 10 scaled with the time taken for a single interface in
which post-treatment fluid displaces the original reservoir fluid for a series of increasing
S as a function of the ratio of the leading viscosity jump to the global viscosity jump.
The optimum viscosity choice balances the adverse viscosity ratios at the leading and
trailing interfaces.
2.6 Conclusions
We have considered the stability of an annulus of fluid spreading from a point source
in a Hele-Shaw geometry, in which there is an increase in viscosity across both the
leading and trailing interfaces, so that both interfaces are potentially unstable. In
a reference calculation of a single interface spreading axisymmetrically from a point
source, we show that the stretching of the interface associated with the radial flow
permits absolute stability for lower azimuthal modes for sufficiently small viscosity
contrasts. We then find a maximum injection rate so that the viscous destabilisation
of the interface is suppressed by the combination of the stretching and the interfacial
tension for all higher modes.
We build on this analysis for a dual-interface system, in which there may be in-
teractions between instabilities at the leading and trailing interfaces. We find that
there is a range of viscosity of the intermediate fluid such that both interfaces are
38 Annular Saffman-Taylor Suppression
absolutely stable and show how this range depends on the ratio of viscosities of the
bounding fluids. We also find the maximum injection rate such that the combination
of stretching and interfacial tension is able to suppress the instability of all modes
in the system. By following the evolution of this critical flow rate with time during
the injection of an annulus of finite volume and a finite volume of the subsequent
fluid inside the annulus, we can calculate the minimum total time for injection. We
examine how this minimum time varies with the viscosity of the fluid in the annulus,
and thereby identify an optimal choice of viscosity that minimises the injection time
while maintaining overall stability of the system. We also show how as the volume of
the annulus as a fraction of the total volume of fluid injected increases, the optimal
viscosity of the annular fluid increases. Also, the reduction in the time required for the
injection of the annular fluid and the following inner fluid relative to the time to inject
an equal volume of the inner fluid becomes progressively larger in order to maintain
stability.
Although our analysis is strictly valid for Hele-Shaw geometries, the analysis sug-
gests that in field operations in which treatment chemicals are added to a production
well for scale management or to suppress sand production, there may be benefit in
making the treatment fluid more viscous in order to minimise the time required for
the treatment but also to ensure a uniform distribution of the treatment fluid around
the well. As a very simplified example, if 10m3 of treatment fluid was added to a well
followed by 10m3 of post-treatment fluid, with an injection rate of 0.0001 m3/s, the
treatment time would be 30-40 hours; by making the treatment fluid more viscous,
this process may be accelerated to a time closer to 10-20 hours, substantially reducing
the non-productive time of the well.
Following the chemical treatment of the well, oil production may begin. During
secondary recovery, polymers may be added to the injected water in order to reduce
adverse viscosity gradients. The potential of deploying a polymer solution that features
time-dependent rheology will now be explored.
Chapter 3
Control of viscous instability by
variation of injection rate in a fluid
with time-dependent rheology
3.1 Summary
Using variational calculus, we investigate the time-dependent injection rate which
minimises the growth of the Saffman-Taylor instability when a finite volume of fluid is
injected in a finite time, tf , into a Hele-Shaw cell. We first consider a planar interface,
and show that with a constant viscosity ratio the constant injection rate is optimal.
When the viscosity of the displacing fluid, µ1(t), gradually increases over time, as
may occur with a slowly gelling polymer solution, the optimal injection rate, U∗(t),
involves a gradual increase in the flow rate with time. This leads to a smaller initial
value of flow rate than the constant injection rate, finishing with a larger value. Such
optimisation can lead to a substantial suppression of the instability as compared to
the constant injection case if the characteristic gelling time is comparable to tf . In
contrast, for either relatively slow or fast gelling there is much less benefit in selecting
the optimal injection rate, U∗(t), as compared to the constant injection rate. In the
case of a constant injection rate from a point source, Q, with a constant viscosity ratio
40 Viscous instability with time-dependent rheology
the fastest growing perturbation on the radially spreading front involves axisymmetric
modes whose wavenumber increases with time. Approximating the discrete azimuthal
modes by a continuous distribution, we find the injection rate that minimises growth,
Q∗(t). We find that there is a critical time for injection, t†f , such that if tf > t
†
f then
Q∗(t) can be chosen so that the interface is always stable. This critical time emerges
from the case with an injection rate given by Q∗ ∼ t−1/3. As the total injection
time is reduced to values tf < t
†
f , the system becomes progressively more unstable
and the optimal injection rate for an idealised continuous distribution of azimuthal
modes asymptotes to a flow rate which increases linearly with time. As for the one-
dimensional case, if the viscosity of the injection fluid gradually increases over time,
then the optimal injection rate has a smaller initial value but gradually increases to
larger values than for the analogous constant viscosity problem. If the displacing fluid
features shear-thinning rheology, then the optimal injection rate involves a smaller flow
rate at early times, although not as large a reduction as in the Newtonian case, and
a larger flow rate at late times, although not as large an increase as in the Newtonian
case.
3.2 Introduction
Viscous fingering occurs when a low viscosity fluid is used to displace a more viscous
fluid through a porous medium. The phenomenon is of considerable importance for the
oil industry since it can lead to injected water bypassing the oil in a reservoir, with
the result that there are large pockets of unswept oil (Lake (1989); Woods (2015)).
The injection of a less viscous fluid into a more viscous fluid in a Hele-Shaw cell is also
subject to viscous fingering. The Hele-Shaw cell is used as a two-dimensional analogue
model of a porous medium where pore-scale physics is neglected. Also, the planform
curvature which acts with surface tension to cut off the wavenumber of instability in a
Hele-Shaw cell is not present in a porous medium. With these problems acknowledged,
we focus on Hele-Shaw flow in order to develop new concepts.
3.2 Introduction 41
A variety of approaches to control the stability of the interface have been explored.
They can be classed as involving modification of (1) the Hele-Shaw cell geometry,
(2) the injection flow rate, (3) fluid properties or combinations thereof. For the first
class, Al-Housseiny and Stone (2013) performed a linear stability analysis for flow in a
tapering cell, and experimentally verified how the tapering angle modified the stability
of the interface whereas Pihler-Puzović et al. (2012) investigated the stabilising effect
of replacing the top plate with an elastic membrane. Zheng et al. (2015) defined a
control parameter that, if held constant in time, permits selection of which azimuthal
mode is manifest throughout the injection. In that study, the control parameter was
held constant by varying the gap thickness, b, according to b(t) ∼ t1/7. In a different
study, the control parameter was held constant by varying the injection flow rate, Q,
according to Q(t) ∼ t−1/3 (Li et al., 2009). Both studies demonstrated good agreement
between experiment and theoretical prediction of the dominant mode of instability
which develops. For the latter case, Dias and Miranda (2010) went on to perform a
weakly nonlinear mode-coupling analysis to explain how the Q ∼ t−1/3 injection flow
rate leads to sharpening, stable fingers.
Modification of the fluids (class three) has also been shown to give interfacial con-
trol. The viscosity of the injected fluid can be increased, for example by adding a
polymer (Sorbie, 1991). Since the field-scale deployment of polymer can be expensive,
in some cases a finite volume of polymer solution is injected between the oil and the
water. This leads to a series of interesting problems concerning the stability of a flood
front with two interfaces. Gorell and Homsy (1983) explored the stability of a uni-
directional constant flow in which the concentration of injected polymer varies with
position in the flow. They developed a variational approach to determine the optimal
concentration of the injected polymer as a function of time in order to to minimise the
viscous instability when a finite mass of polymer and water are injected. Following a
different approach, a stability analysis has been presented for a coupled problem based
on flow in a Hele-Shaw cell (Cardoso and Woods, 1995). Gin and Daripa (2015) devel-
oped this by calculating the stability of a series of discrete layers of different viscosity
42 Viscous instability with time-dependent rheology
injected in sequence. However, there are challenges associated with the direct injection
of polymer solutions into a porous medium, including the risk that the polymer may
block up pore throats near the injector and hence prevent the continued injection of
fluid into the reservoir in order to displace the oil (Sorbie (1991); Woods (2015)). In
order to mitigate such risks, it is possible to use a polymer with a delayed activation so
that the mixture viscosity gradually increases with time. There are a variety of meth-
ods that could achieve this. For example, if the polymer were supplied in a soluble
encapsulant with a finite release time, then the concentration of polymer in solution
would gradually change as the encapsulant dissolves. Gun and Routh (2013) explored
the possibility of using poly(lactic-co-glycolic) acid to encapsulate a gelling agent that
is gradually released. In the context of drug delivery in the body, Makadia and Siegel
(2011) detail how the release profile of this encapsulant can be tuned by varying the
composition. Alternatively, the polymer may be thermally activated, gelling once it
passes through an activation temperature. Tran-viet et al. (2014) explored the ther-
mal response of poly(NIPAM) to this end. The viscosity may also vary on account
of reaction. Polyacrylamide (PAA) solutions are commonly used in oil field displace-
ments (Kamal et al., 2015). In addition to gradual release from encapsulant, time
dependency can be introduced to PAA solution viscosity through tailoring reaction
kinetics. Han et al. (1995) demonstrated that the addition of chromium ions to PAA
led to the formation of a cross-linked gel structure, and that the addition of acetic
acid to this solution delayed the rate of gel formation. Furthermore, we hypothesise
that by limiting the initial concentration of chromium ions, the rheological properties
of the final polymer solution, including the cross-links, can also be controlled. In the
case of PAA polymerisation, Lee et al. (2012) measured the gradual change in solution
viscosity over the course of the reaction and demonstrated that the rate of reaction
is dependent on the concentration of various reagents. With polymer solutions (in-
cluding those of PAA), shear can cause the polymer chains to align which can lead
to a decrease in the effective viscosity. Thus, polymer solutions with time-dependent
viscosity may exhibit shear-thinning behaviour.
3.2 Introduction 43
Drawing on the results of previous work related to time-dependent changes in
viscosity, we now explore how such variations of viscosity may impact the stability
of a moving interface, and given the evolving rheology, we also account for possible
variations in injection rate. Such effects may be key for the class of problem in which
a finite volume of fluid is to be injected in a finite time, and in this case there is the
intriguing possibility that by varying the flow rate with time the final amplitude of
instability can be further reduced. In exploring this class of problem, we are guided
by the pioneering works of Dias et al. (2010) and Dias et al. (2012) who examined the
impact of changes in flow rate on the growth of viscous fingers. In particular, Dias
et al. (2012) considered the problem of injection of a finite volume of fluid in a finite
time, and using variational calculus, developed an expression for the injection flow
rate as a function of time which minimises the final amplitude of instability.
We have arranged the chapter as follows. In §3.3 we consider the impact of time-
dependent rheology in a uni-directional displacement, in which there is a continuous
range of unstable wavenumbers. we assess the relationship between the time-dependent
viscosity of the injected fluid and the optimal injection rate. In §3.4 we extend the
results to consider injection from a central source in which the instability leads to
a discrete series of modes with increasing azimuthal wavenumber, n, of the form
An(t) exp(inθ) (cf. Paterson (1981)). In this case, the analysis is more complex
owing to the continual stretching of the interface which progressively suppresses the
instability of lowest modes as the interface grows, while surface tension suppresses
the instability of the highest modes (cf. Cardoso and Woods (1995)). In order to
make progress, following Dias et al. (2012) we approximate the discrete spectrum of
modes with a continuous spectrum, noting that this will lead to an upper bound on
the actual growth rate. We find that for a given initial radius R0, there is a particular
total injection time, t†f , for which the system is just stable, provided the injection rate
decreases with time according to Q = a(b + ct)−1/3, where the constants a, b, and c
are found in our analysis. If tf decreases to values tf < t
†
f , for this given R0, then
the optimal injection rate gradually evolves towards the simple linearly increasing rate
44 Viscous instability with time-dependent rheology
Q = d+et, where d and e again are found from the analysis. This limit coincides with
the results of Dias et al. (2012), who investigated the optimal injection rate in the limit
tf << t
†
f . We then illustrate how a gradual increase in the viscosity of the injected
fluid, for example resulting from slow activation of gel, modifies this optimal solution
so that, as for the uni-directional flow, the optimal injection rate increases with time
from a smaller to a larger value than for the case in which the injected fluid has a
constant viscosity equal to the initial viscosity. Given that many polymer solutions
exhibit shear thinning rheology, in §3.5 we generalise the analysis to account for such
rheology in the calculation of the optimal flow rate with time. To this end, we draw on
the analysis of Wilson (1990) and Mora and Manna (2009) who developed dispersion
relations for the growth of viscous fingers in non-Newtonian fluids migrating through
a Hele-Shaw cell. We compare our results to Fontana et al. (2014), who generalised
the analysis of Dias et al. (2012) to air displacing a shear-thinning fluid. It was found
that, with radial flow, the optimal injection rate involved a more rapid injection at
early time and slower injection later as compared to the Newtonian case.
3.3 Newtonian Flow in a Rectilinear Geometry
3.3.1 Formulation
The depth-averaged flow of fluid in a Hele-Shaw cell is governed by the following
equations (Saffman and Taylor, 1958):
∇ · u = 0 and u = − b
2
12µ
∇p, (3.1)
where u denotes velocity, b is the plate separation, µ the viscosity of the fluid and p
is the pressure.
By applying the interfacial boundary conditions of continuity of velocity and the
jump in pressure on account of interfacial tension, and exploring the stability of the
3.3 Newtonian Flow in a Rectilinear Geometry 45
interface to small sinusoidal perturbations along the interface, Saffman and Taylor
(1958) obtained the dispersion relation
A˙
A
= σk =
µ2 − µ1
µ2 + µ1
Uk − b
2T
12(µ1 + µ2)
k3, (3.2)
where T is the interfacial tension, A(k, t) is the amplitude of a perturbation of wavenum-
ber k, U(t) is the velocity of the interface and subscript 1 denotes the displacing fluid
whereas subscript 2 is the displaced fluid. With injection of a finite volume of fluid in
a finite time tf , such that the interface migrates a distance xf , this relation may be
expressed in dimensionless form
σˆk =
(1− V )Uˆ kˆ − τ kˆ3
1 + V
, (3.3)
where the hat notation denotes a dimensionless variable, kˆ = kxf , tˆ = t/tf , the
stability parameter τ =
b2Ttf
12µ2x3f
, Uˆ(tˆ) = Utf/xf , σˆ = σtf and V (tˆ) =
µ1(tˆ)
µ2
. We
now drop the hat notation for convenience and henceforth work with dimensionless
variables. The wavenumber, kmax, with maximum growth rate, σmax is given by
kmax =
√
(1− V )U
3τ
, (3.4)
where
σmax =
2
3
√
3τ
(1− V ) 32
(1 + V )
U
3
2 . (3.5)
At each time, the growth rate of any mode, σ(k, t) < σmax(t) and so an upper
bound on the natural logarithm of the final amplitude of any mode is given by
I = log(Af ) =
∫ 1
0
σmaxdt (3.6)
where A is scaled with the starting amplitude of any mode. In order to find the
injection rate, U(t), which minimises I subject to the requirement
∫ 1
0 Udt = 1, we can
follow the Euler-Lagrange framework of variational calculus and seek a solution for
46 Viscous instability with time-dependent rheology
U(t) of the equation
d
dt
(
∂σmax
∂U
)
= 0. (3.7)
This leads to the ordinary differential equation for U(t)
(5 + V (t)) V˙ (t)U(t)−
(
1− V (t)2
)
U˙(t) = 0, (3.8)
with solution
U(t) = U∗(t) = Ω
(1 + V (t))2
(1− V (t))3 ,
where Ω =
(∫ 1
0
(1 + V (t))2
(1− V (t))3dt
)−1 (3.9)
and the star (∗) indicates the variable is of the optimal value. If the viscosity ratio
is constant, V˙ = 0, then (3.8) predicts that the constant flow, U∗ = 1 is optimal. It
follows that
k(t) = k(t)∗max =
√
Ω
3τ
(1 + V (t))
(1− V (t)) (3.10)
and
σmax(t) = σ(t)
∗
max =
2
3
√
Ω3
3τ
(1 + V (t))2
(1− V (t))3 . (3.11)
A key parameter in the present model is τ . Analysis of the upper bound on the
amplitude of perturbations, (3.6), shows that this may be re-expressed in the form
Af = exp
(∫ 1
0
f(U(t), V (t))dt/τ 1/2
)
, (3.12)
where f(U, V ) takes on different functional forms depending on whether the injection
rate is constant or follows the optimal injection rate. It follows that A
√
τ ,
√
τσmax and
√
τkmax are independent of τ .
3.3.2 Effect of a Gradual Increase in Viscosity
As mentioned in §3.2, the time-dependent viscosity could be tailored through choosing
the thickness or composition of polymer micro-encapsulant, altering reagent concen-
3.3 Newtonian Flow in a Rectilinear Geometry 47
trations to control polymerisation or cross-linking reaction kinetics or altering the
temperature to trigger a thermally activated viscosity change. In the case of release
from micro-encapsulant, we assume the viscosity of the displacing fluid is directly pro-
portional to the concentration of released polymer and model this concentration as
following the release profiles shown in Makadia and Siegel (2011). Thus, the viscosity
ratio has the form
V (t) = Vo + (Vf − Vo)
(
1− e−θt
)
, (3.13)
where Vo is the initial viscosity ratio, Vf is the long-time asymptotic viscosity ratio
and θ is the ratio of the injection time, tf , to the characteristic time for the viscosity
change or ’gelling’. In figures 3.1(a) and 3.1(b) we illustrate the optimal solutions for
the cases in which θ = 1, corresponding to slow gelling, and θ = 10, corresponding
to faster gelling. In both cases, the change in viscosity of the injected fluid is also
defined by the ratio Vf/V0 = 90. In the case θ = 10, the viscosity of the injected fluid
reaches this target value early in the flow, whereas for the slow-gelling case, θ = 1,
the viscosity is increasing for the duration of the injection. For the faster-gelling case,
θ = 10, the optimal flow rate strategy involves injection of the majority of the fluid
once the viscosity has reached its maximum value since adverse viscosity differences
are reduced at these later times. In contrast, with slow gelling, the change in viscosity
is smaller and so the optimal flow rate increases gradually with time. In figures 3.1(c)
and 3.1(d) we show how the most unstable wavenumber
√
τkmax evolves in time, for the
optimal flow rate (dashed lines), and, for reference, for the case of a constant injection
rate (solid lines). With a constant flow rate, the most unstable wavenumber decreases
with time, owing to the increasing stability of the system. In contrast, the optimal
injection profile actually leads to a gradual increase in the most unstable wavenumber
with time as a result of the increasing flow rate, even though the viscosity ratio falls
with time. As a result of the different evolution of the maximum growth rate with
time, we find that for the optimal flow, the maximum growth rate actually increases
48 Viscous instability with time-dependent rheology
with time owing to the progressively faster injection, whereas for a constant injection,
the growth rate falls with time (figures 3.1(e), 3.1(f)).
We now explore the final amplitude of the perturbation, once the finite volume
of fluid has been injected. In order to compare the benefit of the optimal injection
rate with the case of constant injection, for different values of θ, it is convenient
to investigate the variation of A
√
τ
f with θ for each case. For clarity, figure 3.2(a)
illustrates the difference between these two amplitudes as a function of wavenumber
in the case Vf = 0.9 and Vo = 0.01. Figure 3.2(b) illustrates the ratio of the two
final largest amplitudes as a function of θ. Both figures illustrate that the maximum
benefit arises when θ is close to unity so that the viscosity is changing over the whole
period of injection. For small or large θ the difference is much smaller, since in either
of these limits, the majority of the injection occurs with either the original or the final
viscosity. In figure 3.2(c), curves are given for V0 = 0.01 and Vf = 0.5, 0.7 and 0.9
corresponding to polymer solutions whose final viscosity is progressively larger. Each
curve has a similar shape, but the magnitude of the reduction in amplitude associated
with using the optimal injection strategy is greater when the change in viscosity of the
polymer gel is greater. This is further illustrated by figure 3.2(d), which shows how
the magnitude of reduction in amplitude diminishes as V0 approaches Vf for θ = 2.
Finally, in order to illustrate the effect of decreasing the overall injection time τ , in
figure 3.2(e) we show contours of the ratio of the final amplitude associated with
constant and optimal injection in θ− τ space, for the case V0 = 0.01 and Vf = 0.9. As
the value of τ is reduced from unity, corresponding to a shorter injection time, it is
seen that the suppression of final amplitude is enhanced. For example, for τ = 10−3,
log
(
Af
A∗
f
)
≈ 2.5, which corresponds to amplitude reduction by a factor of 12. We note
that there is an upper limit on the validity of this analysis; for smaller values of τ and
therefore larger amplitudes, nonlinear effects dominate behaviour and the amplitudes
will be different to those predicted by linear stability analysis.
3.3 Newtonian Flow in a Rectilinear Geometry 49
 
 
t
UV
0 0.5 1
0
1
2
3
0
0.2
0.4
0.6
0.8
1
(a)
t
UV
0 0.5 1
0
1
2
3
0
0.2
0.4
0.6
0.8
1
(b)
 
 
√
τkmax
√
τk∗
max
t0 0.5 1
0
0.2
0.4
0.6
(c)
t
0 0.5 1
0
0.2
0.4
0.6
(d)
 
 
t
√
τσmax
√
τσ∗max
0 0.5 1
10−6
10−3
100
(e)
t
0 0.5 1
10−6
10−3
100
(f)
Fig. 3.1 Left column: θ = 1 (slow gelling). Right column: θ = 10 (rapid gelling).
Top row: viscosity ratio profiles with parameters θ, V0 = 0.01 and Vf = 0.9 (primary
axis) and the optimal velocity profiles (secondary axis). Middle row: the evolution
of
√
τkmax for the optimal (dashed lines) and constant (solid lines) flow rate profiles.
Bottom row: the evolution of
√
τσ.
50 Viscous instability with time-dependent rheology
k
A
√
τ
k,f
θ = 0
θ = 1
0 10 20 30
1
1.25
1.5
(a)
 
 
A
√
τ∗
f
A
√
τ
f
θ
10−2 100 102
1
1.25
1.5
(b)
 
 
0.5
0.7
0.9
θ
(
Af
A∗
f
)√τ
10−2 100 102
1
1.02
1.04
1.06
1.08
(c)
 
 
0.5
0.7
0.9
V0
(
Af
A∗
f
)√τ
0 0.5 0.9
1
1.02
1.04
1.06
1.08
(d)
 
 
τ
log
(
Af
A∗
f
)
θ
10−4 10−3 10−2 10−1
1
2
3
4
5
10−1
100
101
102
(e)
Fig. 3.2 (a) The final amplitude of perturbations of wavenumber k for the viscosity
invariant case (θ = 0, grey) and a gelling case (θ = 1) for a constant flow rate (solid
line) and the optimal flow rate (broken line). (b) The final maximum amplitudes
varying as a function of θ for Vf = 0.9. (c) The ratio of the these two final amplitudes
as a function of θ for a series of Vf . (d) The effect of increasing V0 for a series of Vf
and θ = 2. (e) A contour plot of the final amplitude ratio varying with τ and θ for
V0 = 0.01 and Vf = 0.9.
3.4 Injection from a point source 51
3.4 Injection from a point source
We now turn to the more complex problem of injection from a point source rather
than a line source. Paterson (1981) derived a dispersion relation for the growth of
perturbations on a circular interface of radius R in terms of a series of discrete az-
imuthal modes n, which depends on the viscosity ratio across the interface, V , the
surface tension T and the permeability K, taken to be b2/12 for a Hele-Shaw cell. For
injection of a fixed volume of fluid, πb(R2f−R20) over a time tf , we can scale the growth
rate with 1/tf and the radius with Rf , leading to the dimensionless growth rate
σn(R, R˙) =
R˙
R
(
1− V
1 + V
n− 1
)
− τ
R3(1 + V )
n(n2 − 1) (3.14)
where the stability parameter τ =
b2Ttf
12µ2R3f
. Although n corresponds to a series of
discrete modes, we can find an upper bound on the maximum growth rate at each
time by treating n as a continuous variable (Dias et al. (2012)) leading to
nmax(R, R˙) =
√
1
3
(
1 + 2ΛR˙R2
)
, (3.15)
σmax =
(
1− V
1 + V
)
1
3
√
3ΛR3
(
1 + 2ΛR˙R2
)3/2 − R˙
R
, (3.16)
where Λ = (1−V )
2τ
. Since the amplitude of each mode grows as the exponential of the
integral of the growth rate of that mode, the ultimate amplitude of the instability will
be smaller than the expression
A = exp
(∫ 1
0
σmax(R˙, R, t)dt
)
, (3.17)
If we find a minimum value for A by varying dR/dt through all possible functions
which satisfy the boundary conditions, then this will provide an upper bound on the
amplitude of the instability. To this end we apply the Euler-Lagrange equation
52 Viscous instability with time-dependent rheology
d
dt
(
∂σmax
∂R˙
)
=
∂σmax
∂R
, (3.18)
leading to the following differential equation governing the optimal injection rate:
1 + ΛR2R˙ + Λ2R5R¨− V˙
2τ(1 + V )
(
2R3 + Λ(5 + V )R˙R5
)
= 0. (3.19)
3.4.1 The Constant Viscosity Regime
In the case that the injection fluid has constant viscosity, the evolution of R is given
by the equation
1 + ΛR2R˙ + Λ2R5R¨ = 0, (3.20)
subject to the constraint that R = R0 < 1 at t = 0 and R = 1 at t = 1, as noted by
Dias et al. (2012). Subject to the appropriate boundary conditions, this equation has
an exact solution which is available only in the special case that Λ = Λ†. The solution
is
R =
(
t+ to
1 + to
)1/3
, (3.21)
where to = R
3
0/(1 − R30). The special value Λ† can be found from substituting (3.21)
into (3.20), giving
1 + (1/3)Λ†(1−R30) = (2/9)Λ†2(1−R30)2, (3.22)
which requires
Λ†(R0) = 3
1
(1−R30)
. (3.23)
This solution can be substituted into (3.15) to show that it fixes the most unstable
mode to the mode nmax = 1 for the injection duration. The solution can also be
substituted into (3.16) to show that it corresponds to a maximum growth rate
3.4 Injection from a point source 53
σ†max = −
2
3(t+ t0)
V
(1 + V )
. (3.24)
This solution was recently found independently using a different approach (Batista
et al., 2016). If V = 0 this solution is neutrally stable for the duration of flow (σ†max = 0)
and the result quantitatively matches the flow rate that would be found by setting
σ = 0 in the dispersion relation (3.14), as described in the previous chapter (see §2.3.1).
If 0 < V < 1 then σ†max < 0 ∀ t and the solution (3.23) leads to a decay of each mode.
We have solved (3.20) numerically for Λ = 3Λ† and Λ = 10Λ†, with fixed R0, and
these solutions are shown in figure 3.3 (solid lines). It is seen that as Λ increases to
values much larger than Λ†(R0), the variation of the optimal injection rate, Q∗ = 2RR˙,
with time changes in character from the slowly decaying flow rate for values close to
Λ† (included as a dashed line) to a linearly increasing flow rate in the case Λ >> Λ†
(dashed-dotted line). This occurs since if the length of the injection duration permits,
Λ & Λ† then the destabilising effects can be balanced against interfacial tension which
is achieved by a decreasing flow rate. On the other hand, if the injection duration is
much shorter, Λ >> Λ†, then it is optimal to minimise the amount of time for which
any particular mode is the most unstable by transitioning through the modes quickly,
which is seen to be achieved by a linearly-increasing injection flow rate. The solution
for Λ >> Λ† corresponds to the solution proposed by Dias et al. (2012),
R = R0 + (1−R0)t (3.25)
and the numerical solutions (solid lines) smoothly connect this limit with the solution
(3.21) for Λ = Λ†(R0) which coincides with the analytical solutions in the formulation
of Batista et al. (2016). We only consider values Λ ≥ Λ† since for values less than
this the injection duration is longer than that required to fully stabilise the flow with
interfacial tension. This leads to the solution of (3.20) involving a period of injection
followed by a period of suction. Therefore, Λ = Λ† marks the limit of values for which
this optimisation approach is valid.
54 Viscous instability with time-dependent rheology
 
 
100Λ†
Λ†
Q
∗
t
0 0.2 0.4 0.6 0.8 1
0
1
2
3
Fig. 3.3 The flow rate, Q∗(t), that minimises
∫ 1
0 σmaxdt for Λ = Λ
† (dashed line), Λ =
3Λ† and Λ = 10Λ† (solid lines) and Λ = 100Λ† (dashed-dotted line) while R0 = 0.1.
We now explore the evolution of the instabilities over time, to assess the value of
the predicted optimal injection rate, which is based on minimising the upper bound
of the growth rate as a function of time, and we compare the results with the case of
a simple constant injection rate.
It is useful to recall that the upper bound on the growth rate – σmax – has been
estimated by assuming a continuous distribution of modes, whereas in practice, the
modes are quantised in the azimuthal direction. We anticipate that for cases in which
the most unstable azimuthal wavenumber is large, the model will offer a better bound
than for slower injection rates when only the lowest modes are unstable.
To proceed, we first explore the evolution of the upper bound as a function of time
during the injection, from the case Λ ∼ Λ† to the case Λ≫ Λ†. The only parameters
which govern the optimal flow rate are Λ and R0, however, σmax additionally depends
on the viscosity ratio V . In figure 3.4(a), the bound is shown with broken lines for two
cases in which V = 0 and R0 = 0.2: Λ = Λ
† (black) and Λ = 3Λ† (blue) as labelled.
The solid lines correspond to four cases in which V = 0.1 and R0 = 0.2: Λ = Λ
†
(black), Λ = 3Λ† (blue), Λ = 10Λ† (green) and Λ = 30Λ† (red).
It is seen that in the case V = 0, the upper bound on growth rate is initially zero
and with Λ = Λ† this remains the case throughout the injection duration since it is
the neutrally stable solution. With Λ = 3Λ†, the instability gradually develops with
3.4 Injection from a point source 55
 
 
30
10
3
1
σ
∗ m
a
x
t
0 0.2 0.4 0.6 0.8 1
−1
−0.5
0
0.5
1
(a)
 
 
0.5
0.3
0.1
S
U
R
0
Λs
Λ†
, Λu
Λ†
100 101 102 103
0.2
0.3
0.4
0.5
0.6
0.7
(b)
Fig. 3.4 (a) The evolution of the upper bound of growth rate, σ∗max, for V = 0 (broken
lines), V = 0.1 (solid lines), R0 = 0.2 and a series of values of
Λ
Λ†
that increases from
the bottom of the figure upwards. (b) A regime diagram for the conditions σ∗max < 0 ∀
t (broken lines; region S) and σ∗max > 0 ∀ t (solid lines, region U) for a series of values
of V increasing from left to right (as labelled).
time. In the case V = 0.1, when Λ = Λ† the bound is in fact stable throughout the
injection period, whereas with Λ = 3Λ†, the bound is initially stable, but eventually
becomes unstable just before all the fluid has been injected. As the value of Λ gradually
increases, the bound becomes unstable at progressively earlier times until eventually
it is unstable as soon as the injection commences. If we extrapolate from these results,
it follows that for each initial radius R0 there is a critical value of Λ, above which the
bound is always unstable (Λu, solid lines), and a second critical value of Λ below which
the bound is always stable (Λs, dashed lines). These values are shown in figure 3.4(b).
Curves are given for viscosity ratios of 0.1 (blue), 0.3 (red) and 0.5 (black). For a given
radius, the critical value of Λ below which the modes are always stable increases with
the viscosity ratio owing to the reduced viscous destabilisation of the front. Similarly,
the critical value of Λ for which the modes are always unstable increases with the
viscosity ratio.
The number of discrete modes that become unstable during the course of the
injection can be calculated from the evolution of the continuous bound nmax (3.15).
This is shown for the optimal injection case in figure 3.5(a) where R0 = 0.1. For
Λ = Λ† (black line), the most unstable mode is fixed to nmax = 1 for the duration of
56 Viscous instability with time-dependent rheology
the injection and the line coincides with the x-axis, as discussed. For Λ = 3Λ† (blue
line), initially mode 1 is the most unstable, but at the end of the injection mode 2
is the most unstable, whereas for Λ = 10Λ†, once again initially mode 1 is the most
unstable, but at the end of the injection mode 4 has the largest growth rate. For
Λ = 100Λ† (green line) mode 2 is the most unstable at early times, while the most
unstable mode at the end of the injection is mode 13 (not shown). Similarly, for
Λ = 1000Λ† (purple line), initially mode 4 is the most unstable, whereas the most
unstable mode at the end of injection has now increased to mode 42. As the value
of Λ increases, the discrete azimuthal mode that is initially the most unstable has a
larger wavenumber and furthermore the span of different modes that become the most
unstable during injection increases too. We now investigate the evolution of these
discrete modes (n = 2, 3, ...).
The coloured lines in figures 3.5(b-d) show the evolution of the growth rates of
the discrete modes n = 2, 3, ... during the injection period for the case of injection
at a constant rate. The upper bound σmax is included as the thick black line. The
three figures (b-d) correspond to the cases Λ = 3Λ† (panel b), Λ = 10Λ† (panel c)
and Λ = 100Λ† (panel d). For comparison, in figure 3.5(e-g) we show the evolution
of the growth rates of the discrete modes for the case in which the injection follows
the optimal injection rate, the solution of (3.20). Panels (e-g) correspond to the same
values of Λ as in panels (d-e). In the final three panels, (h-j), we present the amplitude
of each of the discrete modes as a function of time for the case in which the injection
follows the optimal injection rate.
When Λ = 3Λ†, mode 2 is the only mode to become unstable for both the constant
(b) and optimal (e) injection flow rate profiles. In the case of constant injection,
the bound on growth rate is initially large, then decreases to zero, and subsequently
rebounds to a local maximum before decaying away. The bound is initially large on
account of a non-physical mode, n < 1, that causes the interfacial tension term to
become positive in the dispersion relation (3.14). The growth rate of mode 2 becomes
positive during the flow but does not become as large as the bound, even at the end
3.4 Injection from a point source 57
of the injection. In the case of optimal injection, the bound monotonically increases
from near zero to a final value larger than that of the constant injection case. Once
again mode 2 is the only mode to become unstable during the flow, but in the optimal
injection case the magnitude of the growth rate is equal to the bound σmax at the end
of the injection. The amplitude of mode 2 at the end of the flow is small (of the order
10−1) because the late-stage growth of the mode is insufficient to overcome the early
time decay whilst the mode is stabilised by interfacial tension.
Where Λ = 10Λ†, the bound σmax in the constant injection case (c) initially in-
creases, but after reaching a maximum it gradually decreases with time. The maxi-
mum in growth rate is larger and earlier in time than for Λ = 3Λ†. Mode 2 becomes
unstable at an early time, with modes 3 and 4 following over the course of injection.
The growth rate of each of the modes reaches a maximum value then decays away,
leading to a cascade to higher modes. In the case of optimal injection (f), the bound
on growth rate is once again initially small and then increases monotonically. The
onset of instability of mode 2 is delayed relative to the constant injection case and
features a smaller maximum than the constant injection case. The onset of instability
of modes 3 and 4 is also delayed, but the maxima are larger in growth rate than the
constant injection case. In contrast to the constant injection case, mode 5 becomes
unstable during the flow. Figure 3.5(i) illustrates that the amplitude of mode 2 at the
end of the injection phase is greater than the initial value, owing to the dominance
of the instability at the later stages of injection, even though it is initially stable.
However, the amplitudes of modes 3-5 do decay to smaller values.
For the case Λ = 100Λ†, in the constant injection case (d), the bound on growth
rate is initially very large, then monotonically decreases for the duration of the flow
and closely follows the locus of the mode that is the most unstable at that time. Mode
2 (bold orange line) is initially unstable, and also monotonically decreases in growth
rate. Modes 6-16 behave like the lower modes when Λ = 10Λ†, insofar as they become
unstable during the flow, reach a maximum in growth rate and then subsequently decay.
Mode 11 is highlighted with a bold purple line. Curiously, the maximum in growth
58 Viscous instability with time-dependent rheology
rate of a discrete mode can occur before before it becomes tangential to the bounding
curve σmax. When the fluid is injected with the optimal injection rate, (g), the bound
on growth rate is initially non-zero but nonetheless considerably smaller than the
constant injection case. Mode 2 is initially unstable. The growth rate increases to a
maximum coincident with the bound σmax, before decreasing as the radius increases
further. Modes 3-22 all become unstable during the injection. The behaviour of each
mode follows the same pattern of increasing towards the maximum growth rate and
then decaying as the radius increases to larger values. However, owing to the fact that
many of the modes are initially stable, the modes with larger azimuthal wavenumber,
which spend a smaller fraction of the injection period being unstable than stable, finish
with an amplitude smaller than unity. Only modes 2-4 feature overall growth during
the injection.
As the value of Λ increases, the number of different modes that become unstable
increases and the bound – σmax – becomes a better approximation of the locus of the
most unstable mode for both the constant and optimal injection cases. There are
two factors that appear to contribute to stabilisation: i) the large initial growth rates
featured in the case of constant injection are mitigated by slower flow at early times
in the optimal case and ii) a larger number of modes become unstable in the optimal
case, giving each mode relatively less time to grow.
Figure 3.6 is a comparison of the final amplitude as computed from integrating the
bound – σ∗max (grey lines) – and the largest final amplitude of any discrete mode (black
lines) as the total time of injection is reduced, i.e. as the value of Λ is increased for
parameters V = 0 and R0 = 0.1. This figure also compares the amplitudes resulting
from constant injection (solid lines) and optimal injection (broken lines). The discrete-
mode curves are not smooth. The overestimation of the final amplitude as computed
by integrating the bound of all modes as compared to the actual growth rate of each
individual mode can be seen by comparing the grey and black lines. The overall
stabilisation gained by optimal injection can be seen by comparing the broken and
solid lines. The effect of stabilisation is more pronounced when considering each mode
3.4 Injection from a point source 59
 
 
1000
100
10
3
1
t
0 0.2 0.4 0.6 0.8 1
100
101
102
(a) log10(nmax(t))
t
σ
0 1
0
0.3
(b)
t
σ
0 1
0
2.5
(c)
t
σ
0 1
0
8
(d)
t
σ
∗
0 1
0
0.3
(e)
t
σ
∗
0 1
0
2.5
(f)
t
σ
∗
0 1
0
8
(g)
t
lo
g 1
0
(A
∗ )
0 1
10−2
101
(h)
t
lo
g 1
0
(A
∗ )
0 1
10−4
102
(i)
t
lo
g 1
0
(A
∗ )
0 1
10−5
103
(j)
Fig. 3.5 (a) The evolution of the most unstable mode nmax for R0 = 0.1 and a series
of Λ
Λ†
that increases from the bottom of the page upward. (b-g) The growth rate of
discrete modes n for a constant injection rate (2nd row) and the optimal injection
rate (3rd row) for V = 0, R0 = 0.1 and Λ = 3Λ
† (left column), Λ = 10Λ† (centre
column) and Λ = 100Λ† (right column). The mode number increases from the LHS
to the RHS of the figures. (h-j) The amplitude of perturbations following the optimal
injection rate. Initially, the lowest modes have the largest amplitude. Mode 2 (11) is
highlighted with thick line in orange (purple). The growth-rate bound derived from
treating n as a continuous variable is shown as the thick black line.
60 Viscous instability with time-dependent rheology
Λ
Λ†
A
f
101 102 103
100
102
104
106
Fig. 3.6 The final amplitude varying with Λ
Λ†
for V = 0, R0 = 0.1 for constant injection
(solid lines) and optimal injection (broken lines) computed with discrete modes (black)
and a continuous series of modes (grey).
individually, and – as the curves are seen to diverge – the benefit of optimal injection
becomes larger as the total injection time is reduced.
3.4.2 Effect of a Gradual Increase in Viscosity
As for the uni-directional flow problem, when the viscosity of the injected fluid grad-
ually increases with time, we expect that the optimal flow solution will involve an
increase in the flow rate with time relative to the case of a constant viscosity. The
parameter Λ = (1−V )
2τ
is no-longer a constant but varies with the change in viscosity
between Λ0 =
(1−V0)
2τ
and Λf =
(1−Vf )
2τ
. If the duration of the flow tf is small, or
interfacial tension T relatively weak given the duration of the flow, then, assuming
1 − V = O(1), we can write τ → 0. In this case, expanding (3.19) in powers of τ
whilst assuming all other variables remain O(1), leads to
(5 + V )V˙ R˙− (1− V 2)R¨ = 0, (3.26)
which has the solution
R˙ = (1−R0)Ω(1 + V (t))
2
(1− V (t))3 . (3.27)
3.4 Injection from a point source 61
 
 
θ = 2
θ = 1
θ = 0
U
R
0
Λ0,u
Λ†
0
101 102 103
0.2
0.3
0.4
0.5
0.6
0.7
Fig. 3.7 The critical value of Λ0, above which the system is always unstable for the
gelling case in which V0 = 0.3 and Vf = 0.5. A series of values of gelling rate, θ, is
shown.
This closely corresponds to the rectilinear case (3.9). To explore this behaviour, we use
the same example of a gelling process as previously described using equation (3.13).
In an analogous fashion to figure 3.4(b), for a given initial radius there is a particular
value Λ0,u such that if Λ0 > Λ0,u then the bound σ
∗
max is unstable for the duration of
the flow. In figure 3.7, we present this critical value, Λ0,u, for V0 = 0.3 and Vf = 0.5.
The critical value of Λ0 for the non-gelling case, θ = 0 (red), can be recognised from
figure 3.4(b). In the gelling case, θ = 1 (blue), the critical value Λ0,u is larger than for
the non-gelling case since the growth rate can be negative at any stage of the flow. In
the faster-gelling case, θ = 2, the critical value Λ0,u is larger still.
Figure 3.8 shows how the optimal flow rate Q∗(t) varies for a series of increasing
gelling rates, θ, for parameters V0 = 0, Vf = 0.5, Λf = 1250 and R0 = 0.35. The left-
hand panel compares the optimal flow rates when the viscosity change is slow relative
to the injection duration, θ . 1, whereas the right-hand panel corresponds to faster
gelling, θ & 1. For the constant viscosity case, θ = 0, the optimal injection profile
analogous to the result of Dias et al. (2012), equation (3.25), is shown in figures 3.8(a)
and (b) as broken lines. For a very slow gelling rate, θ = 0.1 (red line), the optimal
flow rate is initially smaller than the non-gelling case but then involves injecting with a
larger injection rate than the non-gelling case at later times. This may be understood
in terms of the system optimising the benefits of the higher viscosity at later times in
62 Viscous instability with time-dependent rheology
 
 
1
0.1
Q
∗
t
0 0.2 0.4 0.6 0.8 1
0
1
2
3
(a) θ . 1
 
 
30
10
1
Q
∗
t
0 0.2 0.4 0.6 0.8 1
0
1
2
3
(b) θ & 1
Fig. 3.8 The optimal flow rate, Q∗(t), for a series of different gelling rates θ for R0 =
0.35, V0 = 0, Λf = 1250 and Vf = 0.9 found using (3.27) (solid lines) and the invariant
viscosity case (broken line). In (a), the timescale of viscosity change compared to the
flow is small, whereas in (b) it is large. In (a) the numerical solution to (3.19) is shown
(crosses).
the injection process. As θ increases to the value for which the time scale of gelling and
injection, tf , are matched, θ = 1 (blue line), these features become more pronounced.
We also note that full numerical solution of the ODE (3.19) (crosses in figure 3.8(a))
and the asymptotic solution (3.27) (solid line) agree in this limit of a small total
injection time, or large Λf . Turning to the right-hand panel, for θ = 10 (red line), the
optimal flow rate is more akin to the non-gelling case at late times, however, initially
the flow rate is approximately zero. For θ = 30 (black line; right-hand panel), the
optimal flow rate is again initially approximately zero, and subsequently closely follows
the non-gelling case since the viscosity ratio is effectively the final value, Vf , for the
duration of the flow.
Figure 3.9 explores the benefit of injecting with the optimal flow rates shown in
figure 3.8 at the end of the injection by illustrating the amplitudes of the discrete
modes, n = 2, 3.... In the absence of gelling, θ = 0 (grey lines), the optimal injection
solution (broken line) shows a smaller peak in amplitude of any of these modes than
the constant injection case (solid line); this illustrates the benefit of optimal injection
as discussed in §3.4.1. For θ = 1 (black lines), the gelling leads to a reduction in the
final maximum amplitude with constant injection since the adverse viscosity gradients
3.5 Injection of a Shear-Thinning Fluid 63
are reduced, and deploying the optimal injection strategy leads to further stabilisation.
The variation of the magnitude of these maxima as the rate of gelling, θ, increases
is shown in figure 3.9(b). As the rate of gelling increases, the amplitudes decrease
to a smaller plateau. The amplitude of each mode at the end of the injection period
when deploying the optimal injection rate solution is consistently smaller than those
associated with using a constant injection rate and the ratio of these two values, which
may be interpreted as a measure of the benefit of deploying the optimal injection
strategy, is plotted in figure 3.9(c) (blue line). As with unidirectional flow, it can be
seen that the benefit of using the optimal flow rate is largest when the time scale of
gelling is similar to that of the flow, θ ≈ 1. Also, as the overall injection rate increases,
the benefit of deploying the optimal injection rate is greater still, as illustrated for
example with the case Λf = 2500 (black line).
3.5 Injection of a Shear-Thinning Fluid
In this section, we generalise the results of §3.3 to describe the injection of a fluid
that exhibits shear-thinning rheology. We adopt a power-law model to describe shear-
thinning:
φ = k1γ˙
m, (3.28)
where φ is the shear stress, γ˙ is the strain rate and k1 and m are empirical fitting
parameters known as the consistency index and flow behaviour index respectively. As
pointed out by Mora and Manna (2009), the unperturbed interface velocity features
zero shear at the centreline in Poiseuille flow. However, the power-law model we
have adopted (3.28) neglects the Newtonian plateau in viscosity that might be seen
at low shear rates. We shall discuss why this is necessary in the next paragraph.
The Newtonian plateau is described in more detail by Martyushev et al. (2015). By
ignoring it, they show how the power-law model modifies the Hele-Shaw flow equation
(3.1) to become
64 Viscous instability with time-dependent rheology
A
f
,n
n
0 10 20 30 40
100
102
104
106
(a)
 
 
A∗
A
A
f
θ
10−2 100 102
100
102
104
106
(b)
 
 
1250
2500
A
′
A
θ
10−2 100 102
100
101
102
103
(c) Λf series
Fig. 3.9 (a) The final amplitude of each mode n for the viscosity invariant case (θ = 0,
grey lines), for a constant flow rate (solid lines) and the equivalent optimal flow rate
(broken lines) where V0 = 0, Vf = 0.5, Λf = 1250 and R0 = 0.35. A gelling example
is also shown (θ = 1, black). (b) The variation of both maximum amplitudes with θ.
(c) The ratio of the constant to optimal injection final amplitude plotted against θ for
two values of Λf .
3.5 Injection of a Shear-Thinning Fluid 65
u = − m
2m+ 1
( |∇p|
k1
)1/m (
b
2
)(m+1)/m
∇p
|∇p| . (3.29)
Ghannam and Esmail (1998) have shown that this provides a reasonable representation
of the behaviour of Polyacrylamide (PAA) over a large range of shear rates. If it were
desired to incorporate the low-shear Newtonian plateau, a model such as the Carreau
model (Carreau et al., 1979) could be used. There may additionally exist an infinite-
shear Newtonian plateau, in which case the Cross model (Cross, 1965) could be used.
However, following the analysis of Mora and Manna (2009) for a Carreau fluid, we
note that it is not possible to explicitly write the pressure gradient as a function of
the mean velocity, which leads to it not being possible to write the dispersion relation
in this way. A fully numerical treatment would therefore be required to proceed with
the optimisation and account for the low-shear Newtonian plateau.
To proceed, we now explore the growth rates of sinusoidal perturbations to the
interface. Mora and Manna (2009) derive the dispersion relation for the displacement
of one generalised Newtonian fluid by another (their equation 35), which we re-express
here:
σ =
k (px,1 − px,2 − k2T )√
px,1
U
dpx,1
dU
+
√
px,2
U
dpx,2
dU
, (3.30)
where px,i is the pressure gradient in the displacing (i = 1) and displaced (i = 2) fluids
respectively. Following their analysis, but treating the displaced fluid as Newtonian
(3.1) while using (3.29) to describe the shear-thinning rheology of the injected fluid,
leads to the dispersion relation
σ =
kU(1− k1
µ2
(
2m+1
m
)m
b2
12
(
2
b
)m+1
Um−1)− τk3
1 + k1
µ2
(
2m+1
m
)m
b2
12
(
2
b
)m+1√
mUm−1
. (3.31)
Following the method in §3.3.1, we proceed by making (3.31) nondimensional. Scaling
lengths with xf and time with tf gives
66 Viscous instability with time-dependent rheology
σˆ =
kˆUˆ(1− νˆUˆm−1)− τ kˆ3
1 + νˆ
√
mUˆm−1
, (3.32)
where νˆ = k1
µ2
(
xf
tf
)m−1 (
2m+1
m
)m
b2
12
(
2
b
)m+1
and τ =
b2Ttf
12µ2x3f
. The hat notation is hence-
forth dropped. The most unstable wavenumber, kmax, with the largest growth rate
σmax is given by
kmax =
√
(1− νUm−1)U
3τ
, (3.33)
where
σmax =
2
3
√
3τ
((1− νUm−1)U)3/2
(1 + ν
√
mUm−1)
. (3.34)
Ghannam and Esmail (1998) observed that over the range of concentrations of
PAA 0.25% < CPAA,wt < 1%, the flow behaviour index of the polymer solution, m,
was constant with the value m = 0.5, whereas the consistency index, k1, increased
from 410 to 1800 mPa s0.5. Thus, we investigate the case in which k1(t) varies with
changes in concentration following the release of polymer from an encapsulant but m
is constant. Generally, however, we note that either m or k1 might vary with the
polymer concentration. For direct comparison to §3.3, we suppose that ν(t) has the
same functional form as V (t) in (3.13). In this case, we note that under the condition of
constant flow, U = 1, the pressure gradient, px,1, would be the same if the displacing
fluid was Newtonian – with viscosity µ1 – or non-Newtonian, with the value of k1
prescribed by the functional form of ν. The evolution of the perturbation amplitude,
A
√
τ = e
∫ √
τσmaxdt, during an injection process with a constant flow rate is shown with
the solid lines in figure 3.10 for V0 = 0.01, Vf = 0.9 and θ = 10. For reference, we
include the Newtonian case (black), and a series of shear-thinning cases: m = 0.99
(blue), m = 0.95 (red) and m = 0.78 (yellow). The figure illustrates that the effect of
shear-thinning is to cause a small increase in the amplitude of perturbations when the
unperturbed pressure gradients are the same as in the Newtonian case. We now seek
3.5 Injection of a Shear-Thinning Fluid 67
 
 
A∗
√
τ
A
√
τ
t
0 0.5 1
1
1.01
1.02
1.03
Fig. 3.10 Evolution of amplitudes for a shear-thinning fluid in which the apparent
viscosity ratio varies according to V0 = 0.01, Vf = 0.9 and θ = 10 for constant injection
(U = 1; solid lines) and a series of values of m that decreases from the Newtonian case
(m = 1; black) at bottom of the figure upwards: m = 0.99 (blue), m = 0.95 (red)
and m = 0.78 (yellow). Also included are the evolution of amplitudes in the case of
optimal injection (cf. figure 3.11(a) broken lines).
to minimise the final amplitude of perturbations in the case in which a finite volume
of fluid is injected in finite time.
Following the previous sections, substitution of (3.34) into (3.7) leads to the scalar
ODE for U(t)
ν(t)4U(t)5m−1
5∑
i=1
[
ν(t)−iU(t)−i(m−1)
(
Ξ(m)i1ν˙(t)U(t) + Ξ(m)i2ν(t)U˙(t)
)]
= 0, (3.35)
where Ξ(m) =


−m(2 +m) −m2(2 +m)
−10√m− 5m− 2m3/2 + 2m2 2(m−√m)(m2 − 4m− 10√m− 4)
−√m(−14−√m+ 8m+ 4m3/2) −√m(4 + 9√m− 24m− 7m3/2 + 8m2 + 4m5/2)
(3 + 2
√
m)(1 + 2m) 2
√
m(2m2 + 3m3/2 − 5)
0 −3


.
68 Viscous instability with time-dependent rheology
U
∗
t0 0.5 1
0
0.4
0.8
1.2
1.6
(a)
 
 
A
∗√τ
f
A
√
τ
f
θ
10−1 100 101 102
1
1.05
1.1
(b)
Fig. 3.11 Injection with Vf = 0.9. (a) Optimal injection flow rate profiles. (b) Variation
of the final amplitude with gelling speed parameter θ for constant injection (solid lines)
and optimal injection (broken lines). The values of m in the cases of shear thinning
correspond to those in figure 3.10.
The solution of (3.35) is shown in figure 3.11(a) form = 1 (black broken line, cf. figure
3.1(b)) and a series of decreasing values of m: m = 0.99 (blue), m = 0.95 (red) and
m = 0.78 (yellow) and the variation of ν(t) defined by (3.13).
When the injection rate is optimal, the proportion of fluid injected later in the flow
is reduced compared to the Newtonian case. This is interpreted to be a response to
the shear-thinning that would be enhanced with a faster flow rate and which would
lead to a more unstable front at later times. This is analogous to the modification
of the optimal injection profile for air invading a non-gelling power law fluid in an
axisymmetric geometry described by Fontana et al. (2014). The broken lines in fig-
ure 3.10 show the evolution of the largest amplitude, A
√
τ , resulting from injection
with a variable injection rate as shown in figure 3.11(a). As the flow behaviour index
decreases, the final amplitude increases owing to the greater impact of the shear thin-
ning. Figure 3.11(b) shows the variation of the final amplitude with the gelling speed,
θ, after the optimal injection as shown in figure 3.11(a) (broken lines) in comparison
to the case of constant injection (solid lines). The black lines — which correspond to
Newtonian flow — can be recognised from figure 3.2(b), which should be referred to
in order to see the behaviour at large or small values of θ where the lines converge
3.6 Discussion 69
(since if there is effectively no viscosity change then constant injection is optimal). The
effect of shear-thinning is to reduce the degree of suppression of the amplitude. For
example, for θ = 5 and m = 0.78, the degree to which the amplitude is suppressed is
approximately halved as compared to the Newtonian case.
3.6 Discussion
In this chapter we have explored the use of the Euler-Lagrange variational framework
to find the particular injection rate which leads to a minimisation of the growth of the
modes in the Saffman-Taylor instability. We first explored the problem in a rectilinear
geometry and found that if the viscosity of the injected fluid gradually increases with
time, then it is optimal to increase the injection rate gradually with time, so that more
of the injection occurs during the period when the system is more stable.
In a radial geometry the problem is more complex since the azimuthal modes
are quantized and the curvature of the interface tends to stabilize the system. As
the radius increases, with a constant injection rate, the system becomes unstable to
progressively higher modes leading to a shift in the mode of highest amplitude as a
function of time. This leads to a difference in the optimal injection strategy depending
on the average injection rate over time. With a relatively slow mean injection rate,
the optimal strategy involves a gradual decrease in the injection rate with time, consis-
tently with earlier predictions (Beeson-Jones and Woods, 2015; Cardoso and Woods,
1995), whereas with a faster injection rate, the optimal strategy involves a gradual
increase in the injection rate with time (Dias et al., 2012), so that less time is spent
at conditions near the maximum growth rate of each mode. In this way, the overall
amplitude of the perturbations can be reduced.
In the radial system, if the viscosity of the injected fluid gradually increases with
time, then the variational calculus suggests that the optimal injection rate should
increase at a faster rate with time than the optimal rate in the constant viscosity case,
so that more of the fluid is injected once its viscosity has increased. This involves
70 Viscous instability with time-dependent rheology
commencing the injection more slowly but then gradually building up the injection
rate with time beyond the case of constant viscosity.
Although these results are based on linear stability theory, they point to the signif-
icant benefits of controlling the injection rate in order to control the growth of viscous
instability and of the merits in optimising the deployment of a polymer if the injected
fluid is changing in viscosity with time. There are several interesting developments of
this approach which merit further analysis. First, in some cases it may be that there
is a background level of noise, and in this case, the amplitude of each mode would
remain at least as large as this noise threshold for all time. It would be fascinating
to explore the effect of such a forcing on the present problem. Secondly, the present
analysis only explores the linear phase of the growth of an instability. In many cases,
the instability may grow to have nonlinear effects, and it is not clear how the optimal
injection strategy should evolve once the instability has become nonlinear. To this end,
Huang and Chen (2015) and Chen and Yan (2017) have shown that injection with a
linearly increasing injection rate can suppress instability as the amplitude increases
into the nonlinear regime for flow both in a Hele-Shaw cell and a porous medium.
However, it is not clear whether a variable injection rate will influence the ultimate
development of a fractal pattern of fingering (e.g. Praud and Swinney (2005)), and
this would be of interest to explore in future work.
Chapter 4
The dynamics and control of
late-stage viscous fingering
4.1 Summary
The nonlinear evolution of patterns formed during the injection of an inviscid fluid into
a viscous fluid in a Hele-Shaw cell is analysed through annularly averaged saturation
distributions. The displacement begins from close to the injection point rather than at
a large distance from it. Three injection conditions are considered: constant flow rate,
linearly increasing flow rate and constant inlet pressure. For comparison, we perform
a series of off-lattice Diffusion Limited Aggregation (DLA) numerical simulations. We
show how the saturation distributions can be divided into a growing region, where
pattern development is concentrated, and a frozen region in which there is almost
no growth. For viscous fingering, there is additionally a fully saturated region in
the vicinity of the inlet. In the frozen region, the saturation approximately varies
according to the fractal dimension D = 1.70. In the growing region, the saturation
decreases linearly with radius until reaching zero at the outer edge of the pattern and
also depends on D. For DLA, the maximum radial extent of the pattern is 1.54± 0.08
times larger than that of the frozen region throughout development. For VF, this
ratio is found to be 1.40 ± 0.07 following a period of evolution. This observation
72 Dynamics of Fractal Growth
provides a basis to predict the evolution of the radii of the pattern. The extents of
the overall pattern and the frozen region relative to the extent of the inner-most fully-
saturated region are found to be approximately independent of the injection scheme
and correspond to those values for DLA. However, the extent of the fully-saturated
region does depend on the injection scheme.
4.2 Introduction
4.2.1 Viscous Fingering and the Confinement of Growth
When water is used to displace oil during the recovery of oil from subsurface oil
reservoirs, the interface between the fluids can become subject to viscous fingering
(VF). In VF, fingers of the less viscous fluid invade the more viscous fluid (Hill, 1952;
Saffman and Taylor, 1958). The problem of predicting the behaviour of the interface
has been approached by a variety of methods including the theoretical (Saffman and
Taylor, 1958), the experimental (McCloud and Maher, 1995) and the computational
(Li et al., 2007). Whilst a number of studies have explored porous media flows, the
Hele-Shaw cell is commonly used as a simplified analogue to porous media to make
progress in understanding the basic mechanisms of VF (Homsy, 1987). Whilst the
stability of an interface can be found from linear stability analysis (LSA; Saffman
and Taylor (1958)), as the instability develops the fingers become large in amplitude
and nonlinear calculations are needed in order to determine the shape of fingers. The
behaviour of nonlinear fingers in a Hele-Shaw cell is controlled by the geometry of the
confinement imposed by the cell. Lajeunesse and Couder (2000) explored the impact
of the confinement. If, for example, the finger is a significant width of a rectilinear
channel, confinement is strong and one finger with a stable shape forms (Saffman and
Taylor, 1958). The shape of that finger is well understood (Tanveer, 1987). When
the confinement is weak, the shape of the finger is unstable and the tip-splitting
instability occurs. This tip-splitting behaviour is seen in an axisymmetric geometry
(Couder, 1988), in which the fingers behave as if they are in virtual wedge-shaped
4.2 Introduction 73
cells (Thomé et al., 1989). The nature of heavily branched nonlinear viscous fingering
is too complex for the exact shape of fingers to be determined, thus bulk properties
are used to describe the pattern. As we shall describe, the fractal dimension – D –
is one such bulk property. However, reports of the fractal dimension in the literature
are inconsistent. Before going into more detail on the various reports of the fractal
dimension as applied to viscous fingering, we shall first describe what a fractal object
is and the various methods of determination of the fractal dimension.
4.2.2 Fractal Objects, Diffusion Limited Aggregation and the
Determination of the Fractal Dimension
If one were to measure the area of a square with a number of measuring squares, N ,
then if the side length of the measuring squares, ǫ, is halved, the number of squares
needed to cover that space would quadruple according to the relationship N ∼ ǫ−2.
This is because the object is two dimensional. In contrast, the Sierpinksi gasket set
is comprised of an equilateral triangle that is then recursively subdivided into three
equally-sized smaller equilateral triangles (Meakin, 1998). If this shape’s area were
to be measured with a number of measuring triangles, N , of side length ǫ, then as
the length of the measuring triangle is halved, it would be found that three times the
number are required to most precisely measure the area. This finding would follow
the relationship N ∼ ǫ−log(3)/log(2). In comparison to the square, the spatial dimension
is fractional and is therefore referred to as a fractal dimension. In this context, the
fractal dimension is defined as D in the relationship
N ∼ ǫ−D. (4.1)
The box-counting method is one method of determining the fractal dimension of a
pattern (Foroutan-pour et al., 1999). In the box-counting algorithm, the number of
boxes, N , of a particular size, ǫ, required to cover (or cover the edges of) an object is
determined. Then, the side length of the boxes is decreased, and the new number of
74 Dynamics of Fractal Growth
boxes required to cover the shape is determined. The box-counting fractal dimension
is found by evaluating −dlog(N)/dlog(ǫ). In the cases of the Sierpinski gasket and the
square, this dimension would be measured at all length scales since they are infinitely
self-similar sets. In contrast, in real-world objects, as the measuring boxes become
smaller than the smallest length scale present in the object, such self similarity will
cease. Real-world objects may be described as fractal if they have a consistent fractal
dimension – as measured in this way – over a range of length scales. One example of
a self-similar fractal object is that which results from Diffusion Limited Aggregation
(DLA), a computational process that gives rise to a consistent fractal dimension down
to the smallest length scales present in the simulation, as we shall now discuss.
In DLA, a particle is released and randomly walks until it collides with and then
sticks to a growing cluster of the particle’s predecessors, at which point the next par-
ticle is released (Witten and Sander, 1981). The process can produce patterns that
are strikingly similar in appearance to VF (Praud and Swinney (2005)). Daccord
et al. (1986) describe the mathematical connections between DLA and viscous fin-
gering as follows. In the case that the fluids are incompressible and the displacing
fluid is of negligible viscosity, then the pressure field, p, in the displaced fluid obeys
a Laplace equation, ∇2p = 0. In an analogous fashion, the local probability field of
the randomly-walking particle in DLA obeys this equation. Mathiesen et al. (2006)
further explore mathematical similarities between the two processes, and find them to
belong to the same "scaling universality class". Praud and Swinney (2005) highlight
the various methods that have been used to simulate DLA and list the fractal dimen-
sion measurements. The most precise two results they mention are from the work of
Ossadnik (1991) and Davidovitch et al. (2000), who used an off-lattice simulation and
conformal map theory to find D = 1.712 ± 0.003 and D = 1.713 ± 0.003 respectively.
Witten and Sander (1981) stated that the number of particles, Np, enclosed by a radius
R is
Np ∼ RD. (4.2)
4.2 Introduction 75
This is closely related to (4.1) since it is also a law that describes how the measured
pattern area varies with the scale of observation in self-similar objects. Witten and
Sander (1981) showed that DLA obeys this law and therefore that the aggregate is a
fractal object. On account of the connections to VF, attention has been given to assess
the fractal dimension of VF. We shall now introduce the main methods to evaluate
the fractal dimension, before going on to discuss studies that applied these methods
to VF.
A variety of different methods have been established to estimate the fractal di-
mension, including: (i) the box-counting method (Foroutan-pour et al., 1999), as
mentioned before, (ii) correlation-correlation (used by Witten and Sander (1981)),
(iii) comparing the radius of gyration of the pattern to the area (used by May and
Maher (1989)) and (iv) comparing the density, Np/R2, of the pattern to the total
area (used by Chen (1989)). In the case of the patterns formed in their experiments,
Praud and Swinney (2005) found that methods (i) and (ii) are equivalent, whereas (ii)
and (iii) were found to be approximately equivalent by Daccord et al. (1986). Method
(iv) has not explicitly been checked for consistency with the other methods, and has
only directly been used by Chen (1989). This method follows from (4.2) if a particle
in DLA is equivalent to a small unit of area in VF. We shall proceed by assuming
that these methods are at least approximately equivalent, although there will be more
discussion on this subject later in the chapter. We will now give a brief overview of
three key studies that have investigated the weakly confined VF regime (Chen, 1989;
May and Maher, 1989; Praud and Swinney, 2005). We shall discuss the work done to
characterise the final VF pattern, before going on to discuss dynamical aspects.
4.2.3 The Fractal Dimension of Viscous Fingering Patterns
May and Maher (1989) investigated weakly confined VF patterns in the case of a
constant injection flow rate, whereas Praud and Swinney (2005) investigated the case
of a constant inlet pressure, which leads to a gradually increasing injection flow rate.
In both cases, the value of D was found. Praud and Swinney (2005) injected air at very
76 Dynamics of Fractal Growth
large flow rates, to the extent that an asymptotic pattern formed. Praud and Swinney
(2005) used the box-counting method to measure the patterns’ fractal dimension as
D = 1.70 ± 0.02 across all the experiments over a decade of box sizes. This finding
is consistent with the fractal dimension found by Ossadnik (1991) for DLA. Praud
and Swinney (2005) state that it is only in this asymptotic case of a very large flow
rate that the value of D is expected to converge to this value. Praud and Swinney
(2005) also stipulated that the wavelength of the instability must be constant during
the injection in order that the system converge to a DLA-like fractal dimension. This
condition was met in their experiments since the constant pressure condition gives
rise to an approximately constant tip velocity (Thomé et al., 1989), and therefore a
constant wavelength of instability. To compare these findings to those of May and
Maher (1989), we define the non-dimensional flow rate:
C =
µQ
σb
, (4.3)
where µ is the viscosity of the displaced fluid, Q is the injection flow rate, σ is the
interfacial tension and b is the plate separation. May and Maher (1989) investigated
the range of flow rates 0 < C < 60 whereas Praud and Swinney (2005) investigated
100 < C¯ < 700, where C¯ is the mean flow rate. May and Maher (1989) found the value
of the fractal dimension of the fingering pattern was 1.71± 0.04 for their experiments
with injection flow rates in the range 10 < C < 40. This is inconsistent with the claim
of Praud and Swinney (2005), that very large flow rates, C, are required to bring
about this value of the fractal dimension and surprising since at such modest flow
rates the action of surface tension in shaping the fingers is apparent by inspection of
the patterns. Furthermore, with this constant injection flux, the tip velocity gradually
slows: it is not constant. In the range of injection flow rates 40 < C < 60, May and
Maher (1989) found the fractal dimension slowly increased to the value D = 1.79±0.04.
We note that May and Maher (1989) measured the fractal dimension by comparing
the radius of gyration to the area of the pattern, and assumed that the pattern was
fractal on account of the power-law relationship between these two properties.
4.2 Introduction 77
Like May and Maher (1989), Chen (1989) investigated the case of constant injection
flow rate. He defined the pattern density, at a certain radius, as the fraction of the
total area contained in that radius that is occupied by air. The fractal dimension was
found by measuring the gradient of the region in which there was a linear relation
between log-density and log-radius. It was stated that a number of generations of
tip splitting are required to bring about fractal behaviour, and this explained why
only the experiment with the largest flow rate, C = 19, exhibited a linear region in
the density-radius plot. For this experiment, he found from successive plots that the
fractal dimension decreased from D = 1.9 to D = 1.8 over the course of the injection.
This is inconsistent with the data of May and Maher (1989), who found the fractal
dimension to be 1.71 ± 0.04 for this value of C. Furthermore, the observation that
the dimension evolves with time is inconsistent with the observation of Praud and
Swinney (2005) that growth of the pattern is concentrated at the tips whilst the inner
portion of the pattern is frozen, and provides a channel for the fluid to reach the outer
growing part of the domain. Chen (1989) found a scaling law for the time evolution
of the overall radius of the pattern, but he did not attempt to describe the extent of
the fractal region or the behaviour of the mass distribution at larger radii.
In summary, all these previous studies were centred on the description of viscous
fingers in terms of a fractal dimension of the fingering pattern and the results are
inconsistent. The rate of advance of the leading edge of the fingering pattern, the
variation of the saturation in this leading part of the flow, and the relation between
the radius of the leading edge and the radius of the region in which the fingers follow
a fractal dimension has had much less analysis. This forms the main content of this
chapter. To this end, we re-analyse the published data and we carry out a series
of new experiments. In the case of constant pressure injection, Praud and Swinney
(2005) kindly shared the time series from their experiments with us, and so we have
examined this data as well.
The third interesting class of inlet condition corresponding to a linearly increasing
injection rate is motivated by the work of Dias et al. (2012). By applying the constraint
78 Dynamics of Fractal Growth
of injecting a finite volume of fluid in a finite time, Dias et al. (2012) used variational
calculus to find the injection flow rate profile that minimises the final amplitude of
perturbations as predicted by linear stability analysis. The result of this calculation
was an injection flow rate that linearly increased in time. Injecting in this manner
was demonstrated to be an effective means of suppressing the growth of the instability.
However, the examples shown in Dias et al. (2012) involved the displacement starting
from a radius at some distance from the central inlet. If the starting radius is very
small all modes are predicted to be initially stable by linear stability analysis. It
has been observed that the first mode to become unstable controls the initial growth
whereas if the starting radius is larger a number of modes are linearly unstable and it is
the one with the highest growth rate that controls the shape (Martyushev et al., 2015).
The optimal injection flow rate profile for a very small initial radius is described in the
paper of Dias et al. (2012). Whilst the optimum is not necessarily a linearly increasing
profile in this case, it is nonetheless of interest to compare the mass distribution that
results from such a displacement to the other injection rates to assess the degree to
which the instability is suppressed by this means when the initial radius is very small
and the fingers undergo tip splitting. Bischofberger et al. (2015) provided us with
data to add to our own for this inlet condition.
Finally we have carried out a similar analysis of the mass distribution of the pat-
terns that arise in DLA given the previously mentioned similarities. We perform DLA
calculations and compare results with the new analysis of laboratory experiments.
4.2.4 Chapter Outline
In §4.3, we give an overview of the new laboratory experiments we carried out to
model pattern evolution for the constant and linearly increasing injection flow rates.
We also include the data shared with us by Praud and Swinney (2005). In §4.4 we
illustrate the experimental patterns formed. Subsequently, we find and compare the
mass distributions for the various VF experiments and seek a universal model. We
then go on to perform a series of off-lattice diffusion limited aggregation calculations
4.3 Experiments 79
in §4.4.4 for comparison. In §4.5 we attempt to put our results in the context of other
VF studies. Finally, in §4.6 we draw some conclusions.
4.3 Experiments
4.3.1 Apparatus, Experimental Conditions and Protocol
In our constant flow rate and linearly increasing ("ramped") flow rate experiments,
air was injected into rapeseed oil in a float glass Hele-Shaw cell. The top plate was
1 cm thick and 60 cm in diameter and the bottom plate was 1 cm thick and 55 cm
in diameter. Float glass was chosen as the process of toughening glass deforms it,
and this was seen to affect early (not included) experiments. The physical properties
of the rapeseed oil is shown in table 4.1. Viscosity was measured using a Bohlin
CS-50 rheometer, and surface tension with a Kruss DSA 100 tensiometer Following
the experimental design of May and Maher (1989), the gap width was set with three
positioning micrometers. Air was injected though a 1 cm hole in the centre of the top
plate via two syringe pumps (World Precision Instruments AL1000-220). The Hele-
Shaw cell was immersed in a reservoir of rapeseed oil. Experiments labelled C1-C3 and
R1-R3 were performed by us with this apparatus. The top plate was painted white to
help visualisation, and the pattern was filmed with a Canon D90 digital camera from
underneath such that the view was not obstructed by delivery tubes.
Experiments from Praud and Swinney (2005) are included in table 4.1 (keys P1-P3).
Praud and Swinney (2005) performed constant pressure experiments and achieved non-
dimensional flow rates at least an order of magnitude larger than other studies, by using
6 cm thick optically polished glass to attain a uniform 0.127 mm gap. This thickness
prevented deflection even at the large values of forcing involved in their experiments.
They injected air into silicone oil.
Also included are two previously unpublished experiments provided to us by Bischof-
berger et al. (2015) (keys C4 and R4). Those authors used two glass plates of thickness
1.9 cm and radius 14 cm.
80 Dynamics of Fractal Growth
Table 4.1 Experimental setup. Experiments C1-C3 and R1-R3 were performed by the
authors of this chapter, whereas P1-P3 were performed by Praud and Swinney (2005)
and C4 and R4 were performed by Bischofberger et al. (2015).
Type Key b [mm] ∆P [atm] Fluids µ [Pa s] σ [N/m]
Constant flow C1 1.5 air/rapeseed oil 0.050 0.03
Constant flow C2 1.5 air/rapeseed oil 0.050 0.03
Constant flow C3 1.0 air/rapeseed oil 0.050 0.03
Constant flow C4 0.254 water/silicone oil 0.297 0.027
Constant pressure P1 0.127 0.25 air/silicone oil 0.345 0.021
Constant pressure P2 0.127 0.50 air/silicone oil 0.345 0.021
Constant pressure P3 0.127 1.25 air/silicone oil 0.345 0.021
Ramped flow R1 1.5 air/rapeseed oil 0.050 0.03
Ramped flow R2 1.0 air/rapeseed oil 0.050 0.03
Ramped flow R3 1.0 air/rapeseed oil 0.050 0.03
Ramped flow R4 0.127 water/silicone oil 0.847 0.027
4.4 Results
In this section, we present the measured areal flow rates and then go on to describe
and compare the patterns formed. We then develop a mass distribution model for
these patterns. Finally, we perform a series of DLA simulations and compare the
mass distribution to viscous fingering.
4.4.1 Measured Flow Rates
A MATLAB script was written for this project to detect the boundaries of the evolving
pattern and extract properties such as the area enclosed and furthest distance to a
finger tip. The sequence of MATLAB functions used are included in an appendix to
this chapter (§4.7). The boundary detection was monitored at each step.
Figure 4.1 shows the evolution of the pattern area over the course of the exper-
iments. In each case, we use least-squares regression to fit the data to a power-law
relation, A = A1tA
′
, shown as the broken black lines. The measured values of A1 and
A′ are shown in table 4.2 with 95% confidence intervals shown in brackets. The injec-
tion flow rate is determined from this best-fit line as Q = dA/dt = βtα. The values of
β and α are also shown in table 4.2. Table 4.2 also shows the time-averaged measured
4.4 Results 81
Table 4.2 Measured area and flow rate.
Key A1 A′ β α Q¯ Q¯p C¯
[cm2/sA
′
] [cm2/s(1+α)] [cm2/s] [cm2/s]
C1 112(109,114) 1.00 (0.99, 1.01) 110 0.0 112 110 12
C2 186(184,188) 0.97 (0.96, 0.97) 180 0.0 183 180 21
C3 245(241, 249) 0.95 (0.92, 0.97) 230 0.0 236 230 40
C4 7.97(7.94, 8.00) 1.00 (1.00, 1.01) 8.0 0.0 8.0 8.0 35
P1 1.81(1.69,1.93) 1.50 (1.47, 1.53) 2.7 0.5 6.5 (-) 110
P2 5.20(4.90, 1.39) 1.43 (1.34, 1.46) 7.4 0.4 11 (-) 180
P3 16.9(16.6,17.1) 1.60 (1.57, 1.62) 27 0.6 27 (-) 440
R1 10.2(9.5,10.9) 2.10 (2.05, 2.15) 21 1.1 69 70 7.8
R2 28.6(27.0,30.1) 1.83 (1.78, 1.87) 52 0.8 87 85 14
R3 26.2(24.0,28.4) 2.01 (1.94, 2.08) 52 0.9 95 110 16
R4 0.27(0.28,0.28) 2.02 (2.01, 2.02) 0.55 1.0 2.8 (-) 69
flow rate, Q¯, and the time-averaged flow rate with which the injection pumps were
programmed, Q¯p, in the cases that it is known. Finally, we present the nondimensional
flow rate. The value of C is constant for the constant flow rate experiments, whereas a
mean value – C¯, based on Q¯ – is presented for the experiments that feature a gradually
changing flow rate.
For the experiments with a constant flow rate, the variation of area is approximately
linear. However, for the experiments with larger flow rates, the time exponent becomes
progressively smaller than 1, suggesting there is either a degree of plate lift or air
compression in these cases. Experiments C1 and C4 show excellent agreement with a
linear fit. The agreement between the measured and programmed time-averaged flow
rates is good (within 3%) in all cases.
Figure 4.1(b) relates to the experiments of Praud and Swinney (2005) in which air
was injected at constant pressure. The best-fit lines suggest an approximate power-law
time dependence of flow rate, with α = 0.5 ± 0.1. The programmed flow rates are
unknown.
The experiments featuring a linearly increasing flow rate are shown in figure 4.1(c).
The variation of pattern area is approximately parabolic. In experiment R3 there is
82 Dynamics of Fractal Growth
C1
C2
C3
t[s]
A
[c
m
2
]
0 2 4 6 8
0
250
500
750
1000
(a) Constant flow rate.
P1
P2
P3
t [s]
A
[c
m
2
]
0 5 10 15
0
25
50
75
100
(b) Constant pressure (Praud and Swin-
ney, 2005).
R1
R2
R3
t [s]
A
[c
m
2
]
0 2 4 6
0
250
500
(c) Ramped flow rate.
R4
C4
t [s]
A
[c
m
2
]
0 2 4 6 8 10
0
10
20
30
40
50
(d) Experiments provided by Bischof-
berger et al. (2015)
Fig. 4.1 The evolution of the pattern area, A.
considerable deviation (14%) in the measured flow rates from the programmed value,
suggesting that plate lift is more severe than in other cases.
4.4.2 VF Pattern Formation
Figure 4.2 illustrates the development of the viscous fingering pattern for a selection
of experiments. As was first noted by Paterson (1981), the bases of the fingers – at
the bottom of fjords between fingers – become stationary. Conversely, finger growth
occurs at the tips. Imperfect tip splitting is observed. This is where one branch grows
4.4 Results 83
at the expense of the other branch, rather than permitting the growth of two long-lived
structures (Lajeunesse and Couder, 2000). In the case of the constant flow experiments
(top row), there is a larger number of thinner fingers for C3 as compared to C1 and
C2, which both feature similar patterns. This could be on account of quantization;
that the difference in the value of C is not sufficient to make the next integer mode
the most unstable at the time nonlinear growth begins to dominate.
The patterns for the high forcing, constant pressure experiments of Praud and
Swinney (2005) (middle row) show a highly branched structure. The width of the
fingers appears approximately constant throughout each pattern.
The ramped flow experiments (bottom row) are a return to the less forced patterns;
fingers are wider and there are fewer tip-splitting events. The radial extent of the
patterns is approximately half that of the constant flow series since the large flow rates
involved towards the end of experiments with larger patterns led to plate deformation.
The patterns for R2 and R3 are similar and feature a larger number of thinner fingers
as compared to the R1 case, for which there is also a large central region that has no
fingers.
4.4.3 An Empirical Saturation Model
We define the saturation, S(R), as the fractional area of an annulus at radius R that is
taken up by the pattern. The radial variation of S is plotted as a time series in figure
4.3 for three experiments: the constant flow rate experiment for which C = 12 (C1),
the Praud and Swinney (2005) experiment for which ∆P = 0.25 atm (P1) and the
ramped experiment for which C¯ = 14 (R2). At each time, there is a small inner fully
saturated region, outside of which S gradually decreases. At a particular radius, Rfr,
there is an adjustment from the gradual decrease of S with radius to a region in which
there is a steeper and approximately linear decrease with radius until reaching the
point S = 0. This corresponds to the edge of the growing pattern. For radii smaller
than Rfr the profile is fixed in time, thus Rfr approximately delimits a frozen region
84 Dynamics of Fractal Growth
−20 −10 10 20
−20
−10
10
20
(a) C1, 1s
−20 −10 10 20
−20
−10
10
20
(b) C2, 1s
−20 −10 10 20
−20
−10
10
20
(c) C3, 0.5s
−5 5
−5
5
(d) P1, 3.5s
−5 5
−5
5
(e) P2, 2s
−5 5
−5
5
(f) P3, 1s
−10 10
−10
10
(g) R1, 0.4s
−10 10
−10
10
(h) R2, 0.25s
−10 10
−10
10
(i) R3, 0.25s
Fig. 4.2 VF pattern development time series for the experiment keys and time intervals,
δt, as shown in the captions. Scales in [cm].
4.4 Results 85
R[cm]
S
0 10 20 30
0
0.5
1
(a) C1; (1.75s, 1s).
R[cm]
S
0 5 10 15
0
0.5
1
(b) P1; (1.25s, 2.25s).
R[cm]
S
0 10 20
0
0.5
1
(c) R2; (1.42s, 0.38s).
Fig. 4.3 The evolution of saturation for three different examples of inlet condition with
the model (4.5) fitted (black lines). The time series are evenly-spaced in time and are
shown starting at an initial time ti at intervals ∆t, with these values as shown in the
captions as (ti, ∆t).
B CA
Rb Rfr Rtip
R/Rb
S
0
0.5
1
Fig. 4.4 Schematic of the saturation profile.
from an outer region of growth. Figure 4.4 is a schematic of this profile with the fully
saturated, frozen and growing regions labelled as A, B and C respectively.
In §4.2.2 we introduced (4.2), which describes the number of particles contained
by a radius, R, in a DLA cluster (a fractal object). If – analogous to DLA – viscous
fingering also produces fractal patterns then, given (4.2), one expects the saturation
profile to be given by
Sfr =
(
R
Rb
)D−2
. (4.4)
86 Dynamics of Fractal Growth
In (4.4), Rb is a fitting parameter which approximately corresponds to the radius of
the finger bases, such that S = 1 when R = Rb. The upper black line in figure 4.3 is a
fit of (4.4) with the fractal dimension D = 1.70, the value for viscous fingering found
by Praud and Swinney (2005). The complete empirical DLA saturation model is:
S(R, t) =


1, for 0 < R ≤ Rb(
R
Rb
)D−2
, for Rb < R ≤ Rfr(t)(
Rfr
Rb
)D−2 Rtip(t)−R
Rtip(t)−Rfr(t) , for Rfr(t) < R < Rtip(t)
0, otherwise.
(4.5)
The empirical model appears to fit approximately for each of the experiments
shown in figure 4.3, despite differences in the appearances of the patterns and injection
conditions. Figure 4.5 shows one time profile from every experiment in this chapter: (a)
the constant flow rate experiments, (b) the Praud and Swinney (2005) experiments and
(c) the ramped experiments. Table 4.3 also lists the root-mean-squared-error (RMSE)
between the data and the model, with D = 1.70 in the frozen region (RMSEfr). It is
typically around 5% in this region.
The model appears to provide a good fit across almost all of these experiments.
The largest RMSEfr is experiment R1. It can be seen from figure 4.5(c) that the
approximation is poorer in this case. It is at a lower flow rate and there are fewer
splitting events compared to the other patterns. This suggests that tip splitting gives
rise to the gradual decrease in saturation that fits the frozen region of the model. In
most experiments, there is a dip in S below Sfr at small radii; this could be attributed
to the shape of the fjord bases. In figure 4.5(a) there is a kink in the saturation at the
tip for the C2 and C3 experiments; this could be a result of the finite size of the cell
or tips widening ahead of the next generation of splitting.
We will now explore the values of Rb, Rtip(t) and Rfr(t). They are found by fitting
the model (4.5). The variation of the fitting parameter Rb in the various experiments
is shown scaled with the gap thickness b and varying with C¯ in figure 4.6.
4.4 Results 87
C1
C2
C3
C4
R
Rb
S
1 3 5 7
0
0.5
1
(a) Constant flow rate injection: C1-C4 (from left to right).
P1
P2
P3
R
Rb
S
0 25 50
0
0.5
1
(b) Constant pressure injection: P1-P3 (from left to right).
R1
R3
R2
R4
R
Rb
S
1 5 9
0.5
1
(c) Ramped flow rate injection: R1-R3 (from left to right).
Fig. 4.5 The final saturation profile across experiments. In all figures, the broken line
is (4.4) with D = 1.70.
88 Dynamics of Fractal Growth
R1
R2
R3
C1
C2
C3
P1
P2
P3
C4
R4
C¯
r b
/b
100 101 102 103
10
20
30
40
Fig. 4.6 Variation of the fitting parameter Rb with C¯. The experiment key number
increases from left to right.
The value of Rb/b declines for the ramped experiments and the constant pressure
experiments, whereas for the constant flow rate experiments it appears to be approx-
imately constant. It could be that it is cut-off by the presence of the inlet tube for
these experiments. The value of Rb/b is much greater for the experiment R4, provided
by Bischofberger et al. (2015), than for the other ramped-flow experiments.
Turning to Rtip and Rfr, one measure of the structure of the growing pattern is
the ratio Rtip/Rfr, the evolution of which is shown as a function of Rtip/b in figure
4.7. The left panel corresponds to the constant flow rate experiments, the centre
panel corresponds to the Praud and Swinney (2005) experiments and the right panel
corresponds to the ramped experiments. Praud and Swinney (2005) measured the
largest angles that could be made between points on the growing, leading edge of the
pattern and termed these the unscreened angles. They observed the distribution of
these unscreened angles collapse with Rtip/b across their experiments. Following their
observation, we plot the variation of the ratio Rtip/Rfr with the quantity Rtip/b.
For the Praud and Swinney (2005) data (figure 4.7(b)), there is convergence of
Rtip/Rfr across experiments at around Rtip/b = (Rtip/b)c = 400, which corresponds
to the unscreened angle distribution reaching an asymptotic shape at this radius (c.f.
4.4 Results 89
Rtip/b
R
ti
p
/R
f
r
0 100 200
0
1
2
(a) Constant flow rate
Rtip/b
R
ti
p
/R
f
r
0 300 600 900
0
1
2
(b) Constant pressure
Rtip/b
R
ti
p
/R
f
r
0 200 400
0
1
2
(c) Ramped flow rate
Fig. 4.7 The evolution of the ratio Rtip/Rfr varying with Rtip/b for the injection profiles
as labelled (c.f. figure 4.5 for legend keys).
Praud and Swinney (2005) Fig. 12), although experiment R2 subsequently diverges
somewhat. We define (Rtip/Rfr)∞ as the mean value of Rtip/Rfr after the tip radius
has reached (Rtip/b)c. In figure 4.7(b), the value of this ratio is approximately 1.4
throughout the experiments. The values of these measures are summarised in table
4.3 for the different experiments, where the standard deviation of (Rtip/Rfr)∞ is also
shown.
In the case of constant flow (figure 4.7(a)), the ratio Rtip/Rfr builds throughout
experiments with C = 12 (C1) and C = 21 (C2), but appears to reach a plateau
at approximately (Rtip/b)c = 100 for the experiments with C = 40 and C = 35 (C3
and C4 respectively). It would be useful to have more later stage data to confirm the
asymptotic value since experiments C3 and C4 feature somewhat different values of
(Rtip/Rfr)∞.
Figure 4.7(c) shows that for the ramped flow rate experiments, the ratio Rtip/Rfr
tends to a value at around 1.4 for the C¯ = 12 and C¯ = 15 experiments (R2 and R3),
which collapse with Rtip/b at approximately Rtip/Rfr = (Rtip/b)c = 100, whereas the
C¯ = 7.6 (R1) experiment is approximately at this ratio throughout. This ratio also
converges to a value of approximately 1.4 for the experiment in which C¯ = 68 (R4).
The collapse of data with Rtip/b in figure 4.7 might be because an asymptotic state
is reached in a time that depends on the evolution of the pressure field in the displaced
90 Dynamics of Fractal Growth
Table 4.3 The error in the frozen region, RMSEfr, and shape parameters of the model.
Key RMSEfr (Rtip/b)c (Rtip/Rfr)∞
C1 0.03 100 n/a
C2 0.04 100 n/a
C3 0.05 100 1.53± 0.02
C4 0.04 100 1.37± 0.05
P1 0.04 400 1.37± 0.02
P2 0.03 400 1.45± 0.05
P3 0.05 400 1.29± 0.04
R1 0.07 100 1.44± 0.03
R2 0.03 100 1.33± 0.03
R3 0.03 100 1.35± 0.02
R4 0.03 200 1.45± 0.03
fluid, and thus depends on the radius of the Hele-Shaw cell used. This is supported by
the observation that (Rtip/b)c is consistent across experiments performed in the same
cell with the same fluids, noting that Bischofberger et al. (2015) used the same cell but
different fluids in C4 and R4. The ratio (Rtip/Rfr)∞ is similar across experiments. We
define Γ as the mean value from these experiments and find Γ = 1.40 with a standard
deviation of 0.07. We proceed by exploring the hypothesis that this is a constant,
universal to viscous fingering.
The saturation profile in the linear region can now be rewritten as
S(R, t) =
(
Rfr
Rb
)D−2 Γ−R/Rfr
Γ− 1 , (4.6)
for Rfr(t) < R < ΓRfr(t). The model for saturation can be integrated to give the total
injected volume. If the volume of the small central fully saturated region, 0 < R < Rb,
is neglected, the integration gives
∫ t
0
Qdt′ =
∫ Rfr
0
(
R
Rb
)D−2
2πRdR +
∫ ΓRfr
Rfr
(
Rfr
Rb
)D−2 Γ−R/Rfr
Γ− 1 2πRdR
=
βt1+α
1 + α
(4.7)
4.4 Results 91
We nondimensionalise with the radius and time at which the linear stability analysis
predicts mode 2 first becomes unstable, i.e. when viscous destabilisation first over-
comes the surface tension of the fluid:
t0 =

(π(1 + α)
β
)(
b2σ
µβ
)21/(3α+1) , (4.8)
R0 =
(
β
π(1 + α)
)α/(3α+1) (
b2σ
µβ
)(α+1)/(3α+1)
. (4.9)
Performing the integration in (4.7) gives
Rfr
R0
= K
(
t
t0
) 1+α
D
,
where K =
( DΓRD−2b βt1+α0
2π(Γ +D(Γ− 1)(Γ3 − Γ− 1))RD0 (1 + α)
)(1/D)
.
(4.10)
Equation (4.10) predicts the evolution in time of the frozen radius, Rfr, by assuming
Rtip = ΓRfr. In figure 4.8, we compare the evolution of these radii with the prediction
(4.10), shown as the broken line and displayed as text. The evolution of Rtip is shown
in the left hand column and Rfr is shown in the right hand column. The top row
corresponds to the constant injection flow rate experiments, the middle to the constant
pressure experiments, and the bottom row to the ramped injection experiments.
For the constant flow experiments, (a)-(b), the predictions of Rtip and Rfr are ap-
proximately consistent with the data, whereas the prediction of Rfr is only consistent
for experiments C3 and C4, i.e. the experiments for which the ratio Rtip/Rfr ap-
peared to reach an asymptotic state (c.f. figure 4.7). For the less forced experiments
(C1 and C2), the values of Rtip and Rfr are initially greater than would be predicted
and decreases to below the line of prediction.
Figure (c) shows that Rtip approximately evolves with a time exponent of 1.0 for P1-
P3, consistent with statement that constant pressure injection gives rise to a constant
92 Dynamics of Fractal Growth
0.18t0.59
t/t0
R
ti
p
/R
0
105 106 107 108
102
103
104
(a)
0.12t0.59
t/t0
R
f
r
/R
0
105 106 107 108
102
103
104
(b)
0.65t1.0
t/t0
R
ti
p
/R
0
102 103 104
102
103
104
(c)
0.46t1.0
t/t0
R
f
r
/R
0
102 103 104
102
103
104
(d)
0.76t1.17
t/t0
R
ti
p
/R
0
100 101 102 103
100
101
102
103
(e)
0.54t1.17
t/t0
R
f
r
/R
0
100 101 102 103
100
101
102
103
(f)
Fig. 4.8 The evolution of the radii Rtip (left column) and Rfr (right column). Top row:
constant injection, middle row: constant pressure, bottom row: ramped injection (c.f.
figure 4.5 for legend keys; the experiment key number increases from bottom left to
top right).
4.4 Results 93
tip velocity (Praud and Swinney, 2005). Figure (d) shows that Rfr also approximately
evolves with this time exponent.
For the ramped flow rate experiments, (e)-(f), the data approximately agree with
the prediction for the tip radius, Rtip, whereas for the inner radius, Rfr, the series
for which C¯ = 69 (R4) agrees with the predicted power-law dependence only at later
stages of pattern growth and remains still somewhat low.
These empirical predictions are derived from assuming the saturation profile is
universal. Indeed, it can be seen from (4.7) that the final radii Rtip and Rfr are
functions of the volume injected alone and not the specific injection history if it is
valid that the pattern properties (Rtip/Rfr)∞ and D are consistent across experiments.
Before comparing the final values of these radii across experiments, we will first analyse
the saturation distribution of a series of Diffusion Limited Aggregation simulations in
order to include those results in the comparison too.
4.4.4 Comparison to Diffusion Limited Aggregation
Our aim is to quantitatively compare the growth of viscous instability patterns with
that of a Diffusion Limited Aggregation (DLA) cluster. A DLA cluster was simulated
using the algorithm developed by Kuijpers et al. (2014). The authors of that study
shared their C++ code with us, although it required some bug fixes to compile. We
used an overall grid size of 20,000 px, a "distance matrix" grid size of 50 px and a
minimum step size of 1 px. The particles were 1 px in radius. Ten clusters were grown,
with the number of particles released ranging from 105 to 106 in equal increments. The
final pattern of a 106 particle cluster is shown in figure 4.9(a). The radial saturation
profile was computed by counting the number of particle centres that fell within a
specified range of radii. Figure 4.9(b) shows the resulting distribution, computed
using a radial step size of ten times the radius of particles. The mean of seven repeats
of the saturation profile of each size is shown, as the saturation profile from a single
pattern featured significant noise.
94 Dynamics of Fractal Growth
(a)
R
S
0 2500 5000
0
0.05
0.10
0.15
(b)
Fig. 4.9 (a) A DLA cluster of 105 particles. (b) The radial saturation profiles of
DLA clusters of size ranging from 105 to 106 particles in ten equal increments. The
simulation to produce each cluster size was repeated seven times, and that shown is
the mean saturation of all repeats.
4.4 Results 95
Rtip
R
ti
p
/R
f
r
2000 4000 6000
1
1.5
2
(a)
1.57N0.58p
Np
R
ti
p
105 106
2000
3000
4000
5000
(b)
1.02N0.58p
Np
R
f
r
105 106
1000
2000
3000
4000
5000
(c)
Fig. 4.10 Evolution of (a) the ratio Rtip/Rfr, (b) Rtip and (c) Rfr as a function of the
number of particles released.
As with viscous fingering, the saturation distribution features a frozen region in
which the variation gradually declines and an outer region in which the saturation lin-
early decreases to zero. In viscous fingering, a base region was observed in which the
pattern was fully saturated. This region is absent in these DLA simulations. In the
frozen region, R < Rfr, the saturation distribution varies according to S = (R/Rb)D−2.
The variation shows excellent agreement with the fractal dimension reported by Os-
sadnik (1991), D = 1.71, shown as the broken black line, which validates the DLA
algorithm we used. Interestingly, the value of Rb – the fitting parameter that defines
the base of the saturation distribution in the frozen region – is 0.1. Therefore, it is not
until a radius of ten particle radii from the centre of the cluster that the saturation
distribution follows the fractal variation.
Following our analysis for viscous fingering, we find the ratio and values of the
radii Rtip and Rfr in figure 4.10.
Figure 4.10(a) shows that Rtip/Rfr is approximately constant over the range 100 <
Rtip < 400. The mean of this ratio is ΓDLA = 1.54, with standard deviation 0.08. This
value is greater than the value for viscous fingering, Γ = 1.40 ± 0.07; however, the
values are not inconsistent. Figures 4.10(b) and 4.10(c) show that Rtip and Rfr follow
power laws with the same exponent, N0.58p , which is a rearrangement of (4.2), i.e.
Rtip ∼ N1/Dp . We shall now compare these results to those of viscous fingering. In
96 Dynamics of Fractal Growth
the last section, the evolution of the values of Rtip and Rfr approximately collapsed
when those radii were scaled with R0 and time with t0. These quantities respectively
correspond to the radius and time at which linear stability theory predicts that the
interface becomes unstable. These constants do not exist for DLA, so instead we
choose to scale Rtip and Rfr with the only scale of length present, Rb. The results are
shown in figure 4.11, where the values are shown as a function of the scaled area of
the pattern, A/R2b . The broken black lines correspond to the lines of best fit for DLA
for Rtip (figure 4.11(b)) and Rfr (figure 4.11(a)), whereas the data points correspond
to those radii for various viscous fingering experiments.
The agreement between the experiments and the simulations is remarkably good,
despite the noted differences in the values of D and (Rtip/Rfr)∞. Experiment R4 is
an exception – the radii only agree with the DLA values at later stages in the flow.
That the agreement in figure 4.11 is better then when lengths and time were scaled
with R0 and t0 is probably on account of experimental uncertainty in those constants.
Aside from experiment R4 at early times, to a close approximation the saturation
profile is fully defined by D, (Rtip/Rfr)∞ and Rb, and only the latter – Rb – depends
on the injection history (see figure 4.6). Otherwise, the saturation profiles for viscous
fingering are independent of injection history and approximately the same as for DLA.
4.5 Discussion
In the experiments presented for which the flow rate is sufficiently high, the variation
of saturation in the frozen region is seen to give approximate agreement with that
variation in DLA. We will now discuss the conditions under which such behaviour is
expected and briefly analyse a number of patterns from other studies in the literature.
Finally, we will discuss the methods of measuring the fractal dimension.
Praud and Swinney (2005) suggested a necessary but not sufficient prerequisite for
VF and DLA to agree in fractal dimension is that the smallest length scale – the wave-
length of instability as predicted by linear stability analysis (λLSA = πb
√
σ/(µUtip))
4.5 Discussion 97
C1
C2
C3
P1
P2
P3
R1
R2
R3
C4
R4R4
0.37(A/R2b)
0.58
A/R2b
R
in
/R
b
100 101 102 103
100
101
(a)
R4
0.55(A/R2b)
0.58
A/R2b
R
ti
p
/R
b
100 101 102 103
100
101
(b)
Fig. 4.11 Comparison of the radii Rin (a) and Rtip (b) for the viscous fingering exper-
iments (points) and DLA (broken lines).
98 Dynamics of Fractal Growth
– is constant throughout the pattern during growth, just as the diameter of random
walkers in a growing DLA cluster is unchanging. For VF, they stated that λLSA is
constant if the tip velocity is constant. In fact, if λLSA is on the scale of the plate
separation – b – then the actual wavelength of the VF pattern – λ – becomes a fixed
multiple of b rather than following λLSA (Dias and Miranda, 2013b; Kim et al., 2009;
Maxworthy, 1989). Indeed, in the case of miscible displacements, for which surface
tension is altogether absent so λLSA → 0, the width of fingers is fixed by viscous
dissipation effects to λ ∼ 4b (Paterson, 1985). Figure 4.12(a) shows the evolution of
the tip velocity on the left hand axis and the capillary number, Ca = Utipµ/σ, on
the right hand axis for the Praud and Swinney (2005) data. Figure 4.12(a) shows
that Utip is indeed approximately constant (until the end of the flow when boundary
effects become important). However, it is also true that λLSA/b = π/
√
Ca ranges
from just 4 to 10, so it is quite possible that the pattern is in the regime where the
wavelength of the pattern is fixed by the plate separation and is independent of the tip
velocity. This is further supported by the similarity of the patterns across the Praud
and Swinney (2005) experiments. Indeed, in figure 4.2, the complexity of the patterns
and the width of the fingers do not appear to change considerably for the different
inlet conditions. In this regime, it is quite possible that any inlet condition (constant
or linearly increasing) could give rise to the same pattern; so long as a threshold tip
velocity is maintained it need not be constant for the unstable wavelength, λ, to be
fixed.
We now turn to the constant flow rate experiments of this study. The evolution
of Utip and Ca for the constant-flow experiments is shown in 4.12(b) for experiments
C1 (blue), C2 (red) and C3 (yellow). The value of λLSA/b ranges from 10 to 20, so
the patterns are probably not in a regime where the finger width is controlled by the
plate separation alone. The tip decelerates, violating the condition that λLSA need
be constant. Therefore, DLA-like fractal behaviour is not expected by the arguments
introduced above, and yet figure 4.5 shows that the variation in the frozen region is
reasonably well approximated by D = 1.70, typically within 5% (see table 4.3). Before
4.5 Discussion 99
CaUtip
t[s]10
−1 101
0
0.825
1.65
0
5
10
(a)
C3
C2
C1
C
a
U
ti
p
t[s]
0 4 8
0
0.1
0.2
0.3
0
6
12
18
(b)
Fig. 4.12 The evolution of Utip [cm/s] (left axis) and Ca (right axis). (a) Praud
and Swinney (2005) experiments (bottom to top: P1 to P3). (b) Constant flow rate
experiments C1-C3.
addressing this, we turn to other constant flow rate experiments to see if they too are
well approximated by the saturation profile for D = 1.70. Figure 4.13 compares the
saturation profiles of the final patterns in a variety of constant flow rate studies to the
DLA profile (black broken lines).
The saturation profile for the VF pattern illustrated in May and Maher (1989) is
included in figure 4.13 as the black series. By comparing the radius of gyration to the
area of the pattern, May and Maher (1989) reported a fractal dimension of D = 1.8 for
this experiment. It can be seen from the figure that the variation in the frozen region is
not inconsistent with D = 1.70, however the frozen region is not very well established.
Thus, it could be that the variations in the measured fractal dimension seen in May
and Maher (1989) are found because the patterns are yet to reach an asymptotic state
in which a more substantial proportion of the pattern is frozen, rather than simply
because the injection is of constant flow rate, as suggested by Praud and Swinney
(2005). The Thomé et al. (1989) series (yellow) is for a pattern in which the tip
splitting is perfect (each splitting event leads to two long-lived structures) and it too
gives fair agreement with the DLA saturation profile in the frozen region. Lajeunesse
and Couder (2000) present a low flow rate, perfect tip-splitting pattern (blue) and
a high flow rate, imperfect tip-splitting pattern (red). The former pattern is very
similar to that of May and Maher (1989) – agreement is modest and limited to a
100 Dynamics of Fractal Growth
R
Rb
S
100 101
0
0.5
1
Fig. 4.13 Saturation profiles for the final patterns of May and Maher (1989) (+),
Thomé et al. (1989) (o), and Lajeunesse and Couder (2000) (their fig. 17(a) - (boxes))
and (fig 17.(b) - (diamonds)).
smaller internal area – whereas in the latter the agreement is more compelling and
over a larger range of radii.
One possible conclusion is that the condition for the saturation to approximately
vary according to DLA is simply that there are many tip splitting events. This appears
to be approximately true even in cases where surface tension acts to widen the fingers
as the interface slows.
This suggests that the growing region creates a systematic error in the measurement
of the fractal dimension. In figure 4.14, we show the results of the box counting
algorithm when applied overall (blue), to the shape that results from masking the
growing region (yellow) and to the shape that results from masking the frozen region
(red) for (a) the Praud and Swinney (2005) ∆P = 1.25 atm experiment (key: P3)
and (b) the constant flow C = 21 experiment (C2). In figure 4.14(a), we plot the box-
counting fractal dimension, −dlog(N)/dlog(ǫ), against box size ǫ, where ǫ is scaled
with ǫ0, the resolution of the image (given in the captions). The measurement of D is
approximately constant over the scales 100 < ǫ < 102, with D = 1.70± 0.02, and this
indicates self similarity on those scales for the overall shape (blue). When the growing
region is masked and the box-counting algorithm is applied to the frozen region alone
4.5 Discussion 101
Ov
Sc
Ac
ǫ
D
100 103
1
1.3
1.6
1.9
(a) P3; ǫ0 = 0.27mm.
Ov
Sc
Ac
ǫ
D
100 103
1
1.5
2
(b) C2; ǫ0 = 0.65mm.
Fig. 4.14 The box counting algorithm applied to different regions of two final VF
patterns. From the top series to the bottom: Ov = Overall, Sc = Screened (frozen)
and Ac = Active (growth).
(yellow), the measurement is much the same. This is a naive method, since masking
the growing region will have an effect on the fractal dimension of the frozen shape as
the tips of the fingers near the edge are trimmed, but it nonetheless indicates that for
the patterns presented in Praud and Swinney (2005), the box-counting measurement
of the fractal dimension is not significantly perturbed by the growing region. This is
on account of the large size of the frozen region and the large number of tip-splitting
events within it dominating the overall measurement. In contrast, in the final pattern
of experiment C2, the measurement of the fractal dimension monotonically decreases
indicating that the shape is not self-similar on those length scales and is therefore not
a fractal object. Thus, it is quite remarkable that – despite the lack of self similarity –
the measured values of Rtip approximately agree with the model’s predictions, which
were based on the fractal dimension for viscous fingering found by Praud and Swinney
(2005). Figure 4.14(b) also shows that the result of the box counting method is different
when considering the overall pattern rather than the frozen region alone for this less
complex pattern.
As a final note, the study of Chen (1989) was mentioned in the introduction, as
he also measured the fractal dimension of viscous fingering from the pattern’s mass
distribution (using the patterns’ density, D(R), where D(R) =
∫ R
0 2πRS(R)dR). We
have not further discussed the data of Chen (1989) since his density-radius plots
102 Dynamics of Fractal Growth
evolved in time at all radii, i.e. the frozen region in fact exhibited growth. In that
study, the plate separation was 6/10 that of the Praud and Swinney (2005) experiment
and the glass 11 times thinner, so it is possible that there was some degree of flexion
even though the flow rate is slower than in the Praud and Swinney (2005) experiments.
4.6 Conclusions
In this chapter, we have performed a series of new viscous fingering experiments, re-
analysed viscous fingering experiments in the literature and reproduced simulations of
Diffusion Limited Aggregation (DLA). We have explored a novel approach to analysing
the patterns that form, which reveals features that have not previously been reported.
This approach is to find the fraction of the area covered by the fingers in small annuli
at all radii, i.e. the distribution of saturation of the pattern.
Across all patterns, the evolution of the saturation profile confirms that growth is
concentrated at the extremities of the pattern, as has previously been noted. However,
the analysis of this chapter has found that there is a well-defined evolving radius, the
frozen radius, inside of which the pattern is static and outside of which – in the growth
region – the variation of saturation is approximately linear.
The viscous fingering experiments included a set of experiments involving a moder-
ately large constant flow rate, a set of high flow rate experiments involving a constant
pressure injection performed by Praud and Swinney (2005) and a set of moderate
flow rate "ramped" experiments, in which the flow rate linearly increased. A constant
flow rate experiment and a ramped experiment provided by Bischofberger et al. (2015)
were also analysed. In all cases, the flow rate was sufficient to induce the tip-splitting
instability.
The viscous fingering experiments featured similar saturation profiles to DLA: there
were seen to be a frozen region, in which the saturation decreased approximately ac-
cording to D = 1.70 (to within 5%), and a growth region with linear saturation
variation. However, in addition to these radii, in viscous fingering there is additionally
4.6 Conclusions 103
a fully-saturated inner region at radii smaller than the bases of the fingers, that is un-
seen in DLA. It was generally found that, following a period of development, the outer
and inner radii evolve in approximately constant ratio (1.40 with standard deviation
0.07), except in the two slower experiments with constant injection, in which this ratio
evolved throughout pattern development. This observation was used to predict that
the outer and inner radii of the pattern evolve proportional to t(1+α)/D, where α is the
time exponent of the injection flow rate. The data gave approximate agreement with
this prediction when the radii and time were scaled with R0 and t0. These constants
are the radius and time at which the linear stability theory prediction of the growth
rate of mode 2 first becomes positive.
In the DLA simulations, the saturation in the frozen region was seen to vary in
good agreement with the fractal dimension that has previously been reported for DLA,
D = 1.71 (Ossadnik, 1991). The saturation distribution also reveals that the outer and
inner radii of the growing region evolve in approximately constant ratio (1.54± 0.08)
throughout development.
When the inner and outer radii were scaled with radius of the finger bases and
shown varying with the total area of the pattern, the radii in viscous fingering were
shown to approximately agree with the values from DLA.
These findings lead to the overall conclusion that in unstable situations – where the
nonlinear development of viscous fingering dominates and the tip-splitting instability
arises – when injecting fluid with different flow-rate profiles the only difference in the
saturation distribution is the extent of the fully-saturated region.
Finally, it was noted that in large and complex patterns that feature self similarity,
measuring the fractal dimension of the overall pattern by the box counting method
gives a good estimate of the fractal dimension as can be used to predict the saturation
distribution. In contrast, in simple patterns that feature only a few generations of
tip splitting and lack self similarity, the box-counting method would not be useful for
predicting the saturation distribution.
104 Dynamics of Fractal Growth
4.7 Appendix
Table 4.4 The image processing MATLAB functions used and a brief description
thereof.
VideoReader Extracts each frame from a movie.
im2bw Converts frame to a binary image according to the appropriate
threshold value (found by trial-and-error).
bwareaopen Removes any objects smaller than a threshold number of pixels
(found by trial-and-error).
Ifill Fills any holes in the remaining large object.
bwboundaries Detects the edge of the pattern.
regionprops Determines pattern properties such as area and perimeter.
Chapter 5
Conclusions
This chapter is split into two sections. §5.1 lists the main contributions to research
that have been made from the work in this thesis and §5.2 lists some research ideas
that would draw on this thesis.
5.1 Main contributions
The stability of the interface between two immiscible fluids as one is injected into
the other from a point source in a two dimensional axisymmetric Hele-Shaw cell was
explored. Linear stability analysis of the interface identified that for a specific mode,
n, there is a critical viscosity ratio, Vc, required to ensure the absolute stability of that
azimuthal mode (§2.3.2). For viscosity ratios less adverse than the critical value, Vc,
the azimuthal mode is stable to all injection flow rates, Q. All modes m < n are also
stable under such flow conditions. However, higher order order modes may still be
unstable. We also showed that for a given viscosity ratio, there is a critical flow rate,
Q∗(t), such that all modes are stable. This is possible owing to the balance between
the gradual stretching of the interface, interfacial tension and viscous destabilisation.
Building on this theory, the case in which a second fluid is injected following
injection of the first fluid was also analysed. The two are injected in sequence leading
to the formation of a radially spreading annulus of the first injected fluid. The stability
106 Conclusions
of the annulus depends on the radial extent of the annulus relative to the outer radius
of the annulus. If the radial extent of the annulus is large then the growth rates of
instability on the leading and trailing interfaces can be considered in isolation (§2.4.1).
Generally however, there is coupling of the two interfaces. If the fluid in the annulus is
of intermediate viscosity as compared to the bounding fluids then for each azimuthal
mode, n, analogously to the single interface case there exists a critical viscosity ratio
of the bounding fluids, Vupper(n). For less adverse viscosity ratios the system is stable
to that particular azimuthal mode for all injection flow rates, Q, and at all times.
There exists another critical viscosity ratio, Vlower(n), below which that particular
azimuthal mode is not absolutely stable throughout the injection. For intermediate
viscosity ratios, Vlower(n) < V < Vupper(n), that azimuthal mode may be absolutely
stable initially but at some later time it becomes unstable on account of the coupling
of instability on the interfaces. In an analogous fashion to the single interface case, the
balance between destabilising and stabilising effects gives rise to an azimuthal mode
which features the largest growth rate in the system (§2.4.3) and this most unstable
mode gives rise to an injection flow rate, Q∗(t), such that if Q < Q∗(t) then the
interface is dynamically stable (§2.5). Given the viscosities of the bounding fluids, a
method was found for choosing the viscosity of the intermediate fluid that minimises
the total duration of the dynamically stable injection of the annulus.
In a rectilinear geometry, it may be of interest to determine the injection rate,
U(t), such that if a finite volume of fluid is to be injected in a finite time, then the
magnitude of any instability of the interface is minimised (§3.3.1). This particular
injection rate is found using variational calculus. It is shown that the optimal flow
rate strategy is a constant flow rate in the case of an injected fluid with constant
viscosity. This analysis is generalised to explore the effect of changes in viscosity of
the injected liquid as a function of time. In the case that the viscosity of the injected
fluid gradually increases with time, it is found that the optimal injection strategy
is adjusted so as to lead to the injection of more fluid later in the flow following
the increase in viscosity. It is demonstrated that this strategy can lead to substantial
5.1 Main contributions 107
suppression of instability compared to the case of a constant injection rate (§3.3.2). In
the case of an axisymmetric geometry, there are a series of discrete modes. If instead,
we assume there is a continuous series of modes, then we can follow a similar analysis
and find a time evolving injection rate, Q(t), that minimises the final amplitude of
instability (§3.4.1). If the injection rate is sufficiently low, the hypothetical growth
rate of instability associated with this continuous bound can be dynamically stabilised
by evolving the injection rate according to Q ∼ t−1/3. In contrast, if the injection is
relatively rapid then the optimal injection flow rate is linear in time, Q ∼ t, a limit
found by Dias et al. (2012). In the case of an axisymmetric displacement in which the
injected fluid gradually increases in viscosity, the optimal injection flow rate involves
injecting more of the fluid following the increase in viscosity than in the constant
injection case (§3.4.2).
In many cases of industrial interest, shear-thinning fluids are involved in displace-
ments in thin channels. We have generalised the analysis to consider the minimisation
of instability when a finite volume of shear-thinning fluid is injected in a finite time.
In the case of the displacement of a shear-thinning fluid in a rectilinear geometry, the
optimal injection flow rate involves injecting less fluid later in the flow than in the
Newtonian case (§3.5).
In the final part of the thesis, we have examined the nonlinear growth of the viscous
fingering instability, and in particular tried to assess how the fingering pattern grows
in time. To this end we have analysed the saturation, S, of the fingering pattern as
a function of the radius, R, where the saturation is the azimuthal average of the area
covered by the fingers of injected fluid. We have found that with the development of
a nonlinear fingering regime, there are three main regions: an inner fully saturated
region, an intermediate frozen region which is fixed in time, but whose leading edge
grows in time, and a growing nose region in which the saturation linearly decreases to
zero. The nose region advances radially with time, and in this region the saturation
at each point increases with time. The intermediate region has a saturation which
follows the law S ∼ R−0.30 for the experiments. These observations were found for
108 Conclusions
new experiments in which the injection flow rate was constant or linearly increased and
also for a previously published set of data for which the inlet pressure was held constant.
After a period of evolution, the radius of the intermediate region is approximately 2/3
the radius of the overall pattern in all cases. This allowed prediction of the evolution
of the radii. It is of interest that these results are largely independent of the injection
history. Unlike the linearly stability problem in which we have found optimal injection
rates with time, the nonlinear problem seems to evolve largely independently of the
specific injection rate. We have also performed numerical calculations of the Diffusion
Limited Aggregation (DLA) system and found very similar saturation distributions.
The extents of both the tip radius and the frozen zone radius relative to the fully-
saturated radius was shown to be approximately the same across all viscous fingering
experiments and DLA for a given pattern area.
5.2 Further Work
This thesis has explored the benefits of varying the flow rate during injection of a
continuous stream of polymer that is gradually undergoing a change in viscosity. How-
ever, towards the end of such an injection in the field, it may be desirable to clean up
the well by following the polymer injection with water. Thus, it would be of interest
to combine the principles developed in chapters 2 and 3, and explore the benefits of
modifying the injection flow rate in the case of the injection of a finite mass of poly-
mer laden fluid and a finite mass of water in series in a finite time. The case of the
intermediate polymer laden fluid gradually undergoing a change in viscosity could also
be explored. This problem would be complicated by the presence of both the global
and local modes of instability. Gorell and Homsy (1983) found the optimal spatial
distribution of the concentration of polymer to minimise fingering in the intermediate
fluid of a three fluid rectilinear system when the injection flow rate was constant. It
would be interesting to find this distribution in an axisymmetric geometry and, draw-
5.2 Further Work 109
ing on ideas from chapters 2 and 3, explore what potential there is to further suppress
instability by varying the injection flow rate in this case.
We have made use of the Hele-Shaw flow analogy to flow in porous media, having
acknowledged the limitations of this analogy. Whilst Chen and Yan (2017) have
explored the impact of the Dias et al. (2012) variational calculus framework towards
minimising instability in porous media, it was a computational treatment and it would
be fascinating to explore these principles in a series of porous media experiments.
Chen and Yan (2017) also found that the linearly increasing injection rate did not
lead to suppression of instability in the case of heterogeneity of the porous medium.
A non-uniform gap width in the cross section of a Hele-Shaw cell is analogous to a
heterogeneous porous medium. Viscous fingering in such a cell was explored by Woods
and Mingotti (2016), and it would be interesting to explore the potential of varying
the injection flow rate towards suppressing instability in this case.
The two geometries in which we have investigated flow in this thesis have been the
axisymmetric and the rectilinear. It would be of interest to generalise what has been
found to source-sink flows, such as the five-spot geometry.

References
Al-Housseiny, T. T. and Stone, H. A. (2013). Controlling viscous fingering in tapered
Hele-Shaw cells. Physics of Fluids, 25:092102.
Bataille, J. (1968). Stabilité d’un écoulement radial non miscible. Revue Institute
Petrole, 23:1349–1364.
Batista, C., Dias, E. O., and Miranda, J. A. (2016). Hamiltonian formulation towards
minimization of viscous fluid fingering. Physical Review E, 94(1):013109.
Beeson-Jones, T. H. and Woods, A. W. (2015). On the selection of viscosity to suppress
the Saffman Taylor instability in a radially spreading annulus. Journal of Fluid
Mechanics, 782:127–143.
Beeson-Jones, T. H. and Woods, A. W. (2017). Control of viscous instability by
variation of injection rate in a fluid with time-dependent rheology. Journal of Fluid
Mechanics, 829:214–235.
Bischofberger, I., Ramachandran, R., and Nagel, S. R. (2015). An island of stability
in a sea of fingers: emergent large-scale features of the viscous flow instability. Soft
Matter, 11(345):7428.
Brenner, E. A., Kessler, D., Levine, H., and Rappel, W. (1990). Selection of the
Viscous Finger in the 90° Geometry. Europhysics Letters, 13(2):161–166.
Callan-Jones, A. C., Joanny, J. F., and Prost, J. (2008). Viscous-fingering-like insta-
bility of cell fragments. Physical Review Letters, 100(25):258016.
Cardoso, S. S. S. and Woods, A. W. (1995). The formation of drops through viscous
instability. Journal of Fluid Mechanics, 289:351–378.
Carreau, P., DeKee, D., and Daroux, M. (1979). An analysis of the viscous behaviour
of polymeric solutions. Canadian Journal of Chemical Engineering, 57:135–140.
Casademunt, J. (2004). Viscous fingering as a paradigm of interfacial pattern forma-
tion: Recent results and new challenges. Chaos, 14(3):809–824.
Chen, C. Y. and Yan, P. Y. (2017). Radial flows in heterogeneous porous media with
a linear injection scheme. Computers and Fluids, 142:30–36.
Chen, J.-D. (1989). Growth of radial viscous fingers in a Hele-Shaw cell. Journal of
Fluid Mechanics, 201:223–242.
112 References
Chouke, R. L., van Meurs, P., and van der Poel, C. (1959). The instability of slow, im-
miscible, viscous liquid-liquid displacements in permeable media. Petroleum Trans-
actions, AIME, 216:188–194.
Combescot, R. and Ben Amar, M. (1991). Selection of Saffman-Taylor Fingers in the
Sector Geometry. Physical Review Letters, 67(4):453–456.
Couder, Y. (1988). Viscous Fingering in a Circular Geometry. In Stanley, H. E.
and Ostrowsky, N., editors, Random Fluctuations and Pattern Growth, pages 75–81.
Kluwer Academic Publishers, Dordrecht.
Cross, M. M. (1965). Rheology of non-Newtonian fluids: A new flow equation for
pseudoplastic systems. Journal of Colloid Science, 20(5):417–437.
Daccord, G., Nittmann, J., and Stanley, H. E. (1986). Radial Viscous Fingers and
Diffusion-Limited Aggregation: Fractal Dimension and Growth Sites. Physical Re-
view Letters, 56(4):336–339.
Davidovitch, B., Levermann, A., and Procaccia, I. (2000). Convergent calculation of
the asymptotic dimension of diffusion limited aggregates: Scaling and renormaliza-
tion of small clusters. Phys. Rev. E, 62:R5919–R5922.
Dias, E. O., Alvarez-Lacalle, E., Carvalho, M., and Miranda, J. (2012). Minimization
of Viscous Fluid Fingering: A Variational Scheme for Optimal Flow Rates. Physical
Review Letters, 109(14):144502.
Dias, E. O. and Miranda, J. (2013a). Taper-induced control of viscous fingering in
variable-gap hele-shaw flows. Physical Review E, 87:053015.
Dias, E. O. and Miranda, J. A. (2010). Control of radial fingering patterns: A weakly
nonlinear approach. Physical Review E, 81(1):016312.
Dias, E. O. and Miranda, J. A. (2013b). Wavelength selection in Hele-Shaw flows: A
maximum-amplitude criterion. Physical Review E, 88(1):013016.
Dias, E. O., Parisio, F., and Miranda, J. A. (2010). Suppression of viscous fluid
fingering: A piecewise-constant injection process. Physical Review E, 82(6):067301.
Fontana, J. V., Dias, E. O., and Miranda, J. A. (2014). Controlling and minimizing
fingering instabilities in non-Newtonian fluids. Physical Review E, 89(1):013016.
Foroutan-pour, K., Dutilleul, P., and Smith, D. (1999). Advances in the implemen-
tation of the box-counting method of fractal dimension estimation. Applied Mathe-
matics and Computation, 105(2-3):195–210.
Ghannam, M. and Esmail, M. (1998). Rheological properties of aqueous polyacry-
lamide solutions. Journal of Applied Polymer Science, 69(8):1587–1597.
Gin, C. and Daripa, P. (2015). Stability results for multi-layer radial hele-shaw and
porous media flows. Physics of Fluids, 27:012101.
References 113
Gorell, S. and Homsy, G. M. (1983). A theory of optimal policy of oil recovery by
secondary displacement processes. SIAM Journal on Applied Mathematics, 43:79–
98.
Gun, W. J. and Routh, A. F. (2013). Microcapsule flow behaviour in porous media.
Chemical Engineering Science, 102:309–314.
Han, M.and Shi, L., Ye, M., and Ma, J. (1995). Crosslinking reaction of polyacrylamide
with chromium(Ill). Polymer Bulletin, 35:109–113.
Hele-Shaw, H. (1898). On the motion of a viscous fluid between two parallel plates.
Nature, 58:34–36.
Hill, S. (1952). Channeling in packed columns. Chemical Engineering Science,
1(6):247–253.
Homsy, G. (1987). Viscous Fingering in Porous Media. Annual Review of Fluid
Mechanics, 19:271–311.
Huang, Y. S. and Chen, C. Y. (2015). A numerical study on radial Hele-Shaw
flow: influence of fluid miscibility and injection scheme. Computational Mechan-
ics, 55(2):407–420.
Jha, B., Cueto-Felgueroso, L., and Juanes, R. (2011). Fluid Mixing from Viscous
Fingering. Physical Review Letters, 106(19):194502.
Kamal, M. S., Sultan, A. S., Al-Mubaiyedh, U. A., and Hussein, I. A. (2015). Review
on Polymer Flooding: Rheology, Adsorption, Stability, and Field Applications of
Various Polymer Systems. Polymer Reviews, 55(3):491–530.
Kim, H., Funada, T., Joseph, D. D., and Homsy, G. M. (2009). Viscous potential flow
analysis of radial fingering in a Hele-Shaw cell. Physics of Fluids, 21(7):074106.
Korteweg, D. (1901). Sur La Forme Que Prennent Les Équations Du Mouvements Des
Fluides Si L’on Tient Compte Des Forces Capillaires Causées Par Des Variations
De Densité Considŕables Mais Coninues Et Sur La Thórie De La Capillarité Dans
L’hypoths`e D’une Variation Continue De La Densité. Archives Néerlandaises des
Sciences Exactes et Naturelles, 6.
Kuijpers, K. R., De Martín, L., and Van Ommen, J. R. (2014). Optimizing off-lattice
Diffusion-Limited Aggregation. Computer Physics Communications, 185(3):841–
846.
Lajeunesse, E. and Couder, Y. (2000). On the tip-splitting instability of viscous fingers.
Journal of Fluid Mechanics, 419:125–149.
Lake, L. (1989). Enhanced Oil Recovery. Prentice Hall, Englewood Cliffs, NJ.
Lee, K. E., Morad, N., Teng, T. T., and Poh, B. T. (2012). Kinetics and In Situ
Rheological Behavior of Acrylamide Redox Polymerization. Journal of Dispersion
Science and Technology, 33(3):387–395.
114 References
Li, S., Lowengrub, J., Fontana, J., and Palffy-Muhoray, P. (2009). Control of Viscous
Fingering Patterns in a Radial Hele-Shaw Cell. Physical Review Letters, 102:174501.
Li, S., Lowengrub, J. S., and Leo, P. H. (2007). A rescaling scheme with application
to the long-time simulation of viscous fingering in a Hele-Shaw cell. Journal of
Computational Physics, 225(1):554–567.
Lubkin, S. and Murray, J. (1995). A mechanism for early branching in lung morpho-
genesis. Journal of Mathematical Biology, 34(1):77–94.
Makadia, H. K. and Siegel, S. J. (2011). Poly Lactic-co-Glycolic Acid (PLGA) as
biodegradable controlled drug delivery carrier. Polymers, 3(3):1377–1397.
Martyushev, L. M., Birzina, A. I., Konovalov, M., and Sergeev, A. (2015). Mor-
phological stability of an interface between two non-Newtonian fluids moving in a
Hele-Shaw cell. Physical Review E, 80:066306.
Mathiesen, J., Procaccia, I., Swinney, H. L., and Thrasher, M. (2006). The universality
class of diffusion-limited aggregation and viscous fingering. Europhysics Letters
(EPL), 76(2):257–263.
Maxworthy, T. (1989). Experimental study of interface instability in a Hele-Shaw cell.
Physical Review A, 39(11):5863–5866.
May, S. E. and Maher, J. V. (1989). Fractal dimension of radial fingering patterns.
Physical Review A, 40(3):1723–1726.
McBirney, A. R. (1984). Igneous Petrology. Freeman, Cooper and Co., San Francisco.
McCloud, K. V. and Maher, J. V. (1995). Experimental perturbations to Saffman-
Taylor flow. Physics Reports, 260(3):139–185.
Meakin, P. (1998). Fractals, Scaling and Growth Far from Equilibrium. Cambridge
University Press, Cambridge.
Miranda, J. A. and Widom, M. (1998). Radial fingering in a hele-shaw cell: a weakly
nonlinear analysis. Physica D: Nonlinear Phenomena, 120:315–328.
Mora, S. and Manna, M. (2009). Saffman-Taylor instability for generalized Newtonian
fluids. Physical Review E, 80(1):016308.
Mungan, N. (1971). Improved waterflooding through mobility control. The Canadian
Journal of Chemical Engineering, 49(1):32–37.
Nayfeh, A. H. (1972). Stability of liquid interfaces in porous media. Physics of Fluids,
15:1751–1754.
Ossadnik, P. (1991). Multiscaling analysis of large-scale off-lattice dla. Physica A:
Statistical Mechanics and its Applications, 176(3):454 – 462.
Paraskeva, C. A., Charalambous, P. C., Stokka, L. E., Klepetsanis, P. G., Koutsoukos,
P. G., Read, P., Ostvold, T., and Payatakes, A. C. (2000). Sandbed Consolidation
with Mineral Precipitation. Journal of Colloid and Interface Science, 232:326–339.
References 115
Parisio, F., Moraes, F., Miranda, J. A., and Widom, M. (2001). Saffman-Taylor
problem on a sphere. Physical Review E, 63:036307.
Park, C. and Homsy, G. (1984). Two-phase displacement in Hele Shaw cells : theory.
Journal of Fluid Mechanics, 139:291–308.
Paterson, L. (1981). Radial fingering in a Hele Shaw cell. Journal of Fluid Mechanics,
113:513–539.
Paterson, L. (1985). Fingering with miscible fluids in a Hele Shaw cell. Physics of
Fluids, 28(1):26–30.
Perkins, T. and Johnston, O. (1969). A Study of Immiscible Fingering in Linear
Models. Society of Petroleum Engineers Journal, 9(1):39–46.
Pihler-Puzović, D., Illien, P., Heil, M., and Juel, A. (2012). Suppression of Complex
Fingerlike Patterns at the Interface between Air and a Viscous Fluid by Elastic
Membranes. Physical Review Letters, 108(7):074502.
Pojman, J. A., Whitmore, C., Liria, M., Liveri, T., Lombardo, R., Marszalek, J.,
Parker, R., and Zoltowski, B. (2006). Evidence for the Existence of an Effective
Interfacial Tension between Miscible Fluids : Isobutyric Acid-Water and 1-Butanol-
Water in a Spinning-Drop Tensiometer. Langmuir, 22(6):2569–2577.
Praud, O. and Swinney, H. L. (2005). Fractal dimension and unscreened angles mea-
sured for radial viscous fingering. Physical Review E, 72(1):011406.
Saffman, P. G. and Taylor, G. I. (1958). The penetration of a fluid into a porous
medium or hele-shaw cell containing a more viscous liquid. Journal of Fluid Me-
chanics, 245:312–329.
Sorbie, K. (1991). Polymer-Improved Oil Recovery. CRC Press, Boca Raton, Florida.
Talaghat, M. R., Esmaeilzadeh, F., and Mowla, D. (2009). Sand production control by
chemical consolidation. Journal of Petroleum Science and Engineering, 67(1-2):34–
40.
Tanveer, S. (1987). Analytic theory for the linear stability of the Saffman–Taylor
finger. Physics of Fluids, 30(8):2318–2329.
Thomé, H., Rabaud, M., Hakim, V., and Couder, Y. (1989). The Saffman–Taylor
instability: From the linear to the circular geometry. Physics of Fluids A, 1(2):224–
240.
Tran-viet, A., Routh, A., and Woods, A. (2014). Control of the Permeability of a
Porous Media Using a Thermally Sensitive Polymer. American Institute of Chemical
Engineers, 60(3):1193–1201.
Wilson, S. D. R. (1990). The Taylor-Saffman problem for a non-Newtonian liquid.
Journal of Fluid Mechanics, 220:413–425.
Witten, T. and Sander, L. (1981). Diffusion-Limited Aggregation, a Kinetic Critical
Phenomenon. Physical Review Letters, 47(19):1400–1403.
116 References
Woods, A. W. (2015). Flow in Porous Media. Cambridge University Press, Cambridge.
Woods, A. W. and Mingotti, N. (2016). Topographic viscous fingering : fluid – fluid
displacement in a channel of non-uniform gap width. Philosophical Transactions of
the Royal Society A, 374:20150427.
Zheng, Z., Kim, H., and Stone, H. A. (2015). Controlling Viscous Fingering Using
Time-Dependent Strategies. Physical Review Letters, 115(17):174501.