# On the fractal dimension of non-Newtonian Hele-Shaw flow subject to Saffman-
Taylor instability

J. Adriazola  Department of Mathematics, Southern Methodist University,
Dallas, Texas, 75205  B. Gu  Department of Mathematical Sciences, Worcester
Polytechnic Institute, Worcester, Massachusetts, 01609  L.J. Cummings and L.
Kondic  Department of Mathematical Sciences, New Jersey Institute of
Technology, Newark, New Jersey 07102, USA

###### Abstract

We introduce a discrete numerical method based on the diffusion-limited
aggregation (DLA) approach to simulate two-fluid Hele-Shaw flow subject to the
Saffman-Taylor interfacial instability, in the case where the displaced fluid
is non-Newtonian. Focusing on fluids for which the most relevant non-Newtonian
aspect of the thin-gap flow is shear-thinning, we introduce a history-
dependent aspect into the algorithm, modeling shear-rate-dependent fluid
viscosity. The main finding is that the morphology of the emerging patterns,
characterized by the fractal dimension, is modified in a nontrivial manner by
the shear-thinning nature of the displaced fluid. In particular, we
consistently find that shear-thinning leads to the formation of patterns
characterized by a smaller fractal dimension, compared to the corresponding
Newtonian fluid.

The classical Hele-Shaw (H-S) experiment [1] is simple: inject a fluid into an
immiscible, more viscous fluid confined between two closely-spaced plates.
Such a setup is known to be subject to the Saffman-Taylor (S-T) instability
[2], resulting in complex, fractal-like pattern formation. Extensive reviews
of experiments [3], and of various aspects of theoretical and computational
findings [4, 5, 6, 7] are available. Figure 1(a) shows a typical example of
experimental results obtained in NJIT’s Capstone Laboratory [8] – many such
examples can be found in the literature; see [3] and the references therein.

A key reason for the interest in H-S flow and related instabilities, in
addition to easy visualization, is that its mathematical description is
similar to Darcy’s formulation modeling porous media flow, which is highly
relevant to applications that vary from secondary oil recovery to injection
molding [9] and many other natural and man-made setups. Furthermore, there is
a close connection between the S-T problem and the Mullins-Sekerka instability
relevant to the propagation of a solidification front in an undercooled liquid
[10]. Therefore, there are multiple reasons to work towards understanding the
nature of pattern formation and emergent interface morphologies.

While a significant body of experimental work has been carried out [11, 12,
13, 14, 15, 16] (some with non-Newtonian fluids), it is difficult to reach
precise understanding of the morphology of the emerging patterns, for either
Newtonian or non-Newtonian fluids. Such morphology is often characterized by
the fractal dimension, $D_{\rm f}$, for which values in the range 1.2 - 2.0
have been reported; see, e.g., [17, 15, 3] for discussion. Most of the
(relatively) recent experiments seem to converge to a value close to 1.8 [3],
but the differences between results are too large to attempt to identify the
influence of non-Newtonian effects in the displaced fluid. This influence is
crucial to numerous applications; however, due to inherent limitations of the
experimental setups, it is difficult to quantify such differences based on
experiments alone [14, 3].

From the theoretical perspective, non-Newtonian free surface flows are
nontrivial to model due to often complicated rheological properties. However,
the H-S geometry once again comes to the rescue since, in the thin gap limit,
for a significant class of non-Newtonian fluids, the most important aspects of
their rheology can be reduced to a shear-thinning model involving pressure-
gradient dependent viscosity [18, 19, 16, 20]. Such a simplification was
carried out systematically some years back [21] showing that, even if elastic
effects may be important in general, the thin gap limit reduces the problem to
the Darcy-type law ${\bf u}=-b^{2}\nabla p/[12\mu(|\nabla p|^{2})]$, where
${\bf u}$ is the gap (of thickness $b$) averaged fluid velocity, $p$ is the
pressure, $\nabla=(\partial_{x},~{}\partial_{y})$, with $(x,y)$ the in-plane
coordinates, and $\mu(|\nabla p|^{2})$ is the pressure-gradient dependent
viscosity. We point out that $\mu(|\nabla p|^{2})$ is obtained by gap
averaging shear-rate dependent viscosity $\bar{\mu}(|{\bf u}_{z}|^{2})$; see
[18] for details of how to translate between the two viscosity
representations. Such modified Darcy formulations can be obtained by
systematic asymptotic expansions applicable to a wide class of non-Newtonian
fluids, even when the Weissenberg number quantifying the relevance of elastic
effects is $O(1)$. When combined with the incompressibility condition,
$\nabla\cdot{\bf u}=0$, a nonlinear elliptic problem for the fluid pressure is
obtained, which is significantly easier to deal with than the original non-
Newtonian formulation involving elastic effects. Despite this simplification,
it is still challenging and computationally intensive to solve for the
complicated interfaces that ultimately develop due to the S-T instability; the
problem is of free-boundary type involving in-plane curvature through the
Young-Laplace boundary condition $\bar{p}=-\gamma\kappa$, where $\bar{p}$ is a
suitably renormalized pressure, and $\gamma$ the surface tension; see e.g. [7]
for details. Thus, since carrying out large simulations is difficult,
quantifying the emerging pattern morphology and its dependence on non-
Newtonian fluid behavior remains elusive.

An alternative approach is modeling using discrete methods, particularly
algorithms based on diffusion-limited aggregation (DLA) [22]. DLA is a version
of Monte-Carlo simulation, which uses a growing aggregate of particles to
simulate solutions of Laplace’s equation with free boundaries. In the context
of H-S flow, DLA-type algorithms have been modified to include surface tension
effects [23] and used extensively to discuss the pattern formation that
emerges in the S-T instability of displaced Newtonian fluids [17]: it was
found that surface tension modifies the fractal dimension $D_{\rm f}$ of the
emerging patterns from the “pure” DLA value of 1.67 to larger values. Exact
values of $D_{\rm f}$ are often unclear, however, and there is a significant
body of work discussing the asymptotic value of $D_{\rm f}$ in the limit of
large aggregate size, as well as the fractal structure of the emerging
patterns [24, 25, 26]. In any case, to the best of our knowledge, whether and
to what extent the non-Newtonian character of the displaced fluid influences
the value of $D_{\rm f}$, either in experiments or in simulations, is
uncertain.

(a)

(b)

Figure 1: Examples of (a) Newtonian Hele-Shaw flow and (b) shear-thinning
flow.

In this Letter, we formulate a DLA-type method to model the Hele-Shaw flow of
a non-Newtonian fluid. The motivation for this study comes in part from
experiments that we briefly report in what follows. The main modeling idea is
to incorporate the fluid viscosity’s shear rate dependence via a history-
dependent pattern growth. We then discuss how such modification influences the
morphology of the emerging patterns, focusing for brevity on the fractal
dimension of the emerging patterns only.

#### Experiments

Figure 1 shows two examples of experimental patterns obtained by injecting, in
a controllable manner, water into an H-S cell initially containing more
viscous fluid, that is either Newtonian (glycerol) (a) or non-Newtonian
(polyethylene oxide, PEO) (b). The H-S cell consists of two plexiglass plates
of approximate thickness of 1/2 inch, placed typically at a distance 100-800
$\mu$m apart, with a 2mm wide hole in the top plate. A syringe is connected to
the hole via a plastic tube and water is injected by controlled applied
pressure. Both applied pressure and the spacing between the plates were varied
in the experiments; the video recordings are available [8]. While the
experiments were conducted carefully and consistent results obtained, we note
that they are of limited scope and intended only to illustrate visual clues
suggesting that the non-Newtonian response of PEO modifies the morphology of
the evolving patterns; the interested reader is directed to the extensive
review [3] and numerous references therein showing consistent findings. We
have computed the fractal dimensions, $D_{\rm f}$ of the emerging patterns
(using methods discussed later in the text), finding systematically lower
values for the non-Newtonian fluid. However, since computing $D_{\rm f}$
involves consideration of different scales that should span orders of
magnitude, it is necessary to consider patterns that are much larger than
those shown in Fig. 1. For this purpose, we resort to simulations.

#### Simulation Methods.

The DLA algorithm that provides the basis for this work is relatively simple,
particularly the on-lattice version we consider. The Brownian motion required
by the DLA model is simulated by a random walk along the four possible
directions of the grid. One starts with a seed particle at the center of a
lattice, generates a particle (‘walker’) far away, allows the walker to
perform a random walk on the lattice until it arrives at a lattice site
adjacent to the seed, where it sticks (with unit probability), transforming
this particular lattice site from free to occupied (in the algorithm that we
use this transformation is permanent; see [27] for more details and references
to alternative formulations).

To simulate the S-T instability in the H-S setup, the basic DLA algorithm was
modified by Vicsek [23] via the sticking probability. The modification
promotes sticking in neighborhoods that contain a large number of occupied
sites, since such regions, on average, are characterized by a smaller
curvature. Within this model, the sticking probability is defined by

$p_{\rm Newt}=A\left({n_{\cap}\over n_{\rm{total}}}-n_{0}\right)+B,$ (1)

where $n_{\cap}$ denotes the number of occupied sites in some neighborhood of
linear dimension $l$ (measured in the number of grid points) centered at the
considered sticking location, and $n_{\rm total}=l^{2}$ (Vicsek [23] used
$l=17$ and $n_{0}=(l-1)/(2l)$). The original DLA model is recovered by setting
$A=0$ and $B=1$; choosing $B\in(0,1)$ and increasing the value of $A$ models a
stronger surface tension since this encourages walkers to stick in areas of
high local particle density, hence, smaller local curvature.

We wish to simulate both Newtonian and non-Newtonian flows, which requires
changes to the algorithm (discussed in some detail in [27]; here we provide
only an outline). In our Newtonian simulations, we improve Vicsek’s algorithm,
modifying it to avoid the formation of ‘holes’ (unoccupied sites surrounded by
occupied ones), and by promoting azimuthal symmetry and minimizing grid
anisotropy.

Building on this improved Newtonian algorithm, we further adapt it to model
interface growth for a non-Newtonian (and in particular, shear-thinning)
displaced outer fluid. To this end, we need to modify the sticking probability
so as to model shear-thinning behavior. This can be done by realizing that,
instead of a shear-rate-dependent viscosity as discussed in the introduction,
we may instead think of a viscosity that depends on the velocity magnitude of
the displaced fluid, since the two quantities (shear rate and velocity
magnitude) are related (fluid velocity is obtained by gap averaging the fluid
shear rate). In the context of DLA-based simulations, we can then model shear
thinning behavior by making the sticking probability velocity-magnitude
dependent. This is the basic idea underpinning the selected approach.

In physical experiments, the dependence of viscosity on shear-rate may vary
widely. Here, we do not attempt to model any particular fluid; instead (with
the above discussion in mind) we assume a simple velocity-dependent non-
Newtonian correction to the sticking probability of the form

$p_{\rm NN}=CV^{s}.$ (2)

Here, $C$ and $s$ are two free parameters, chosen in such a way that $p_{\rm
NN}$ reaches values comparable to $p_{\rm Newt}$ (one could think of these
parameters as modeling the shear-thinning properties of the displaced fluid);
while $V$ is the velocity magnitude (speed), inversely proportional to the
“age” of a particle that has since become part of the aggregate (see [27] for
details). This means that younger (more recently attached) particles are more
likely to attract walking particles. The local age of the aggregate has two
components: a “horizontal” age $T_{H}$ and a “vertical” age $T_{V}$, which are
determined pointwise on each lattice site. To illustrate the algorithm,
suppose the aggregate comprises $k-1$ particles. Depending on how the $k$th
particle attaches, its age is determined differently. If the particle $k$
attaches horizontally, then its vertical age is inherited directly from the
occupied lattice site responsible for the attachment (neighbor), and the
horizontal age is defined by $T_{H}=k-T_{H}^{\rm neighbor}$. For a vertical
attachment, the two ages are calculated analogously. The velocity is then
given by the vector $\mathbf{v}=\left[T_{H}^{-1},T_{V}^{-1}\right]$, and the
speed by $V=||\mathbf{v}||$. There are additional considerations if a walker
attaches to more than one neighbor at once; see [27] for modifications of our
algorithm in such cases. The probability of attachment that we use for non-
Newtonian simulations includes Newtonian and shear-thinning effects through
the following expression,

$p_{\rm total}=\begin{cases}0,&\mathrm{if}\ p_{\rm Newt}+p_{\rm NN}<0,\\\
1,&\mathrm{if}\ p_{\rm Newt}+p_{\rm NN}>1,\\\ p_{\rm Newt}+p_{\rm
NN},&\mathrm{otherwise},\end{cases}$ (3)

with $p_{\rm NN}$ as defined in (2).

(a)

(b)

Figure 2: (a) Newtonian and (b) non-Newtonian simulations of pattern growth.
Both simulations use $A=0.4$ and $B=0.4$, while the non-Newtonian simulation
also uses $C=0.3$ and $s=0.2.$ Both plots are on $1600\times 1600$ grids; the
color palette shows the number of attached walkers.

#### Results

Figure 2 shows two examples of our Newtonian (a) and shear-thinning (b)
simulations. Visually, these figures suggest the formation of skinnier
patterns and possibly different morphology for the shear-thinning case;
however, the differences are not simple to quantify, and more importantly,
they are dependent on the pattern size, as observed in previous experiments
[3] and simulations [25]. Our finding, consistent with Newtonian simulations
[25], is that to obtain convergence of the fractal measures quantifying the
emerging patterns we require millions of walkers (the simulation examples in
Fig. 2 show $\sim 10^{5}$ walkers, two-to-three orders of magnitude less than
required to obtain converged results for $D_{\rm f}$). We therefore proceed to
discuss the results of much larger simulations, and characterize the
morphologies of the resulting aggregates by computing $D_{\rm f}$ using box-
counting and correlation algorithms; see [27] for details. We note that
additional measures, such as radius of gyration, were considered but did not
provide additional insight and are not reported here for brevity.

Figure 3 shows $D_{\rm f}$ as a function of the number of walkers (measured in
millions) as the parameters $A,~{}B$ (entering the definition of $p_{\rm
Newt}$) are modified. First, Figure 3(a) shows the results for $A=0$, without
surface tension effects. We observe that large simulations are indeed needed
to reach both convergence and an agreement between the two calculations; for
less than $\sim$10 million walkers, the two measures produce differing
results. For sufficiently many walkers, we find the value $D_{\rm f}\approx
1.67$, in agreement with values reported in the literature for on-lattice DLA
simulations [4, 25, 28]. We have confirmed that our results are realization-
independent; repeating a simulation with different random seed generators
produces consistent $D_{\rm f}$ values (for more than 10 million walkers).
Regarding the influence of the non-Newtonian correction, recall that within
the DLA algorithm the probability of sticking is unity, and therefore, the
non-Newtonian contribution provided by $p_{\rm NN}$ is irrelevant in this case
(see (3)). Assuming that (Newtonian) DLA represents well the physical S-T
problem with vanishing surface tension, the (perhaps not so obvious)
conclusion is that $D_{\rm f}$ is not influenced by the non-Newtonian nature
of the displaced fluid in this case.

(a)

(b)

(c)

(d)

Figure 3: Fractal dimensions using Box Counting (“$\circ$”) and Correlation
Dimension (“$\times$”), for Newtonian (red) and non-Newtonian (blue)
simulations for 1 to 25 million attached walkers, as the surface tension
parameter $A$ is varied ($A=0$ (a), 1 (b), 2 (c), 4 (d)). For $A=0$ (DLA) we
use $B=1$, and for $A\neq 0$, $B=0.4$; also for the non-Newtonian cases,
$C=0.9,~{}s=0.07$.

We now consider the effect of changing the surface tension parameter, $A$.
Figure 3(b - d) shows the results for $A=1,~{}2,~{}4$. As $A$ increases,
$D_{\rm f}$ increases for both Newtonian and non-Newtonian simulations, but
differently. We identify a gap between the corresponding Newtonian and non-
Newtonian fractal dimensions, $\Delta D_{\rm f}=D_{\rm f}^{\rm Newt}-D_{\rm
f}^{\rm nNewt}$ and find that the size of $\Delta D_{\rm f}$ depends non-
monotonously on the value of $A$, with the maximum value observed for $A=2$,
see Fig. 3(c). For yet larger values of $A$, $\Delta D_{\rm f}$ is found to
disappear completely (figures not included for brevity).

The results of Fig. 3 lead naturally to the following question: for a given
surface tension parameter $A$, is $\Delta D_{\rm f}$ non-negative for all
choices of shear-thinning parameters $(C,s)$? To investigate this question, we
carried out additional simulations with several different $(C,s)$ pairs,
simulating fluids with different shear-thinning properties. The results, given
in Table III in [27], suggest that the answer to this question is affirmative:
while the size of $\Delta D_{\rm f}$ depends on $(C,s)$ and on $A$, our
numerical experiments always find $\Delta D_{\rm f}\geq 0$.

#### Conclusion

In this work, we discuss the fractal dimension, $D_{\rm f}$, of the patterns
that form due to the Saffman-Taylor instability of Hele-Shaw flow. We find
that the shear thinning property of the more viscous fluid leads to an overall
decrease of $D_{\rm f}$, compared to its Newtonian counterpart. This result is
obtained from large-scale simulations of our proposed model, which modifies
the well-known approach based on the DLA algorithm for simulating unstable
Hele-Shaw two-fluid flow. Such a result would be difficult to obtain either
experimentally or via continuous PDE-based simulations, since extremely large
patterns are needed to produce accurate and converged results.

The manner in which $D_{\rm f}$ changes due to non-Newtonian behavior is
nontrivial. First, for vanishing surface tension, non-Newtonian behavior is
not relevant since our model reduces to classical DLA simulations. For very
large surface tension, the difference in $D_{\rm f}$ between Newtonian and
non-Newtonian simulations becomes small, which may be expected in hindsight
since the local curvature effects are dominant. For intermediate surface
tension values, however, we find consistently smaller values of $D_{\rm f}$
for non-Newtonian fluids, compared to Newtonian ones. We hope that our
findings will encourage new research, both experimental and theoretical, that
will make progress in identifying the general connection between fluid
rheological properties and fractal dimension of the emerging patterns.

Supplementary Material:
On the fractal dimension of non-Newtonian Hele-Shaw flow subject to Saffman-
Taylor instability

##  I Monte Carlo simulations: modified Vicsek algorithm

The main text provides the basic description of the Monte Carlo-based model
that we use for simulating aggregate growth. The algorithm we use can simulate
both Newtonian and shear-thinning non-Newtonian Hele-Shaw flow. We will first
discuss the Newtonian aspects at the core of this algorithm and delay a
discussion of the non-Newtonian modifications until Sections II and III.

The algorithm is based on an irreversible growth model known as diffusion
limited aggregation (DLA) [22]; see also [4] for a concise review. Briefly,
DLA is a random walk on a lattice that leads to the formation of an aggregate
from a seed particle. To align this with the Hele-Shaw geometry, the seed
particle is placed at the origin of a planar rectangular grid. Subsequent
particles are initialized one at a time at infinity (from a relatively large
distance in the practical implementation) and randomly walk until they come in
contact with the seeded aggregate, at which point they stick irreversibly.
This process is then repeated until an aggregate of desired size is formed.

To simulate the viscous stress forces present in Newtonian flow, Vicsek uses
this DLA approach but with a probability of sticking dependent on the local
density of the aggregate [23]. This sticking probability reflects the Laplace-
Young boundary condition, relating the jump in pressure to the local curvature
of the interface, in the DLA framework. When in the vicinity of the aggregate,
the random walker can detect the aggregate’s local density by counting the
number of occupied lattice sites within a suitably defined neighborhood, taken
to be a square of side length $l\in\mathbb{Z}^{+}$ and odd, centered at the
walker’s lattice site (the length is measured in units of grid points).

Vicsek’s sticking probability is given by the three-parameter expression

$p_{\rm
Vicsek}(n_{\cap})=A\left(\frac{n_{\cap}}{n_{\rm{total}}}-n_{0}\right)+B,$ (S1)

where $n_{\cap}$ denotes the number of occupied sites in the walker’s
neighborhood. The three free parameters are $A\in\mathbb{R}^{+},\ B\in[0,1],$
and the odd positive integer $l$, which enters via the $l-$dependent
parameters $n_{\rm total}=l^{2}$ and $n_{0}=(l-1)/(2l)$. To ensure that this
sticking probability is normalized, we restrict $p_{\rm Vicsek}(n_{\cap})$ to
the interval $[0,1].$ Vicsek showed that this heuristic approach works well in
practice to simulate Newtonian Hele-Shaw flow with suitably chosen parameters
$(A,B,l)$ [23].

We have made two modifications to Vicsek’s approach, which are still within
the domain of Newtonian flow. First, the aggregate’s topology may become
multiply-connected, i.e., holes may form. We implement a simple rule to avoid
this physically unrealistic (in the context of viscous flow) scenario. We use
the following criteria: if a free walker is found adjacent to the aggregate,
we check the eight lattice sites that enclose it. If an adjacent site is
occupied and its counter-clockwise neighbor site is not, or vice-versa, we
call this a flip. If the number of flips is greater than 3 after checking the
entire adjacent perimeter of the particle, then a potential hole has formed;
the walker does not stick; it moves back to its immediately previous location,
and we perform another random walk until the particle is adjacent to the
aggregate once again.

The second modification we make is to use a different metric to define the
local neighborhood of the walker. Instead of the rectangular metric used by
Vicsek, we use a circular one. In other words, the walker’s neighborhood is
now a circle instead of a square. This further requires us to reinterpret the
parameters $n_{0}$ and $n_{\rm total}$ defining Eq. (S1). In the original
formulation, Vicsek motivates the choice of these parameters by the fact that
$n_{0}=(l-1)/(2l)$ corresponds to the number of lattice points contained by a
rectangle of size $l\times(l-1)/2$ divided by the total number of lattice
points $n_{\rm total}=l^{2}$. In other words, $n_{0}$ is the number of lattice
points that would be occupied if the walker were in contact with a flat
interface. As an example, if $l=5$, then the neighborhood is a $5\times 5$
square with $n_{\rm total}=25$ and $n_{0}=4/10=10/25.$

In our modified approach, instead of a rectangle, we use a circle of diameter
$l$ centered at the walker position and once again count the number of
occupied sites within this circle. Since this circle is defined on top of the
rectangular geometry specified by the grid, finding a general formula that
expresses $n_{0}$ in terms of $l$ is difficult. For example, for $l=5$,
$n_{\rm total}=13$ and $n_{0}=4/13$. In the case of $l=17,$ a value we use in
all simulations throughout our work, we find by direct counting that
$n_{0}=89/{n_{\rm total}}$, with $n_{\rm total}=193$ here and in Eq. (1) in
the main text. Figure S1 illustrates the typical effects of using the modified
versus the original algorithm. We find that the modified approaches encourage
azimuthal symmetry of the emerging pattern for larger values of the surface
tension parameter $A$.

Figure S1: A comparison of modifying the neighborhood used in Newtonian
simulations. All four simulations have $B=0.4$ and $l=17.$ The top two panels
have $A=4$ while the bottom two panels have $A=16.$ The left panels use a
square neighborhood, while the right use a circular one. All plots are on
$1600\times 1600$ grids.

##  II Discussion of simulation parameters

Equations (2), (3) and the accompanying description in the main text explain
how we further modify Vicsek’s algorithm to account for non-Newtonian effects
due to a shear-thinning displaced fluid. We choose the parameters $C$ and $s$
in Eq. (2) (main text) in a manner that leads to an appreciable change in the
sticking probability. This is done to mimic typical shear-thinning rheology:
larger velocity magnitude of the interface leads to larger sticking
probability, modeling smaller viscosity. Different non-Newtonian fluids have
different responses to shear, thus, we can think of different values of these
parameters as modeling different fluids.

To aid in deciding appropriate parameter values to use in our implementation,
we plot the maximum value of the non-Newtonian contribution to the sticking
probability, $p_{\rm NN}$, recorded over successively larger periods (labeled
by integer $k$) of newly attached walkers (the length of each period is the
largest integer of $1.3^{k+1}-1.3^{k}$, with the data points in Fig. S2
corresponding to $k\in[16,45]$). The choice of $1.3$ and the range of $k$ is
arbitrary and was chosen simply for illustrative purposes.

We see that, for the chosen values of parameters $C$ and $s$, the desired
qualitative behavior is observed, namely, for small numbers of walkers while
the aggregate is increasing rapidly in size, the non-Newtonian correction is
large, providing a sizeable non-Newtonian effect. This mimics the large shear-
thinning effect anticipated at early times in an experiment, where velocity
and shear-rate are large. At later times however, the (now large) aggregate is
growing more slowly; one would anticipate a low shear-rate in the analogous
physical fluid experiment, and the non-Newtonian correction to the sticking
probability is accordingly much smaller.

Table 1 summarizes the parameter values used in our simulations.

Figure S2: The maximum non-Newtonian probability probability $p_{\rm NN}$ of Eq. (2) in the main text, running over each period of increasing length shown by the data points. The Newtonian parameters used are $A=4$ and $B=0.4$. Panel (a) uses $C=0.3,$ $s=0.2$ while panel (b) uses $C=0.5,$ $s=0.15$.  Table 1: List of parameters used in simulations. Parameter | Description | (Range of) Value(s)
---|---|---
$A$ | Surface tension sticking probability | $\left[1,4\right]$
$B$ | Constant sticking probability | $[0.4,1]$
$l$ | Radius of walkers’s neighborhood (nbd) | 17
$n_{0}$ | Occupied sites in a semi-circular nbd | 89/193
$\left(C,s\right)$ | Parameters defining $p_{\rm NN}$ | $(0,1)\times(0,1)$

##  III Further discussion of the attachment rules modeling shear thinning
behavior

In the main text, we discuss the concept of ‘age’. Here, we discuss the
necessary algorithm modification relevant when the walker is in contact with
exactly three occupied lattice sites when it attaches. Depending on how the
walker is attached, the average of either their horizontal or vertical
neighbors’ ages are used (e.g., for such a horizontal attachment, the
horizontal age of the newly-attached walker is inherited from the single
neighbor in the aggregate in horizontal contact, while vertical age is
obtained by averaging over the two vertical neighbors). Figures S3 and S4
illustrate the algorithm for a few typical cases.

Figure S3: Illustration of attachment rules modeling shear thinning behavior.
White squares represent occupied lattice sites, while black ones represent
unoccupied sites. In all panels, age is written in the form
$\left(T_{H},~{}T_{V}\right)$ discussed in the main text, and the $k=1$ seed
particle located at the origin has an age of (1,1). Panels (a) and (b) show
simple examples of the horizontal and vertical age evolution of the aggregate
in one dimension. Panels (c) and (d) show more elaborate examples where both
the horizontal and vertical ages of the aggregate must be accounted for.
Figure S4: An example where a walker, $(k=6),$ wants to attach at a lattice
site in contact with three occupied sites. The red arrow indicates that this
walker intends to attach horizontally from the left. In this case, the
vertical age $T_{V}$ is the average of the vertical ages resulting from the
$k=3$ and $k=5$ occupied sites, found from $((6-1)+(6-3))/2=4$, while the
horizontal age $T_{H}$, found from $6-1=5$, is updated solely from the $k=1$
site.

##  IV Calculations of fractal dimension

We characterize the morphologies of the resulting aggregates by computing
their fractal dimensions. Since the fractal dimension is the crucial measure
that we use for the quantification of the emerging patterns, we describe the
algorithms implemented in detail here: box-counting (Sec. IV.1) and
correlation length (Sec. IV.2) based.

###  IV.1 Box-counting Dimension

The box-counting dimension of a bounded set, also known as the Minkowski-
Bouilgand dimension [29], is defined as

$\displaystyle D_{{\rm f,box}}:=\lim_{\epsilon\rightarrow
0}\frac{\log{N\left(\epsilon\right)}}{\log{1/\epsilon}},$ (S2)

where $N\left(\epsilon\right)$ is the number of boxes of side length
$\epsilon$ required to cover the bounded set in a metric space.

In practice, the limit specified by Eq. (S2) is difficult to compute. To
approximate it, we consider a nonzero limiting value of $\epsilon$ in Eq.
(S2), and multiply both sides by $\log\left(1/\epsilon\right)$ to obtain a
rough estimate

$N\left(\epsilon\right)\approx c\epsilon^{-D_{{\rm f,box}}},$ (S3)

where $c$ is a constant, and $N\left(\epsilon\right)$ can be further
understood as the number of $\epsilon$-boxes required to cover the aggregate.
Then, we calculate the slope of the log-log linear fit of
$N\left(\epsilon\right)$ versus side length $\epsilon$. To facilitate this
calculation, we use boxes $B_{\epsilon}(\vec{x})$ centered about a data point
$\vec{x}$, defined by the $l^{\infty}$-norm,

$B_{\epsilon}\left(\vec{x}\right)=\left\\{\vec{y}:\left\|\vec{x}-\vec{y}\right\|_{\infty}\leq\epsilon\right\\}.$
(S4)

The numerical approach is summarized by the Algorithm 1. For the domain size
$N\times N,$ we use intermediate box sizes ranging from $N/4$ to $3N/4.$

Computing the Box Counting Dimension

  1. 1.

Begin with an on-lattice box $B_{\epsilon}$ of side length $\epsilon$ and with
a corner that coincides with the top left corner of the domain.

  2. 2.

If the box has any particles in it or on the boundary, count it as a cover.

  3. 3.

Sweep across the domain from left to right, top to bottom, without
overlapping. Record the total number of covers $N\left(\epsilon\right)$.

  4. 4.

Decrease the size of the box $B_{\epsilon}$.

  5. 5.

Repeat Steps 1-4 until enough data pairs
$\left(\epsilon,N\left(\epsilon\right)\right)$ are collected.

  6. 6.

Perform a linear fit of the logarithm of the number of covers versus logarithm
of the size of the corresponding cover. The magnitude of the best-fit line’s
slope is the box-counting dimension.

Algorithm 1 Table 2: Key algorithm parameters used in Algorithms 2 and 3. Parameter | Description
---|---
$D_{{\rm f,box}}$ | Box counting fractal dimension
$D_{{\rm f,corr}}$ | Correlation fractal dimension
$\epsilon$ | Box side length
$N\left(\epsilon\right)$ | Number of boxes of side length $\epsilon$ to cover the aggregate
$C\left(\epsilon\right)$ | The correlation integral
$C_{N}\left(\epsilon\right)$ | An estimator of the correlation integral

###  IV.2 Correlation Length-based Dimension

The second approach we use to calculate the fractal dimension is through the
correlation dimension [30]. First, we form the correlation integral
$C(\epsilon)$ (adopting the nomenclature used by [31]), understood as the
average probability of finding two particles within a ball of radius
$\epsilon$, defined via a limit,

$C\left(\epsilon\right)=\lim_{N\rightarrow\infty}\frac{1}{N\left(N-1\right)}\sum_{\scriptstyle
i,j=1\atop\scriptstyle i\neq
j}^{N}\Theta\left(\epsilon-\left\|\vec{x}\left(i\right)-\vec{x}\left(j\right)\right\|_{l^{2}}\right),\quad\vec{x}\left(i\right)\in\mathbb{R}^{2},$
(S5)

where $\Theta$ is the Heaviside step function. The correlation dimension
$D_{{\rm f,corr}}$ is then defined by

$D_{{\rm f,corr}}=\lim_{N\rightarrow\infty}\lim_{\epsilon\rightarrow
0^{+}}\frac{\log C\left(\epsilon\right)}{\log\epsilon}.$ (S6)

In practice, we estimate $C\left(\epsilon\right)$ by the finite correlation
sum

$C_{N}\left(\epsilon\right)=\frac{1}{N\left(N-1\right)}\sum_{\scriptstyle
i,j=1\atop\scriptstyle i\neq
j}^{N}\Theta\left(\epsilon-\left\|\vec{x}\left(i\right)-\vec{x}\left(j\right)\right\|_{2}\right),\quad\vec{x}\left(i\right)\in\mathbb{R}^{2},$
(S7)

which is a viable approach due to the large number of particles we simulate,
thus producing a very large $N$ and a reasonable approximation of Eq. (S5).
Therefore, to estimate the correlation dimension using
$C_{N}\left(\epsilon\right)$, we use

$D_{{\rm f,corr}}\approx\lim_{\epsilon\rightarrow 0^{+}}\frac{\log
C_{N}\left(\epsilon\right)}{\log\epsilon}.$ (S8)

The approach to computing the correlation sum $C_{N}\left(\epsilon\right)$ is
outlined by the following:

Computing the Correlation Dimension

  1. 1.

Choose the largest search radius $\epsilon$ from a preset range.

  2. 2.

Select a random particle $\vec{x}\left(i\right)$.

  3. 3.

Count the number of unique pairs of particles within a Euclidean $l^{2}$
distance of $\epsilon$ from the particle at $\vec{x}\left(i\right)$.

  4. 4.

Remove $\vec{x}\left(i\right)$ from the aggregate (to avoid double counting)
and go to step 2 until all particles are used.

  5. 5.

Compute a total count $C_{N}\left(\epsilon\right)$.

  6. 6.

Restore aggregate and go to step 1 using a smaller value of $\epsilon$.

Algorithm 2

A subtlety in the implementation of the correlation integral calculation is to
avoid double counting in the double-indexed sum in Eq. (S7); double-counting
arises from the fact that the intersection of $\epsilon$-neighborhoods of
$\vec{x}_{i}$ and $\vec{x}_{j}$, $B_{\epsilon}\left(\vec{x}_{i}\right)\cap
B_{\epsilon}\left(\vec{x}_{j}\right)$ is in general non-empty. To avoid
double-counting, after we finish counting for $\vec{x}_{i}$, we remove this
particle from the aggregate since we have counted all its unique pairs within
its $\epsilon$-neighborhood (step 4 of Algorithm 2). A technicality of this
implementation is that, since this removal procedure alters the aggregate, we
must return to the original aggregate every time we change the search radius
$\epsilon$.

##  V Additional Results

###  V.1 Reproducibility

In this section, we discuss the reproducibility of the fractal dimension
results. To show that the results are realization independent, we set the
parameters $A=1$ and $B=0.4$, and carry out two independent realizations of
Newtonian and non-Newtonian simulations each (non-Newtonian parameters $C=0.9$
and $s=0.07$). Figure S5 shows the results of this numerical experiment for
both fractal dimension calculation methods. The results, at sufficiently large
aggregate sizes, are found to be consistent.

Figure S5: A consistency test for 2 sample simulations with Box Counting
(circle “$\circ$”) and Correlation dimension (“$\times$”). Panels (a) and (b)
are Newtonian simulations (red). Panels (c) and (d) are non-Newtonian
simulations (blue) with $C=0.9$ and $s=0.07$. In all cases, we use $A=1$ and
$B=0.4$.

###  V.2 Influence of non-Newtonian parameters $(C,s)$

The parameters $(C,s)$ specify shear-thinning response, therefore modeling the
rheology of the displaced fluid. Since we find that $D_{\rm f}$ differs
between Newtonian and non-Newtonian simulations, it is not surprising that for
different values of $(C,s)$ one could expect to find different values of
$D_{\rm f}$. Table 3 shows that this is indeed the case; in addition, this
table shows that that the change of $D_{\rm f}$, $\Delta D_{\rm
f,\dagger}=D_{{\rm f,\dagger}}^{{\rm Newt}}-D_{{\rm f,\dagger}}^{{\rm nNewt}}$
(where $\dagger$ is ${\rm box}$ or ${\rm corr}$) is always non-negative, as
pointed out also in the main text.

$\left(D_{{\rm f,box}}^{{\rm Newt}},D_{{\rm f,corr}}^{{\rm
Newt}}\right)=\left(1.7132,1.7033\right)$ $C$ $s$ $D_{\rm f,box}^{\rm nNewt}$
$D_{\rm f,corr}^{\rm nNewt}$ 0.9000 0.07000 1.6807 1.6795 0.3750 0.06250
1.6912 1.6872 0.2250 0.03750 1.6902 1.6887 0.3750 0.1250 1.6978 1.6932 0.2250
0.1000 1.7004 1.6943 0.2250 0.1625 1.7038 1.7005 0.3000 0.1875 1.7039 1.6995
0.2250 0.1750 1.7038 1.7005 0.3000 0.2000 1.7030 1.7002 0.2250 0.1875 1.7052
1.7026 0.2250 0.2000 1.7056 1.7043 0.3000 0.2250 1.7058 1.7040 0.2250 0.2125
1.7045 1.7030 0.2250 0.2250 1.7048 1.7012 0.3000 0.2500 1.7060 1.7020

$\left(D_{{\rm f,box}}^{{\rm Newt}},D_{{\rm f,corr}}^{{\rm
Newt}}\right)=\left(1.7506,1.7613\right)$ $C$ $s$ $D_{\rm f,box}^{\rm nNewt}$
$D_{\rm f,corr}^{\rm nNewt}$ 0.9000 0.07000 1.7079 1.7002 0.3750 0.06250
1.6885 1.6850 0.2250 0.03750 1.6917 1.6841 0.3750 0.1250 1.7208 1.7159 0.2250
0.1000 1.7260 1.7256 0.2250 0.1625 1.7349 1.7405 0.3000 0.1875 1.7366 1.7419
0.2250 0.1750 1.7357 1.7397 0.3000 0.2000 1.7376 1.7400 0.2250 0.1875 1.7342
1.7426 0.2250 0.2000 1.7350 1.7406 0.3000 0.2250 1.7365 1.7420 0.2250 0.2125
1.7370 1.7434 0.2250 0.2250 1.7387 1.7435 0.3000 0.2500 1.7369 1.7439

$\left(D_{{\rm f,box}}^{{\rm Newt}},D_{{\rm f,corr}}^{{\rm
Newt}}\right)=\left(1.7609,1.7759\right)$ $C$ $s$ $D_{\rm f,box}^{\rm nNewt}$
$D_{\rm f,corr}^{\rm nNewt}$ 0.9000 0.07000 1.7526 1.7581 0.5000 0.05000
1.6833 1.6792 0.3000 0.03000 1.6867 1.6820 0.5000 0.1000 1.7043 1.6922 0.3000
0.08000 1.6960 1.6883 0.5000 0.1500 1.7083 1.6968 0.3000 0.1300 1.7189 1.7143
0.3000 0.1400 1.7263 1.7293 0.4000 0.1600 1.7280 1.7283 0.3000 0.1500 1.7352
1.7423 0.4000 0.1800 1.7407 1.7479 0.3000 0.1600 1.7444 1.7516 0.5000 0.2000
1.7417 1.7542 0.3000 0.1700 1.7456 1.7561 0.3000 0.1800 1.7487 1.7580

Table 3: Values of $\left(C,s,D_{\rm f,box}^{\rm nNewt},D_{\rm f,corr}^{\rm
nNewt}\right)$ for $A=1$ (left), $A=2$ (middle), and for $A=4$ (right). Each
table title provides the corresponding box counting and correlation fractal
dimension obtained from Newtonian simulations.

## References

  * [1] Henry Selby Hele-Shaw.  Flow of water.  Nature, 58(1509):520–520, 1898.
  * [2] P.G. Saffman and G.I. Taylor.  The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid.  Proc. R. Soc. Lond., 245:312–329, 1958.
  * [3] K.V. McCloud and J.V. Maher.  Experimental Perturbations to Saffman-Taylor Flow.  Phys. Rep., 260:139, 1995.
  * [4] T. C. Halsey.  Diffusion‐limited aggregation: A model for pattern formation.  Phys. Today, 53:36–41, 2000.
  * [5] J. Casademunt.  Viscous fingering as a paradigm of interfacial pattern formation: Recent results and new challenges.  Chaos, 14:809–824, 2004.
  * [6] A. Vasil’ev.  From the Hele-Shaw experiment to integrable systems: a historical overview.  Complex Anal. Operator Theory, 3:551–585, 2009.
  * [7] G.M. Homsy.  Viscous Fingering in Porous Media.  Ann. Rev. Fluid Mech., 19:271, 1987.
  * [8] Capstone Laboratory: Hele-Shaw Project.  https://cfsm.njit.edu/capstone/projects/2016/main.php, 2016.
  * [9] C.A. Hieber.  In A.I. Isayev, editor, Injection and Compression Molding Fundamentals. Marcel Dekker, 1987.
  * [10] J.S. Langer.  In H.E. Stanley, editor, Statistical Physics. Noth-Holland, 1986\.
  * [11] L. Paterson.  Radial fingering in a hele-shaw cell.  J. Fluid Mech., 113:513–529, 1981.
  * [12] H. van Damme, F. Obrecht, P. Levitz, L. Gaineau, and C. Laroche.  Fractal viscous fingering in clay slurries.  Nature, 320, 1986.
  * [13] E. Ben-Jacob, G. Deutscher, P. Garik, Nigel D. Goldenfeld, and Y. Lareah.  Formation of a dense branching morphology in interfacial growth.  Phys. Rev. Lett., 57:1903–1906, 1986.
  * [14] S. E. May and J. V. Maher.  Fractal dimension of radial fingering patterns.  Phys. Rev. A, 40:1723–1726, Aug 1989.
  * [15] H. Zhao and J. V. Maher.  Viscoelastic effects in patterns between miscible liquids.  Phys. Rev. A, 45:R8328–R8331, 1992.
  * [16] A. Lindner, D.Bonn, E. C. Poiré, M. Ben Amar, and J. Meunier.  Viscous fingering in non-Newtonian fluids.  J. Fluid Mech., 469:237–256, 2002.
  * [17] G. Daccord, J. Nittmann, and H. E. Stanley.  Radial viscous fingers and diffusion-limited aggregation: Fractal dimension and growth sites.  Phys. Rev. Lett., 56:336–339, 1986.
  * [18] L. Kondic, P. Palffy-Muhoray, and M. J. Shelley.  Models of non-Newtonian Hele-Shaw flow.  Phys. Rev. E, 54:4536, 1996.
  * [19] L. Kondic, M. J. Shelley, and P. Palffy-Muhoray.  Non-Newtonian Hele-Shaw flow and the Saffman-Taylor instability.  Phys. Rev. Lett, 80:1433, 1998.
  * [20] M. Ben Amar and E. C. Poiré.  Pushing a non-Newtonian fluid in a Hele-Shaw cell: From fingers to needles.  Phys. Fluids, 11:1757–1767, 1999.
  * [21] P. Fast, L. Kondic, M. J. Shelley, and P. Palffy-Muhoray.  Pattern formation in non-Newtonian Hele-Shaw flow.  Phys. Fluids, 13:1191, 2001.
  * [22] T. A. Witten and L. M. Sander.  Diffusion-limited aggregation.  Phys. Rev. B, 27:5686–5697, 1983.
  * [23] T. Vicsek.  Pattern formation in diffusion-limited aggregation.  Phys. Rev. Lett., 53:2281–2284, 1984.
  * [24] M. Blunt and P. King.  Scaling structure of viscous fingering.  Phys. Rev. A, 37:3935–3941, 1988.
  * [25] B. B. Mandelbrot, B. Kol, and A. Aharony.  Angular gaps in radial diffusion-limited aggregation: Two fractal dimensions and nontransient deviations from linear self-similarity.  Phys. Rev. Lett., 88:055501, 2002.
  * [26] J. Lee, A. Coniglio, and H. E. Stanley.  Fractal-to-nonfractal crossover for viscous fingers.  Phys. Rev. A, 41:4589–4592, 1990.
  * [27] J. Adriazola, B. Gu, L.J. Cummings, and L. Kondic, 2024.  See Supplemental Material at [URL will be inserted by publisher].
  * [28] D. N. Bankar, P. M. Gade, A. V. Limaye, and A. G. Banpurkar.  Segregation of fractal aggregates grown from two seeds.  Phys. Rev. E, 75:051401, 2007.
  * [29] K. Falconer.  Fractal geometry: mathematical foundations and applications.  John Wiley & Sons, 2004.
  * [30] S. H. Strogatz.  Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.  CRC press, 2018.
  * [31] P. Grassberger and I. Procaccia.  Measuring the strangeness of strange attractors.  Physica D, 9:189–208, 1983.