# Dense Suspension Flow in a Penny-Shaped Crack
Part I : Theory

George R. Wyatt1 & Herbert E. Huppert2 1Emmanuel College, St. Andrew’s Street,
Cambridge, CB2 3AP<EMAIL_ADDRESS>
2King’s College, King’s Parade, Cambridge, CB2 1ST<EMAIL_ADDRESS>
###### Abstract.

We study the dynamics of proppants carried by fluid driven into an evolving
penny-shaped fracture. The behaviour of the slurry flow is investigated in two
phases: pressurised injection and elastic closure. During injection the slurry
is modelled using a frictional rheology that takes into account the shear-
induced migration and jamming of the proppants. Making pragmatic assumptions
of negligible toughness and cross-fracture fluid slip, we find self-similar
solutions supporting a range of proppant concentration profiles. In
particular, we define an effective viscosity, which equates the fracture
evolution of a slurry flow with a given proppant volume fraction, to a
Newtonian flow with a particular viscosity. Using this framework, we are able
to make predictions about the geometry of the growing fracture and the
significance of tip screen-out. In the closure phase, proppants are modelled
as incompressible and radially immobile within the narrowing fracture. The
effects of proppant concentration on the geometry of the residual propped
fracture are explored in full. The results have important applications to
industrial fracking and geological dike formation by hot, intruding magma.

###### Key words and phrases:

Hydraulic fracture, suspension flow, rheology, proppant transport, elastic,
tip screen-out, penny-shaped, cavity flow

##  1\. Introduction

Receiving a patent for his ‘exploding torpedo’ in 1865, US Civil War veteran
Col. Edward Roberts established the practice of fracturing bedrock to
stimulate oil wells [1]. A technique, known as hydraulic fracturing, which
uses pressurised fluid rather than explosives to develop fracture networks,
only came into practice much later, in 1947 [2], and is the topic of this
paper. In particular, we will concentrate on the convective transport of
proppants within an evolving cavity. These are small particles added to the
fracturing fluid in order to prop open the developed fracture, which closes
under far-field stress once the fluid pressure is released. Aside from its use
in hydrocarbon recovery, hydraulic fracturing, or fracking, has uses including
the measurement of in-situ stresses in rocks [3], generation of electricity in
enhanced geothermal systems [4] and improvement of injection rates in CO2
sequestration [5]. Hydraulic fracturing processes are also ubiquitous in
geology: dikes and sills arise from cracks whose growth is driven by magma,
with magmatic crystals taking the place of synthetic proppants. Phenomena such
as crystallisation and gas exsolution in the cooling magma mean models of dike
propagation vary widely, as is summarised in [6]. Notably, Petford & Koenders
[7] utilise granular flow theory to model the ascent of a granitic melt
containing solids.

This paper combines two significant, but often disconnected, fields of
fracking study, cavity flow and suspension flow:

  * •

The study of (elastohydrodynamic) cavity flow focusses on the interplay
between hydrodynamic properties of the fracturing fluid and material
properties of the medium being fractured. In the zero-proppant case, the
problem of a fluid-driven, penny-shaped crack requires the joint solution of a
nonlinear Reynold’s equation, which governs flow within the crack, and a
singular integral boundary condition, which takes into account the elastic
properties of the surrounding medium. The general strategy used in this paper
takes inspiration from the work of Spence & Sharp [8], who in 1985,
restricting to the two-dimensional case, were the first to solve these
integro-differential equations. In particular, we will focus on cavities that
keep the same shape in some evolving coordinate system, using series
expansions to represent both the width and pressure profiles within the
fracture. More recently, in 2002, Savitski & Detournay [9] solved similar
three-dimensional versions of these equations, allowing them to find fracture
evolutions with simple time dependence in both the viscous and toughness
dominated regimes. In the former, the principal energy dissipation is by
viscous flow, and in the latter, energy dissipation is mostly by creating new
fracture surfaces. Notably, the same paper [9] verifies that industrial
fracking occurs in the viscous regime; this assumption makes the problem
considered in this paper tractable to a semi-analytical approach.

  * •

The mathematical study of suspension flow dates back to 1906, when Einstein
used properties of suspensions to estimate the size of a water molecule [10].
In particular, he showed that very dilute particle-laden flows are Newtonian,
with a viscosity which increases with the concentration of particles. However,
during hydraulic fracturing it is necessary to model a full range of proppant
volume fractions, which we denote by $\phi$. It is typical to have both dilute
flow near the crack walls, as well as plug flow at the centre of the cavity,
where the slurry behaves as a porous granular medium. More recent experiments
by Boyer et al. in 2011 [11] investigate dense suspension rheology. They show
that particles in suspension, subject to a constant normal particle pressure
that is applied by a porous plate, expand when a shear is applied to the
mixture. As a result, it is possible to write $\phi=\phi(I)$, where the
dimensionless parameter, $I$, is the ratio between the fluid shear stress,
which is proportional to the shear rate, and the particle normal stress.
Likewise, fixing the solid volume fraction, they showed that the normal
particle pressure is proportional to the mixture shear stress. It is also
shown that the constant of proportionality, $\mu$, can be expressed as a
decreasing function of $\phi$. In the same paper [11], forms of the
rheological functions $I$ and $\mu$ are suggested, showing good agreement with
experimental data. Since then, several papers have suggested slightly
different rheological models and are reviewed by Donstov et al. in [12]. These
all feature a jamming limit, $\phi_{m}$, which is the volume fraction at which
the flowing slurry transitions into a granular solid. We will utilise the
frictional rheology given by Lecampion & Garagash [13], which is unique in
allowing packings with $\phi>\phi_{m}$. These denser packings form due to ‘in-
cage’ particle rearrangements caused by velocity and pressure fluctuations in
the surrounding flow.

The endeavours of this paper may be condensed into three main objectives. The
first is to establish a mathematical framework that captures the behaviour of
the proppant suspension as it interacts with the growing cavity. Here we will
utilise a lubrication model, along with the assumption that the proppant flow
is fully developed; equivalently, that the transverse fluid slip is
negligible. Crucially, we will try to justify these assumptions using typical
parameters from industrial fracking. We will also make a zero-toughness
assumption, which is validated in [9]. Once we have developed this framework,
an important step will be to compare its features to those derived in the
zero-proppant, viscosity dominated case by Savitski & Detournay [9],
particularly because we utilise a frictional rheology fitted to the dense
regime. The second objective is to find and examine accurate numerical
solutions modelling the developing cavity, given a range of proppant
concentrations. We will explore the empirical effects of changing proppant
concentration on the geometry of the developing fracture, as well as the
distribution of proppants. Where possible, we will evaluate the consistency of
our model and forecast potential shortfalls such as proppant screen-out near
the crack tip. The third, and final, objective is to leverage our results to
make predictions about the geometry of the fracture after the fluid pressure
is released. By assuming the remaining proppants are immobile and
incompressible, we aim to establish simple formulae predicting the width and
radius of the developed fracture. Since these relate directly to the
conductivity of the formation, this third objective is potentially the most
significant.

Aside from the availability of semi-analytical solutions, the problem of
proppant flow in a penny-shaped crack is particularly appealing because of the
potential of practical verification. Recent experiments by O’Keeffe, Huppert &
Linden [14] have explored fluid-driven, penny-shaped fractures in transparent,
brittle hydrogels, making use of small particle concentrations to measure in-
crack velocities. This paper is the first of two; the second of which will be
a practical treatise on slurry driven-fractures in hydrogels, aiming to verify
the predictions made here by repeating the experiments of [14] including
proppant concentrations.

##  2\. Injection: Problem Formulation

Figure 1. Schematic of the penny-shaped crack.

###  2.1. Fracture Mechanics

We model the propagation of a penny-shaped crack similar to that shown in
Figure 1, using the framework of Detournay & Savitski [9]. We will make the
following assumptions:

  * •

The crack is axisymmetric and has reflectional symmetry in $z=0$, with half
width $w(r,t)$ and total radius $R(t)$, so $w(R,t)=0$.

  * •

The fluid is injected from a point source, with the wellbore radius negligible
compared to the fracture radius.

  * •

The lag between the fracture tip and the fluid front is negligible compared to
the fracture radius.

  * •

The fracture propagates in continuous mobile equilibrium.

  * •

The normal stress on the fracture walls due to proppants is negligible
compared to the fluid pressure.

The third assumption is validated by Garagash & Detournay [15] and introduces
a negative pressure singularity at the tip of the crack ($r=R$). The fourth
and fifth assumptions lead to the following integral equations from linear
elastic fracture mechanics. These relate the net fluid pressure, $p(r,t)$, to
the opening of the fracture and the toughness of the surrounding rock.

(1) $\displaystyle w(r,t)$ $\displaystyle=\frac{4R}{\pi
E^{\prime}}\int_{r/R}^{1}\frac{y}{\sqrt{y^{2}-(r/R)^{2}}}\int_{0}^{1}\frac{xp(xyR,t)}{\sqrt{1-x^{2}}}dxdy,$
(2) $\displaystyle K_{Ic}$ $\displaystyle=\frac{2}{\sqrt{\pi
R}}\int_{0}^{R}\frac{p(r,t)r}{\sqrt{R^{2}-r^{2}}}dr,$

where $E^{\prime}$ is the plane strain modulus, given by the Young modulus,
$E$, and the Poisson ratio, $\nu$, as $E^{\prime}=E/(1-\nu^{2})$. $K_{Ic}$ is
the material toughness. These equations can be attributed to Sneddon [16] and
Rice [17] respectively. We note that $p$ represents the fluid pressure minus
the in-situ stress of the surrounding rock, which is assumed to be isotropic.
We write $p$ with radial spatial dependence only; this will be validated
later, along with the fifth assumption, using a lubrication argument.

###  2.2. Frictional Rheology

We model the injected flow as a Newtonian fluid containing identical spherical
particles. Recent approaches in modelling dense slurry flow are characterised
by empirical relations originally proposed by Boyer et al. [11]. The first of
these relates the fluid shear stress to the normal stress required to confine
the particles; the second gives the ratio of the mixture shear stress to the
particle confining stress,

(3) $\displaystyle I(\phi)$
$\displaystyle=\eta_{f}\dot{\gamma}/\sigma_{n}^{s},$ $\displaystyle\mu(\phi)$
$\displaystyle=\tau/\sigma_{n}^{s}.$

Here $\eta_{f}$ is the carrying fluid’s dynamic viscosity, $\phi$ is the
volume fraction of the proppants, $\dot{\gamma}$ is the solid shear rate and
$\sigma_{n}^{s}$ is the normal particle stress, which we will sometimes refer
to as the particle pressure. The second ratio is given the symbol $\mu$, not
to be confused with dynamic viscosity, because it resembles a friction
coefficient. These relations are given a clear experimental grounding in [11],
which is discussed in the introduction. Various forms of the dimensionless
functions $I(\phi)$ and $\mu(\phi)$ have been compared to experimental results
in [12] using the equivalent formulation:
$\tau=\eta_{s}(\phi)\eta_{f}\dot{\gamma}$ and
$\sigma_{n}=\eta_{n}(\phi)\eta_{f}\dot{\gamma}$, where
$\eta_{s}=\mu(\phi)/I(\phi)$ and $\eta_{n}=1/{I(\phi)}$.

In our calculations we will utilise the frictional rheology provided by B.
Lecampion & D. I. Garagash [13], which is unique in allowing packings with
volume concentrations greater than $\phi_{m}$. Here $I(\phi)=0$, meaning the
proppants have zero shear rate and effectively resemble a permeable solid.
Explicitly, we use the expressions

(4)
$\displaystyle\mu=\mu_{1}+\frac{\phi_{m}}{\delta}\left(1-\frac{\phi}{\phi_{m}}\right)$
$\displaystyle+\left(I(\phi)+\left[\frac{5}{2}\phi_{m}+2\right]I(\phi)^{0.5}\right)\left(1-\frac{\phi}{\phi_{m}}\right)^{2},$
(7) $\displaystyle I(\phi)$
$\displaystyle=\left\\{\begin{array}[]{rl}\left(\phi_{m}/\phi-1\right)^{2}&\textrm{
if }\phi<\phi_{m}\\\ 0&\textrm{ if }\phi\geq\phi_{m},\end{array}\right.$

where $\phi_{m}=0.585$, $\mu_{1}=0.3$ and $\delta=0.158$; these are plotted in
Figure 2. We might have used a different rheology, but this model shows good
agreement with the data of Boyer et al. [11] and Dagois-Bohy et al. [18] for
$0.4<\phi<\phi_{m}$. Furthermore, owing to its linear extension beyond
$\phi_{m}$, $\mu$ is a simple monotonic function, meaning we can invert it
easily to find $\phi$. In other models $\phi(\mu)$ is constant for
$\mu<\mu(\phi_{m})$; this means that $\phi_{m}$ is the maximum volume
fraction, regardless of how small shear stresses in the jammed slurry become.
An important observation is that $\mu=0$ implies
$\phi=\phi_{m}+\delta\mu_{1}\approx 0.63\approx\phi_{rcp}$. Here $\phi_{rcp}$
is the random close packing limit, the maximal observed volume fraction due to
random packing. This reflects the fact that, for a given confining stress, as
the shear stress tends to zero, the particles pack to this maximal density.

This rheology uses a continuum model that requires particles to be small
compared to the size of the fracture. This is in order to well-define the
proppant volume concentration, $\phi$. In our model the relevant ratio is that
of the particle diameter to the typical crack width, the smallest cavity
length scale. In [13], good results are obtained using the same rheological
model, with this ratio taking values as large as $1/10$. However, as the ratio
approaches unity we have to consider non-local effects, such as proppant
bridging across the crack width. This is particularly important near the
fracture tip, where $w$ approaches zero. These effects will be discussed in
greater detail in Section 7, once we have formed a model of the evolving
fracture. We must also be cautious applying these rheological models to dilute
flows, since they are fitted to experimental data from the dense regime, where
$\phi>0.4$. This difficulty is somewhat inevitable, since the determination of
$I$ and $\mu$ requires measurement of the particle normal stress, or particle
pressure, which becomes very small in the dilute regime.

(a) $I$

(b) $\mu$

(c) $I/\mu$

(d) $I/\mu$ data

Figure 2. Plots of the rheological functions $I$, $\mu$ and $I/\mu$ given by
Lecampion & Garagash [13]. Also plotted is the experimental data of Boyer et
al. [11] using polystyrene spheres of diameter 580$\mu$m in $2.15$Pa s fluid
(red), as well as poly(methyl methacrylate) spheres of diameter 1100$\mu$m
suspended in $3.10$Pa s fluid (orange); and of Dagois-Bohy et al. [18] using
polystyrene spheres of diameter 580$\mu$m suspended in $2.27$Pa s fluid
(purple). All experiments are carried out with a fixed particle pressure,
applied by a porous plate.

###  2.3. Fluid Slip

We define $\mathbf{u}$ as the slurry velocity, $\mathbf{v}$ as the particle
velocity and $\mathbf{q}=\mathbf{u}-\mathbf{v}$ as the slip velocity. We then
employ the slip relation

(8) $\displaystyle\mathbf{q}$
$\displaystyle=\frac{a^{2}\kappa(\phi)}{\eta_{f}}\nabla\cdot\sigma^{f},$ (9)
$\displaystyle\kappa(\phi)$ $\displaystyle=\frac{2(1-\phi)^{5.1}}{9\phi},$

where $a$ is the particle radius and $\sigma_{f}$ is the fluid stress tensor.
Since fluid and particle shear rates are often similar, we ignore fluid shear
stresses and take $\sigma^{f}=-pI$; this is typical in the analysis of porous
media flow. This simplifies (8) to Darcy’s law. However, the effect of fluid
shear stress is taken into account in the frictional rheology, where it is
included as part of the solid shear stress. $\kappa$ is a normalised form of
the permeability of the solid particles; we use the function suggested by
Garside & Al-Dibouni [19], which is based on the phenomenology first described
by Richardson & Zaki [20]. This choice of permeability function shows
excellent agreement with the experimental results of Bacri et al. [21].

###  2.4. Conservation Equations

We consider the effective Reynolds number,

(10) $\displaystyle\textrm{Re}_{\textrm{eff}}=\frac{\rho
u_{r}w^{2}}{\eta_{f}R},$

to be negligible. We also neglect the effect of gravity, since we are mainly
concerned with small or neutrally buoyant proppants, which settle slowly.
Hence, our momentum balance becomes

(11) $\displaystyle\nabla\cdot\sigma=0,$

where $\sigma=\sigma^{s}+\sigma^{f}$ is the mixture stress tensor, composed of
the particle and fluid stresses respectively. We also note that, subtracting
the hydrostatic pressure term, we write $\sigma=\tau-pI$. Since we assumed
$\sigma^{f}=-pI$ in deriving the fluid slip equation, we deduce
$\sigma_{s}=\tau$. This is a notational quirk arising from the frictional
rheology because $\tau$ does include shear stress originating from the viscous
carrier fluid. Herein we will refer to $\sigma^{s}_{zz}$ and $\tau_{rz}$,
since the former generally arises from the proppants and the latter stems from
both the proppants and the carrier fluid. The assumption of axisymmetry gives

(12) $\displaystyle\frac{1}{r}\frac{\partial(r\tau_{rr})}{\partial
r}+\frac{\partial\tau_{rz}}{\partial z}-\frac{\partial p}{\partial r}$
$\displaystyle=0,$
$\displaystyle\frac{1}{r}\frac{\partial(r\tau_{rz})}{\partial
r}+\frac{\partial\sigma^{s}_{zz}}{\partial z}-\frac{\partial p}{\partial z}$
$\displaystyle=0.$

We also have the continuity equations

(13) $\displaystyle\nabla\cdot(\mathbf{v}+\mathbf{q})$ $\displaystyle=0,$
$\displaystyle\frac{\partial\phi}{\partial t}+\nabla\cdot(\phi\mathbf{v})$
$\displaystyle=0.$

The first of these can be integrated over the fracture volume to give
$Qt=4\pi\int_{0}^{R}rw(r,t)dr.$ Here, $Q$ is the rate at which the slurry is
pumped into the crack, which we will assume is constant. We will also assume
that the proppants are injected at a constant rate, meaning the average
concentration at the wellbore is constant.

##  3\. Injection: Scalings

To help implement the assumptions of a lubrication model, where the crack
width is far smaller than the crack radius, we introduce the scaled
coordinates,

$\displaystyle T$ $\displaystyle=T(t),$ $\displaystyle r$
$\displaystyle=L(t)\Gamma(T)\xi,$ $\displaystyle z$
$\displaystyle=\epsilon(t)L(t)\eta.$

Here $T(t)$ is the internal time scale, a monotonic function to be specified
later; $\epsilon(t)$ is a small number; and $\Gamma(T)$ is the crack radius,
measured in the scaled coordinates, so $\xi=1$ implies $r=R$. We multiply the
variables accordingly,

$\displaystyle w(r,t)$ $\displaystyle\to\epsilon Lw(\xi,T),$ $\displaystyle
p(r,z,t)$ $\displaystyle\to\epsilon E^{\prime}p(\xi,\eta,T),$ $\displaystyle
R(t)$ $\displaystyle\to L\Gamma(T),$ $\displaystyle v_{z}(r,z,t)$
$\displaystyle\to-\dot{\epsilon}Lv_{z}(\xi,\eta,T),$ $\displaystyle
v_{r}(r,z,t)$
$\displaystyle\to\frac{-\dot{\epsilon}L}{\epsilon}v_{r}(\xi,\eta,T),$
$\displaystyle q_{r}(r,z,t)$
$\displaystyle\to\frac{\epsilon}{L}\frac{a^{2}E^{\prime}}{\eta_{f}\Gamma}q_{r}(\xi,\eta,T),$
$\displaystyle q_{z}(r,z,t)$
$\displaystyle\to\frac{1}{L}\frac{a^{2}E^{\prime}}{\eta_{f}}q_{z}(\xi,\eta,T),$
$\displaystyle\tau(r,z,t)$
$\displaystyle\to-\frac{\dot{\epsilon}}{\epsilon^{2}}\eta_{f}\tau(\xi,\eta,T),$
$\displaystyle\sigma^{s}(r,z,t)$
$\displaystyle\to-\frac{\dot{\epsilon}}{\epsilon^{2}}\eta_{f}\sigma^{s}(\xi,\eta,T).$

The appearance of minus signs reflects the fact that $\epsilon$, the ratio of
the characteristic radius to the characteristic width of the fracture, is
decreasing. We also assume the scaling is suitable so that all the scaled
variables are $\mathcal{O}(1)$. Herein, we will use $(\dot{})$ for derivatives
with respect to $t$ and $(^{\prime})$ for those with respect to $T$.

In the new, rescaled coordinates the equations describing the frictional
rheology become $I(\phi)=\dot{\gamma}/\sigma_{n}^{s}$ and
$\mu(\phi)=\tau/\sigma_{n}^{s}$. The slip equation becomes
$\mathbf{q}=-\kappa(\phi)\nabla p,$ where $\nabla$ is now with respect to
$(\xi,\eta)$. The integral equations become

(14) $\displaystyle w(\xi,T)$
$\displaystyle=\frac{4\Gamma}{\pi}\int_{\xi}^{1}\frac{y}{\sqrt{y^{2}-\xi^{2}}}\int_{0}^{1}\frac{xp(xy,T)}{\sqrt{1-x^{2}}}dxdy,$
$\displaystyle\aleph\equiv\frac{K_{Ic}}{\epsilon E^{\prime}\sqrt{L}}$
$\displaystyle=2\sqrt{\frac{\Gamma}{\pi}}\int_{0}^{1}\frac{p(\xi,T)\xi}{\sqrt{1-\xi^{2}}}d\xi.$

The momentum equations are

(15)
$\displaystyle\frac{\epsilon}{\Gamma\xi}\frac{\partial(\xi\tau_{rr})}{\partial\xi}+\frac{\partial\tau_{rz}}{\partial\eta}+\frac{\epsilon^{3}E^{\prime}t}{\eta_{f}}\frac{\epsilon}{\dot{\epsilon}t\Gamma}\frac{\partial
p}{\partial\xi}$ $\displaystyle=0,$
$\displaystyle\frac{\epsilon^{2}}{\Gamma\xi}\partialderivative{(\xi\tau_{rz})}{\xi}+\epsilon\partialderivative{\sigma^{s}_{zz}}{\eta}+\frac{\epsilon}{\dot{\epsilon}t}\frac{\epsilon^{3}E^{\prime}t}{\eta_{f}}\partialderivative{p}{\eta}$
$\displaystyle=0.$

Since we expect the radial pressure gradient to be comparable to the shear
stress, $\tau_{rz}$, we choose $\epsilon$ so that the dimensionless quantity
$\epsilon^{3}E^{\prime}t/\eta_{f}=1$. Finally, the global volume conservation
equation then becomes $Qt/(\epsilon L^{3})=4\pi\Gamma^{2}\int_{0}^{1}\xi
w(\xi,T)d\xi,$ so in a similar manner we choose the dimensionless quantity
$Qt/\epsilon L^{3}=1.$ These choices mean

(16) $\displaystyle\epsilon(t)$
$\displaystyle=(\eta_{f}/E^{\prime})^{\frac{1}{3}}t^{-1/3},$ $\displaystyle
L(t)$ $\displaystyle=(E^{\prime}Q^{3}/\eta_{f})^{\frac{1}{9}}t^{4/9}.$

We will repeatedly use the relations $\dot{\epsilon}t/\epsilon=-1/3$ and
$\dot{L}t/L=4/9$. Using this choice of $\epsilon$ we note that, before
scaling, $\sigma^{s}/p=\mathcal{O}(\epsilon)$; this validates the assumption
that particle pressure is negligible compared to hydrostatic pressure at the
crack walls. Also, by the scaled momentum equations,

(17) $\displaystyle\frac{\partial\tau_{rz}}{\partial\eta}$
$\displaystyle=\frac{3}{\Gamma}\frac{\partial
p}{\partial\xi}+\mathcal{O}(\epsilon),$ $\displaystyle\frac{\partial
p}{\partial\eta}$
$\displaystyle=\frac{\epsilon}{3}\frac{\partial\sigma^{s}_{zz}}{\partial\eta}+\mathcal{O}(\epsilon^{2}),$

the second of which verifies the assumption that $p$ has spatial dependence in
the radial direction only. Because of the $\eta=0$ reflectional symmetry, we
note that $\tau_{rz}(\xi,0)=0$. So, ignoring $\mathcal{O}(\epsilon)$ terms and
integrating (17.1), we see that

(18) $\displaystyle\tau_{rz}=\frac{3\eta}{\Gamma}\frac{\partial
p}{\partial\xi},$

and, using the scaled equations from the frictional rheology,

(19) $\displaystyle\sigma_{zz}^{s}$
$\displaystyle=\frac{3|\eta|}{\Gamma}\frac{1}{\mu(\phi)}\frac{\partial
p}{\partial\xi},$ $\displaystyle\frac{\partial v_{r}}{\partial\eta}$
$\displaystyle=\frac{3\eta}{\Gamma}\frac{I(\phi)}{\mu(\phi)}\frac{\partial
p}{\partial\xi}.$

Then, using the condition $v_{r}(\xi,\pm w)=0$, we deduce that

(20) $\displaystyle
v_{r}(\xi,\eta)=-\frac{3}{\Gamma}\partialderivative{p}{\xi}\int_{\eta}^{w}\frac{I(\phi)\eta}{\mu(\phi)}d\eta.$

##  4\. Injection: Time Regimes

In this choice of scaling, the slurry conservation equation becomes

(21) $\displaystyle\frac{1}{3\Gamma\xi}\frac{\partial(\xi
v_{r})}{\partial\xi}+\frac{1}{3}\frac{\partial
v_{z}}{\partial\eta}+\left(\frac{a}{L\Gamma}\right)^{2}\frac{1}{\epsilon^{2}\xi}\frac{\partial(\xi
q_{r})}{\partial\xi}+\left(\frac{a}{L}\right)^{2}\frac{1}{\epsilon^{4}}\frac{\partial
q_{z}}{\partial\eta}=0.$

Combining this with the scaled slip equation, noting (17), we obtain

(22) $\displaystyle\frac{1}{3\Gamma\xi}\frac{\partial(\xi
v_{r})}{\partial\xi}+\frac{1}{3}\frac{\partial
v_{z}}{\partial\eta}-\frac{\epsilon\lambda}{\Gamma^{2}\xi}\partialderivative{\xi}\left[\xi\kappa(\phi)\partialderivative{p}{\xi}\right]-\frac{\lambda}{3}\partialderivative{\eta}\left[\kappa(\phi)\partialderivative{\sigma^{s}_{zz}}{\eta}\right]=0.$

Here $\lambda=a^{2}/(L^{2}\epsilon^{3})$ is a constant; we will later identify
it as the ratio of the fracture length scale to the development length scale,
over which we expect proppant flow to stabilise.

According to Shiozawa & McClure [22], Chen Zhixi et al. [23] and Liang et al.
[24], we utilise the following constants, relevant to hydraulic fracturing, as
given in Table 1.

Constant | Typical Value
---|---
$Q$ | $0.04\textrm{m}^{3}\textrm{ s}^{-1}$
$E^{\prime}$ | $40\textrm{ GPa}$
$\eta_{f}$ | $0.01\textrm{ Pa s}$
$\rho_{f}$ | $1000\textrm{ kg m}^{-3}$
$K_{Ic}$ | $0.5\textrm{ MPa m}^{0.5}$
$a$ | $5\times 10^{-5}\textrm{m}$
Table 1. Typical values of constants, given by Shiozawa & McClure [22], Chen
Zhixi et al. [23] and Liang et al. [24].

The choice of $a$ represents a typical diameter for the finer proppants
commonly used at the initiation of fracturing [24]. This gives us the
following estimates

$\displaystyle\epsilon$ $\displaystyle\approx 6\times 10^{-5}\cdot t^{-1/3},$
$\displaystyle L$ $\displaystyle\approx 9\times 10^{0}\cdot t^{4/9},$
$\displaystyle\textrm{Re}_{\textrm{eff}}$ $\displaystyle\approx 1\times
10^{-2}\cdot t^{-7/9},$ $\displaystyle\aleph$ $\displaystyle\approx 4\times
10^{-2}\cdot t^{1/9},$ $\displaystyle\lambda$ $\displaystyle\approx 1\times
10^{2}\cdot t^{1/9},$ $\displaystyle a/(\epsilon L)$ $\displaystyle\approx
1\times 10^{-1}\cdot t^{-1/9}.$

The value of $\textrm{Re}_{\textrm{eff}}$ is calculated using formula (10),
substituting each term with its typical scaling.

Considering the same problem in the zero-proppant case, Detournay & Savitski
[9] show that when $1.6\aleph<1$, the fracture evolution is well approximated
by taking the dimensionless toughness $\aleph=0$. Also, the choice $T=\aleph$
is taken, reflecting the dependence of the scaled solution on this
monotonically increasing parameter; assuming $\aleph$ is negligible it is
possible to neglect any $T$ dependence. We will also use these assumptions,
since toughness plays its greatest role near the fracture tip, where the crack
is typically too narrow for proppants to interfere. Given our estimate for
$\aleph$, this means we must take $t<1.5\times 10^{7}$.

In general we will assume $t>250$, so we may ignore $\epsilon$ and
$\textrm{Re}_{\textrm{eff}}$ terms. This also means $2a/(\epsilon L)<1/10$, so
the fracture is typically more than 10 particles wide. Lecampion & Garagash
[13], conclude that non-local phenomena such as proppant-bridging aren’t
important in such cases; however we can still expect to see these effects near
the narrow crack tip. The significance of this behaviour will be discussed in
greater detail in Section 7.

We also note that $\lambda$ is large; so in an effort to remove time
dependence from our equations, we may neglect the first three terms in the
continuity equation (22),

(23)
$\displaystyle\partialderivative{\eta}\left[\kappa(\phi)\partialderivative{\sigma^{s}_{zz}}{\eta}\right]=0.$

By the assumption of reflectional symmetry, the particle pressure gradient
must vanish at $\eta=0$. Because $\kappa$ is generally non-zero, we deduce
that the particle pressure is constant with $\eta$; and, by (19), so is
$|\eta|/\mu(\phi)$. Hence,

(24)
$\displaystyle\phi(\xi,\eta)=\mu^{-1}\left(\mu_{w}(\xi)\frac{|\eta|}{w(\xi)}\right),$

where $\mu_{w}$ is an undetermined function of $\xi$, which we recognise as
the value of $\mu$ at the crack wall. Noting that $\mu$ is a decreasing
function, we see that $\mu_{w}$ also describes the rate at which the
concentration drops from the centre to the wall of the cavity. We also notice
that, in accordance to Donstov et al. [25], we have plug flow in the centre of
the channel, where concentrations are greater than $\phi_{m}$. Because the
slurry flows away from the wellbore, the distribution of proppants, which is
described by $\mu_{w}$, depends on the concentration of proppants in the
injected mixture and how that changes with time. Hence, an important step in
the determination of $\mu_{w}$ will be implementing the assumption that the
average concentration at the wellbore is constant. This will be discussed in
greater detail in Section 7.

It is interesting to note that [13] verifies a length scale of
$\epsilon^{3}L^{3}/a^{2}$ for proppant flow in a channel, or pipe, to become
fully established. This means the particle pressure gradient becomes
negligible, and the cross fracture concentration profile becomes independent
of the distance from the channel, or pipe, entrance. As a result, the constant
$\lambda=a^{2}/(L^{2}\epsilon^{3})$ can be interpreted as the ratio of the
fracture length to the development length. Because this is large, an
alternative route to (24) would have been to assume the transverse particle
pressure is constant, reflecting the full development of the flow.

##  5\. Injection: Governing Equation for fracture width

In scaled coordinates, the governing equation for the conservation of proppant
mass becomes

(25)
$\displaystyle\frac{\xi\dot{L}t}{L}\frac{\partial\phi}{\partial\xi}+\left[\frac{\dot{\epsilon}t}{\epsilon}+\frac{\dot{L}t}{L}\right]\eta\frac{\partial\phi}{\partial\eta}=-\frac{\dot{\epsilon}t}{\epsilon\Gamma\xi}\partialderivative{(\xi\phi
v_{r})}{\xi}-\frac{\dot{\epsilon}t}{\epsilon}\partialderivative{(\phi
v_{z})}{\eta}.$

Then, implementing our choices of $\epsilon$ and $L$, we obtain

(26)
$\displaystyle\frac{4\xi}{3}\frac{\partial\phi}{\partial\xi}+\frac{\eta}{3}\frac{\partial\phi}{\partial\eta}=\frac{1}{\Gamma\xi}\partialderivative{(\xi\phi
v_{r})}{\xi}+\partialderivative{(\phi v_{z})}{\eta}.$

Integrating from $-w$ to $w$ with respect to $\eta$, leaving details to
Appendix A for brevity, we obtain

(27) $\displaystyle
4\xi\partialderivative{\xi}\left[w\Pi\circ\mu_{w}(\xi)\right]-w\Pi\circ\mu_{w}(\xi)=-\frac{9}{\Gamma^{2}\xi}\partialderivative{\xi}\left[\frac{\xi
w^{3}}{\mu_{w}(\xi)^{2}}\partialderivative{p}{\xi}\Omega\circ\mu_{w}(\xi)\right].$

Here we have defined the rheological functions

(28) $\displaystyle\Pi(x)$
$\displaystyle=\frac{1}{x}\int_{0}^{x}\mu^{-1}(u)du,$ $\displaystyle\Omega(x)$
$\displaystyle=\frac{1}{x}\int_{0}^{x}[\Pi(u)I\circ\mu^{-1}(u)u]du,$

which we plot in Figure 3.

(a) $\Pi(x)$ as a function of $x$

(b) $\Omega(x)$ as a function of $x$

(c) $x^{2}\Pi/\Omega$ as a function of $\Pi(x)$

Figure 3. Plots of the rheological functions $\Omega$, $\Pi$ and
$x^{2}\Pi/\Omega$.

Multiplying by $\xi$ and integrating from $\rho$ to $1$, we obtain

(29) $\displaystyle\int_{\rho}^{1}\xi
w\Pi\circ\mu_{w}(\xi)d\xi+\frac{4}{9}\rho^{2}w\Pi\circ\mu_{w}(\rho)=-\frac{\rho
w^{3}}{\Gamma^{2}\mu_{w}^{2}}\partialderivative{p}{\rho}\Omega\circ\mu_{w}(\rho),$

which lends itself more easily to computation. Here we have taken
$w^{3}\partial p/\partial\xi\to 0$ as $\xi\to 1$; this is physically motivated
by the fact that this term is proportional to the radial flux, which vanishes
at the crack tip. Moreover, Spence & Sharp [8] show that, in the zero-
proppant, zero-toughness regime, near the crack tip, $p\propto(1-\xi)^{-1/3}$
and $w\propto(1-\xi)^{2/3}$.

In order to compare this equation to the zero-proppant case, we assume
$\mu_{w}$ is independent of $\xi$ and take $\mu_{w}\to\infty$, to obtain

(30) $\displaystyle\int_{\rho}^{1}\xi
w(\xi)d\xi+\frac{4}{9}\rho^{2}w=-\frac{\rho
w^{3}}{\Gamma^{2}}\partialderivative{p}{\rho}\lim_{\mu_{w}\to\infty}\left[\frac{\Omega(\mu_{w})}{\mu_{w}^{2}\Pi(\mu_{w})}\right].$

From Figure 3(c) we deduce the right hand limit is approximately $2/5$, which
is confirmed exactly in Appendix B. Modelling the fluid as Newtonian, also
leaving the details to Appendix B, we obtain the same equation, with a factor
of $1/3$ instead. We conclude that the equations governing Newtonian flow are
not the same as those in the zero-proppant slurry flow limit. This is clearly
a limitation of our approach, which arises from using a dense-fitted rheology
in the dilute regime. However, the fact that the equations share a nearly
identical form is promising, as we expect the qualitative behaviour of slurry
flow to be similar to that of Newtonian flow.

##  6\. Injection: Numerical Solution

We implement the numerical method first used by Spence & Sharp [8], with the
adaptions of Detournay & Savitski [9], to solve the equations we have derived
so far. It will be useful to introduce $h(\xi)=w(\xi)/\Gamma$. The lubrication
equation derived above, the elasticity equations and the global volume
conservation equation become

(31) $\displaystyle\int_{\rho}^{1}(\xi h\Pi\circ\mu_{w})d\xi$
$\displaystyle+\frac{4}{9}\rho^{2}h\Pi\circ\mu_{w}=-\rho
h^{3}\partialderivative{p}{\rho}\frac{\Omega\circ\mu_{w}}{\mu_{w}^{2}},$ (32)
$\displaystyle h(\xi)$
$\displaystyle=\frac{4}{\pi}\int_{\xi}^{1}\frac{y}{\sqrt{y^{2}-\xi^{2}}}\int_{0}^{1}\frac{xp(xy)}{\sqrt{1-x^{2}}}dxdy,$
(33) $\displaystyle 0$
$\displaystyle=\int_{0}^{1}\frac{p(\xi)\xi}{\sqrt{1-\xi^{2}}}d\xi,$ (34)
$\displaystyle 1$ $\displaystyle=4\pi\Gamma^{3}\int_{0}^{1}(\xi h)d\xi.$

These equations alone do not give unique solutions for $\\{p,h,\mu_{w}\\}$, so
we will prescribe $\mu_{w}$ as part of the problem data. This allows us to
uniquely determine a solution for $\\{p,h\\}$. We seek series approximations
of the form

(35) $\displaystyle p(\xi)$ $\displaystyle=\sum_{i=-1}^{N-1}A_{i}p_{i}(\xi),$
$\displaystyle h(\xi)$ $\displaystyle=\sum_{i=-1}^{N}B_{i}h_{i}(\xi),$

where we define

(42) $\displaystyle p_{i}(\xi)$
$\displaystyle=\left\\{\begin{array}[]{ll}-\ln\xi+\ln 2-1&(i=-1)\\\ &\\\
(1-\xi)^{-1/3}J_{i}(\frac{4}{3},2,\xi)+\omega_{i}&(i\geq
0)\end{array}\right\\},$ $\displaystyle
h_{i}(\xi)=\left\\{\begin{array}[]{ll}\frac{4}{\pi}\left[(1-\xi^{2})^{1/2}-\xi\cos^{-1}(\xi)\right]&(i=-1)\\\
&\\\ \\\ (1-\xi)^{2/3}J_{i}(\frac{10}{3},2,\xi)&(i\geq
0)\end{array}\right\\}.$

Here the $i=-1$ terms are used to account for the logarithmic singularity in
pressure at the inlet, expected as a result of the point source injection; the
other terms allow for a general solution of (32). Importantly, we note that
the $p_{i}$ terms have a $(1-\xi)^{-1/3}$ singularity near the crack tip and
the $h_{i}$ terms are proportional to $(1-\xi)^{2/3}$ (for $i\geq 0$). This
deliberately matches the asymptotic calculations from Spence & Sharp [8],
which arise from the assumptions of zero-lag and zero-toughness in an
expanding hydraulic fracture. This allows the numerical method to converge
accurately with few terms. The $J_{i}(p,q,\xi)$ are Jacobi Polynomials of
order $i$ defined on the interval $[0,1]$, in the sense defined by Abramowitz
& Stegun [26], normalised to satisfy the orthonormality condition,

(43)
$\displaystyle\int_{0}^{1}(1-\xi)^{p-q}\xi^{q-1}J_{i}(p,q,\xi)J_{j}(p,q,\xi)d\xi=\delta_{ij}.$

This means that the $h_{i}$ ($i\geq 0$) are orthonormal with respect to an
inner product weighted by $\xi$. The $\omega_{i}$ are simply constants to
ensure each of the $p_{i}$ obey the zero-toughness equation; adding these
constants means that the $p_{i}$ lose their orthonormality properties, however
this doesn’t affect the solution finding process.

Because of its linearity, these series approximations reduce (32) to a linear
equation,

(44) $\displaystyle B_{i}=\sum_{j=-1}^{N-1}P_{ij}A_{j}.$

Here $(P)_{ij}$ is an $(N+2)\times(N+1)$ matrix whose entries we only have to
calculate once by using the orthogonality relation given above, along with the
fact that $\\{p_{-1},\theta_{-1}\\}$ are a solution pair to (32). The entries
of $M$, which can be found in [9], are listed in Appendix C for $N=4$. The
subtleties of calculating elements of $P_{ij}$, in the face of strong singular
behaviour, are important and described in depth in [9]. Finally, using the
values of $B_{i}$ given above, we assign a cost to each choice of $A$ given by

(45)
$\displaystyle\Delta(A)=\sum_{\xi\in\\{0,1/M,...,1\\}}\left(\frac{\textrm{RHS}(\xi;A)}{\textrm{LHS}(\xi;A)}-1\right)^{2}.$

This is calculated by considering the discrepancies between the left and right
hand sides of (31), calculated at M+1 equally spaced control points. We then
minimise $\Delta$ with respect to $A$ using the Nelder-Mead Simplex method
[27].

##  7\. Injection: Solutions for a constant $\mu_{w}$

For most monotonic choices of $\mu_{w}$, the numerical method above shows good
convergence. We see that the coefficients $A_{i}$ and $B_{i}$ drop off quickly
with $i$, and the final value of $\Delta$ tends to zero rapidly as we increase
$N$. If $\mu_{w}$ is a more complicated function, like in the case of Figure
4, we may need to use a larger value of $N$, but good convergence is still
possible.

Figure 4. Plot of cavity width profile and proppant distribution in the case
where $\mu_{w}$ is sinusoidal. Here $N=8$ is used.

This leads us to consider which choices of $\mu_{w}$ are most likely to appear
in reality. We note that by (24),

(46)
$\displaystyle\Pi\circ\mu_{w}(\xi)=\frac{1}{2w}\int_{-w}^{w}\phi(\xi,\eta)d\eta,$

so we may view $\Pi\circ\mu_{w}(\xi)$ as the average proppant concentration at
a given value of $\xi$. Since $\Pi\circ\mu_{w}$ is independent of time, we
automatically satisfy the condition that the injection rates of the proppants
and the fluid are constant. However this condition also means that the average
concentration at the wellbore, $\Pi\circ\mu_{w}(0)$, must equal the average
concentration taken by integrating over the entire crack volume. For a
monotonic choice of $\mu_{w}$ this implies that $\mu_{w}$ must be independent
of $\xi$. Herein we will make the assumption that $\mu_{w}$ is a constant and,
as a result, so is $\Pi=\Pi(\mu_{w})$. This is a natural assumption: at early
times we don’t expect significant concentration differences along the crack
because radial length scales are small.

A great advantage of a constant $\Pi$ is that we can define an ‘effective
viscosity’, which we can absorb into our scaled variables the same way as we
did with fluid viscosity. Under the assumption that $\mu_{w}$ is constant,
(31) becomes

(47) $\displaystyle\int_{\rho}^{1}\xi
h(\xi)d\xi+\frac{4}{9}\rho^{2}h=-\frac{\rho
h^{3}}{\eta_{e}}\partialderivative{p}{\rho},$

where $\eta_{e}=\mu_{w}^{2}\Pi/\Omega$ is what we call the effective
viscosity. It is plotted in Figure 3(c), and is best thought of as a function
of the average concentration, $\Pi$. Making the transformations

(48) $\displaystyle h$ $\displaystyle=\eta_{e}^{1/3}\tilde{h},$ $\displaystyle
p$ $\displaystyle=\eta_{e}^{1/3}\tilde{p},$ $\displaystyle\Gamma$
$\displaystyle=\eta_{e}^{-1/9}\tilde{\Gamma},$

our governing equations become

(49) $\displaystyle\int_{\rho}^{1}\xi\tilde{h}d\xi$
$\displaystyle+\frac{4}{9}\rho^{2}\tilde{h}=-\rho\tilde{h}^{3}\partialderivative{p}{\rho},$
$\displaystyle\tilde{h}(\xi)$
$\displaystyle=\frac{4}{\pi}\int_{\xi}^{1}\frac{y}{\sqrt{y^{2}-\xi^{2}}}\int_{0}^{1}\frac{x\tilde{p}(xy)}{\sqrt{1-x^{2}}}dxdy,$
$\displaystyle 0$
$\displaystyle=\int_{0}^{1}\frac{\tilde{p}(\xi)\xi}{\sqrt{1-\xi^{2}}}d\xi,$
$\displaystyle 1$
$\displaystyle=4\pi\tilde{\Gamma}^{3}\int_{0}^{1}(\xi\tilde{h})d\xi.$

We will solve them using the numerical method described before, except with
(49) in the place of (31-34).

Figure 5 plots $\tilde{h}$ and $\tilde{p}$, calculated using $N=4$ and
$M+1=501$ control points. Promisingly, we note that $\tilde{h}>0$ and $p$
shows the expected asymptotic behaviour. The value $\tilde{h}(0)=1.36$ will be
important in later discussion. The first column of table 3 shows the
coefficients $A_{i}$ and $B_{i}$, as well as the calculated value of
$\tilde{\Gamma}=0.598$. Significantly, we see that $A_{i}$ and $B_{i}$
decrease rapidly with $i$, suggesting that a solution with higher order terms
is unnecessary. This is supported by the small value of $\Delta\approx 5\times
10^{-5}$, with evenly spread contributions from control points along the
radius of the crack. This suggests that we have found a genuine solution, and
that the tip asymptotics are indeed suitable.

Figure 5. $(\xi,\eta)$ plots of $\tilde{h}$ and $\tilde{p}$, the scaled width
and pressure solutions to the absorbed effective viscosity system.

We now focus on finding numerical solutions for different concentrations in
order to consider features such as the velocity profile and proppant
distribution within the cavity. We consider the case of four different values
of the average concentration, $\Pi$. These are given in table 2, along with
the corresponding values of $\mu_{w}$ and $\eta_{e}$.

$\Pi$ | $\mu_{w}$ | $\eta_{e}$
---|---|---
0.05 | 487.3 | 2.74
0.20 | 23.35 | 3.92
0.40 | 3.93 | 10.37
0.55 | 1.06 | 96.60
Table 2. Test values of $\Pi$, $\mu_{w}$ and $\eta_{e}$.

The latter columns of table 3 show the values of $A$, $B$ and $\Gamma$
calculated using the exact method suggested in Section 6. Again we use
$M+1=501$ control points and $N=4$. Happily, the same values are observed by
using the values of $A$, $B$ and $\Gamma$ listed in the first column,
calculated after absorbing the effective viscosity, and using the relations
(48) to return to the concentration-specific values. We calculate the same
value of $\Delta\approx 5\times 10^{-5}$ each time; this is to be expected as
the equations are equivalent once the solutions have been scaled.

| $\Pi$ |  | Absorbed | 0.05 | 0.20 | 0.40 | 0.55
---|---|---|---|---|---|---|---
| $A_{-1}$ |  | 0.14786 | 0.20710 | 0.23326 | 0.32238 | 0.67830
| $A_{0}$ |  | 0.53529 | 0.74974 | 0.84444 | 1.16709 | 2.45559
| $A_{1}$ |  | 0.01929 | 0.02702 | 0.03043 | 0.04206 | 0.08849
| $A_{2}$ |  | 0.00402 | 0.00563 | 0.00634 | 0.00877 | 0.01844
| $A_{3}$ |  | 0.00035 | 0.00049 | 0.00055 | 0.00076 | 0.00159
| $B_{-1}$ |  | 0.14786 | 0.20710 | 0.23326 | 0.32238 | 0.67830
| $B_{0}$ |  | 0.53805 | 0.75361 | 0.84879 | 1.17311 | 2.46825
| $B_{1}$ |  | 0.05435 | 0.07612 | 0.08573 | 0.11849 | 0.24931
| $B_{2}$ |  | 0.00012 | 0.00016 | 0.00019 | 0.00026 | 0.00054
| $B_{3}$ |  | 0.00081 | 0.00114 | 0.00128 | 0.00177 | 0.00373
| $B_{4}$ |  | 0.00029 | 0.00041 | 0.00046 | 0.00064 | 0.00134
| $\Gamma$ |  | 0.59812 | 0.534579 | 0.513799 | 0.461261 | 0.359968
Table 3. Values of $A_{i}$, $B_{i}$ and $\Gamma$ obtained using (49) with
effective viscosity absorbed into the scaling and (31-34) with
$\Pi\in\\{0.05,0.20,0.40,0.55\\}$. We use $M=500$ and $N=4$ throughout.

Figure 6 shows the distribution of proppants within the fracture for each
value of $\Pi$. They are overlaid with an arrow plot of the proppant velocity
profile, $\mathbf{v}$, scaled by $\xi$ to show the equivalent two-dimensional
flux. The calculation of $\mathbf{v}$ is omitted since it is lengthy and
similar to the derivation of (27) in Appendix A. As $\Pi$ increases we see a
growing disk of plug flow where $\phi>\phi_{m}$, marked with a magenta
contour. We also see a tendency towards proppant velocity across the crack,
rather than along it; this is because the shape of the crack becomes shorter
and wider as the effective viscosity increases.

(a) $\Pi=0.05$

(b) $\Pi=0.20$

(c) $\Pi=0.40$

(d) $\Pi=0.55$

Figure 6. Concentration-specific $(\Gamma\xi,\eta)$ plots of developing
fractures with total solid volume fraction, $\Pi$, taking the values $0.05$,
$0.20$, $0.40$ and $0.55$. These are presented with filled contours displaying
proppant concentration; arrows showing $\xi$-scaled velocity; and magenta
contours indicating the transition into plug flow at the centre of each
cavity.

Drawing on calculations we have made so far, we are now in a position to
assess the significance of tip screen-out in our model, something we have
neglected so far by adopting a continuum model of proppant transport. This is
where, near the crack tip, the narrowing crack aperture causes proppants to
jam and block the fracture, significantly affecting the development of the
evolving formation and the convective transport of proppants. In [28] this
problem is addressed using a ‘blocking function’ which reduces proppant flux
to zero in apertures smaller than three times the average particle’s diameter.
We will use this threshold to weigh the significance of ignoring screen-out in
our model. Figure 7(a) shows the volume-proportion of proppants predicted in
fracture regions of width less than this threshold, dependant on the time,
$t$, and the average proppant concentration, $\Pi$. We see that for early
times and low concentrations, our model predicts a significant proportion of
proppants in these regions, where the fracturing fluid is clear in reality.
However, in concentrations greater than $0.3$ this proportion is relatively
small; this means our model, which ignores tip screen-out, is self-consistent.
This difference arises from the effective viscosity, which increases with
$\Pi$ and causes the ratio of fracture width to length to decrease.

Lecampion & Garagash [13] conclude that their rheology, which is employed
throughout this paper, agrees very well with experimental results when the
predicted width of plug flow is greater than a particle’s width. In figure
7(b), we see this condition holds for moderate times when $\phi>0.4$. It does
not for $\phi<0.4$. Therefore, in this regime we can expect slight mismatches
between predicted and practical concentration profiles; this arises from a
breakdown of the continuum model in the jammed part of the flow [13].

(a) $w<6a$

(b) Plug width $<2a$

Figure 7. Proportion of proppants by volume, predicted in fracture regions
where $w<6a$, or plug width $<2a$, given average concentration, $\Pi$, and
time, $t$.

##  8\. Crack Closure: Problem Formulation

In the zero-proppant case, Lai et al [29] have confirmed experimentally that
for late times after the fluid pressure is released, the crack radius is
constant and volume scales as $t^{-1/3}$. It is tempting to repeat our
previous work in order to find an asymptotic solution with a generalised total
fracture volume $Qt^{\alpha}$. We would then let $\alpha=-1/3$ to model the
case of closure. This approach leads us to

(50) $\displaystyle\alpha\int_{\rho}^{1}\xi
h(\xi)d\xi+\beta\rho^{2}h=-\frac{\rho
h^{3}}{\eta_{e}}\partialderivative{p}{\rho},$

in the place of (47). Here $\beta=(3\alpha+1)/9$ is the exponent for $L$,
giving the radial growth of the fracture. However, we see that attempts to
solve (50) using the previous numerical method fail as
$(\alpha,\beta)\to(-1/3,0)$, corresponding to the case in [29]. This is
because the tip asymptotes $w\propto(1-\xi)^{2/3}$ and
$p\propto(1-\xi)^{-1/3}$ are a result of an advancing fracture in a zero-
toughness medium. Spence & Sharp [8] note that $h\sim C(1-\xi)^{\tau}$ implies
$p\sim C\tau(\cot\pi\tau)(1-\xi)^{\tau-1}$. Balancing terms in (50), we are
forced with $C\leq 0$ if $\beta\leq 0$ which clearly can’t lead to physical
solutions, given the constraint $h\geq 0$. In the same paper, solutions for
$\beta=0$ are shown to exist without the assumption of zero-toughness; these
have $h\sim(1-\xi^{2})^{1/2}$. However, this causes difficulties in the case
of an evolving fracture, since a non-zero toughness parameter, $\aleph$,
brings time dependence to the scaled equations we have derived. An alternative
solution would be the addition of a non-zero fluid lag, providing a region of
negative pressure between the fluid front and the crack tip. Such a region
exists in reality, containing either vapour from the fracturing fluid or, if
the surrounding medium is permeable, pore fluid [30, 31]. Zero-toughness
solutions using this formulation are explored in [32]. Schematics of each
possible solution type are shown in Figure 8.

Figure 8. Possibilities for modelling the crack tip.

Any model utilising a time independent concentration profile is likely to fail
in describing fracture closure at late times. This is because the width of the
crack is decreasing as $t^{-1/3}$, so it is bound to become comparable to the
proppant diameter. At the point where $\epsilon L/a\approx 6$, the proppants
begin to bridge across the fracture, effectively fixing them in position [28];
therein, concentrations will increase as the carrier fluid is forced from the
cavity. For this reason, we will instead address the problem of finding the
residual crack shape, given some axisymmetric initial distribution of
proppants; we will assume these are radially immobile from the moment pressure
is released. This method has been used with success to model the closure of a
bi-wing fracture by Wang et al. [33, 34].

##  9\. Crack Closure: Residual Width Profiles

We model the residual shape of the fracture using $w_{p}(r)$, defined as the
close packed width of proppants. That is to say, after packing the proppants
as tightly as possible in the z direction, so $\phi=\phi_{rcp}$, this is the
residual width. Given some radial distribution of proppants described by the
average concentration, $\Pi$, and un-scaled width profile, $w$, we deduce that
$w_{p}=w\Pi/\phi_{rcp}$. This description is compatible with the frictional
rheology of Lecampion & Garagash [13], used previously, which asserts that a
non-zero normal force on the proppants, along with vanishing shear stress,
causes compression up to the random close packing limit. We then assume that
the surrounding fracture simply collapses around the proppant pack. Our
primary interest will be in using proppant distributions, arising from the
injection phase described previously, to predict the geometry of the residual
formation.

In [34] a more complicated model is offered; this considers stress from the
contact of opposing crack asperities, proppant embedment into the fracture
walls, and compression of proppants. Since we will be concerned with cases
where $w_{p}$ is non-zero along the entire crack radius; the contact term
arising from the crack asperities, which is significant in the un-propped
case, will not be necessary. Furthermore, in the same paper [34] the depth of
proppant embedment is shown to be of the order
$K_{e}=a(3/4E^{\prime})^{2}(16mE^{\prime 2}/9c_{p})^{2/3}$. Here, $m\approx
2\sqrt{3}$ is a constant which depends on the packing of proppants. Using the
value of $c_{p}=3.9\times 10^{-8}\textrm{Pa}^{-1}$ [34], as well as the
typical values of $a=50\mu\textrm{m}$ and $E^{\prime}=40\textrm{GPa}$
mentioned earlier, we note that $K_{e}\approx 1\mu\textrm{m}$, around 100
times smaller than the given proppant diameter. Since we will generally model
proppant packs which are several times the size of the proppant diameter in
width, we will ignore this phenomenon. Finally, we note that, according to our
previous estimates, more than $10\textrm{s}$ into the injection phase we
should expect pressures of less than $1\textrm{MPa}$. In [34] the compressive
stress required to reduce the width of the closely packed proppant bed from
$w_{p}$ to $w$ is given by $1/c_{p}\ln(w_{p}/w)$; using this, the same stress
would only cause a $4\%$ reduction in width. Since typical stresses involved
in the closure phase are much smaller than this, we will model the proppants
as incompressible.

This model of crack closure leads to a simple description of the residual
crack profile. We have two parameters: one for average concentration, $\Pi$,
and another for the time that injection ceases, $t_{0}$. Herein we will denote
$\\{\tilde{h},\tilde{p},\tilde{\Gamma}\\}$ as the solution to the system of
equations given in (49); $\tilde{h}$ and $\tilde{p}$ are plotted in Figure 5
and we use the value $\tilde{\Gamma}=0.598$. Then, using (48) and the original
scaling arguments, we deduce that

(51) $\displaystyle w_{p}(\xi;t_{0},\Pi)$
$\displaystyle=\frac{\Pi}{\phi_{rcp}}\epsilon(t_{0})L(t_{0})\eta_{e}(\Pi)^{2/9}\tilde{\Gamma}\tilde{h}(\xi),$
(52) $\displaystyle R(t_{0},\Pi)$
$\displaystyle=L(t_{0})\eta_{e}(\Pi)^{-1/9}\tilde{\Gamma}.$

From Figure 5 we notice that $\max(\tilde{h}_{1})\approx 1.35$. Using this, we
may plot Figure 9(a), which shows the effect of average concentration on the
maximum residual width of the formation. It is interesting to note that the
propped width doesn’t grow proportional to the proppant concentration, as one
may expect from the close packing of the suspended proppants. Instead, the
dependance is superlinear, because greater proppant concentrations lead to a
higher effective viscosity; this causes the fracture to take a wider shape
before the release of injection pressure. We can also see that $t_{0}$ has
relatively little effect on the maximum crack width. This is because the
$t_{0}$ dependent term, $\epsilon L$, grows with $t_{0}^{1/9}$. By contrast,
in Figure 9(b) we see a greater time dependence in the final radius, which
grows with $L\propto t^{4/9}$. As the proppant concentration increases, with
$t_{0}$ fixed, we see a decrease in the final radius of fracture achieved,
arising from an increase in the effective viscosity.

(a) Maximum fracture width.

(b) Fracture radius.

Figure 9. Plots showing the effect of average concentration on the maximum
residual fracture width and radius for $t_{0}\in\\{100,500,1000\\}$.

##  10\. Conclusions

We have established a mathematical framework that captures the behaviour of a
slurry within a pressure driven cavity. Using typical parameters from
industrial fracking, we predict that the development length, required to
establish stable proppant flow away from the wellbore, is negligible compared
to the typical radius of the penny-shaped fracture generated. As a result, we
may assume the flow is fully developed, reducing the in-fracture distribution
of proppants to a function of the radial distance from the wellbore. A further
assumption of constant proppant injection rate allows us to describe the
proppant distribution with one parameter, the total solid volume fraction. In
the zero-concentration limit, our model becomes similar to one derived using
Newtonian flow, with some disagreement arising from our choice of a dense
frictional rheology.

Within this framework, we are able to define an effective viscosity, which we
may absorb into our equations using a suitable choice of scaling. This is a
particularly striking result because it establishes an equivalence between
slurry flow of a given solid fraction and simple Newtonian flow with some
particular viscosity, at least in the sense of fracture development. Solving
the resulting set of equations numerically, we may then return to our original
scaling to investigate concentration-specific solutions. Unsurprisingly, we
predict width and pressure profiles with the tip-asymptotic behaviour
described in [9]. As the proppant concentration increases we expect shorter
and wider fractures with steeper fluid pressure gradients. In the centre of
the fracture, where shear rate vanishes, we predict the formation of a disk of
plug flow with width, in relation to the crack, increasing with the average
proppant concentration. Evaluating our model, we see that the unaccounted
effect of tip screen-out is likely to be significant in the low concentration,
low effective viscosity case, particularly at early times. Here, the cavity
formed is narrow, so near its tip, particle bridging is likely. Moreover, we
observe that for typical fracturing timescales, if $\Pi<0.4$, our model
predicts plug flow thinner than one particle width: suggesting that our use of
a continuum model may not be appropriate. Otherwise, the plug flow is broader
than a particle’s width, meaning it is physically realisable and the results
of [13] suggest we should have good experimental agreement.

Lastly, we have adopted a simple model of crack closure which regards the
remaining proppants to be immobile and incompressible. This allows us to
predict the shape of the residual crack, based on two parameters: the average
proppant concentration within the injected fluid and the length of time
between the initiation of fracking and the release of pressure. Simple
formulae show that the residual fracture width increases significantly with
proppant concentration, and grows very slowly with time; fracture radius
however, decreases with proppant concentration and increases with time.

The results established here have important applications in both contexts of
industrial fracking and geological dike formation. Diagnostics of tip screen-
out and forecasts of residual fracture geometry are relevant to the formation
of conductive fractures, whilst predictions about the shape and particle
distribution of a slurry driven crack relate more to a cooling magma. The
discovery of an effective viscosity may also provide a foothold in
understanding slurry driven fractures, particularly given the bounty of
literature surrounding cracks generated by Newtonian fluid. In spite of all
this, experimental investigation is necessary to bolster the predictions we
have made. We hope this will form the basis of a second article, with
tentative title: ‘Proppant flow in a penny-shaped crack. Part II :
Experimental Investigation’.

##  11\. Acknowledgements

The authors would like to thank Derek Elsworth (Pennsylvania State
University), Elisabeth Guazzelli (Centre National de la Recherche
Scientifique) and Emmanuel Detournay (University of Minnesota) for their
support and guidance in the drafting of this paper; with special gratitude to
Elisabeth for providing the data used in Figure 2. We would also like to thank
John Willis (University of Cambridge) for his support in the publication of
the paper.

##  Appendix A Integrating the $\phi$,phionservation equation over the crack
width

In this Appendix we integrate equation (25) over $(-w,w)$ to yield (27); we
will take a term-by-term approach. First, we note that by (24),

(53) $\displaystyle\int_{-z}^{z}\phi(\xi,\eta)d\eta$
$\displaystyle=2\int_{0}^{z}\mu^{-1}\left(\mu_{w}(\xi)\frac{\eta}{w}\right)d\eta,$
(54) $\displaystyle=2z\Pi\left(\mu_{w}(\xi)\frac{z}{w}\right).$

Hence, we see that

(55) $\displaystyle\int_{-w}^{w}\partialderivative{\phi}{\xi}d\eta$
$\displaystyle=\partialderivative{\xi}\int_{-w}^{w}\phi
d\eta-2\phi(\xi,w)\partialderivative{w}{\xi},$ (56)
$\displaystyle=2\partialderivative{\xi}\left[w\Pi\circ\mu_{w}(\xi)\right]-2\phi(\xi,w)\partialderivative{w}{\xi}.$

Then, integrating by parts, we find

(57)
$\displaystyle\int_{-w}^{w}\eta\partialderivative{\phi}{\eta}d\eta=2\left[w\phi(\xi,w)-w\Pi\circ\mu_{w}(\xi)\right].$

Furthermore, utilising the expression of $v_{r}$ given in (20) and the
condition $v_{r}(\xi,\pm w)=0$ we determine

(58) $\displaystyle\int_{-w}^{w}\partialderivative{(\xi\phi v_{r})}{\xi}d\eta$
$\displaystyle=\partialderivative{\xi}\left[\xi\int_{-w}^{w}\phi
v_{r}d\eta\right],$ (59)
$\displaystyle=-\frac{6}{\Gamma}\partialderivative{\xi}\left[\xi\partialderivative{p}{\xi}\int_{0}^{w}\phi(\xi,\eta)\int_{\eta}^{w}\frac{I(\phi(\xi,z))z}{\mu(\phi(\xi,z))}dzd\eta\right],$
(60)
$\displaystyle=-\frac{6}{\Gamma}\partialderivative{\xi}\left[\xi\partialderivative{p}{\xi}\int_{0}^{w}\int_{0}^{z}\phi(\xi,\eta)\frac{I(\phi(\xi,z))z}{\mu(\phi(\xi,z))}d\eta
dz\right],$ (61)
$\displaystyle=-\frac{6}{\Gamma}\partialderivative{\xi}\left[\xi\partialderivative{p}{\xi}\int_{0}^{w}z^{2}\Pi\left(\frac{\mu_{w}z}{w}\right)\frac{I(\phi(\xi,z))}{\mu(\phi(\xi,z))}dz\right].$

However, by (24), $\mu(\phi(\xi,z))=\mu_{w}z/w$, so

(62) $\displaystyle\int_{-w}^{w}\partialderivative{(\xi\phi v_{r})}{\xi}d\eta$
$\displaystyle=-\frac{6}{\Gamma}\partialderivative{\xi}\left[\frac{w\xi}{\mu_{w}}\partialderivative{p}{\xi}\int_{0}^{w}z\Pi\left(\frac{\mu_{w}z}{w}\right)I\circ\mu^{-1}\left(\frac{\mu_{w}z}{w}\right)dz\right],$
(63) $\displaystyle=-\frac{6}{\Gamma}\partialderivative{\xi}\left[\frac{\xi
w^{3}}{\mu_{w}(\xi)^{2}}\partialderivative{p}{\xi}\Omega\circ\mu_{w}(\xi)\right].$

Finally, we know that

(64) $\displaystyle\int_{-w}^{w}\partialderivative{(\phi
v_{z})}{\eta}d\eta=2\phi(\xi,w)v_{z}(\xi,w).$

In the original scaling we have the boundary condition
$v_{z}(x,w)=\partialderivative{w}{t}(x,t)$; in the lubrication scaling this
becomes

(65) $\displaystyle-\dot{\epsilon}Lv_{z}(\xi,w)$
$\displaystyle=\left[\dot{\epsilon}L+\epsilon\dot{L}\right]w(\xi,T)-\epsilon
L\xi\left[\frac{\dot{L}}{L}+\frac{\Gamma^{\prime}\dot{T}}{\Gamma}\right]\partialderivative{w}{\xi}+\dot{T}\partialderivative{w}{T}.$

Hence,

(66) $\displaystyle
v_{z}(\xi,w)=\frac{w}{3}-\frac{4\xi}{3}\partialderivative{w}{\xi},$

and so

(67) $\displaystyle\int_{-w}^{w}\partialderivative{(\phi
v_{z})}{\eta}d\eta=2\phi(\xi,w)\left[\frac{w}{3}-\frac{4\xi}{3}\partialderivative{w}{\xi}\right].$

Adding these terms together and making various cancellations, we derive
equation (27).

##  Appendix B Zero-Concentration Limit

In this Appendix, we will compare the properties of equation (29) to the
equivalent zero-proppant equation. Modelling the flow as Newtonian instead, we
would have used the relation $\tau=\eta_{f}\dot{\gamma}$. In our choice of
scaling this becomes $\tau=\dot{\gamma}$. Hence (19.2) is replaced by

(68)
$\displaystyle\partialderivative{v_{r}}{\eta}=\frac{3\eta}{\Gamma}\partialderivative{p}{\xi},$

where $\mathbf{v}$ is the fluid velocity. With the assumption that
$\nabla\cdot v=0$, our scaled continuity equation is simply

(69) $\displaystyle\frac{1}{\Gamma\xi}\partialderivative{(\xi
v_{r})}{\xi}+\partialderivative{v_{z}}{\eta}=0.$

Integrating first over $(-w,w)$ as in Appendix A, making use of (66), (68) and
$\tau=\dot{\gamma}$, we obtain

(70)
$\displaystyle\frac{w}{3}-\frac{4\xi}{3}\partialderivative{w}{\xi}=\frac{1}{\xi\Gamma^{2}}\partialderivative{\xi}\left[\partialderivative{p}{\xi}\xi
w^{3}\right].$

Then, multiplying by $\xi$ and integrating from $\rho$ to 1, we use the
$w^{3}\partial p/\partial\xi\to 0$ limit employed to derive (29),

(71) $\displaystyle\int_{\rho}^{1}\xi wd\xi+\frac{4}{9}\rho^{2}w=-\frac{\rho
w^{3}}{3\Gamma^{2}}\partialderivative{p}{\rho}.$

In order to compare (29) and (71), we are required to find the limit of
$\Omega/(x^{2}\Pi)$ as $x\to\infty$. Explicitly we see that

(72) $\displaystyle\lim_{x\to\infty}\frac{\Omega(x)}{x^{2}\Pi(x)}$
$\displaystyle=\lim_{x\to\infty}\frac{1}{x^{3}\Pi(x)}\int_{0}^{x}\Pi(u)I\circ\mu^{-1}(u)udu,$
(73)
$\displaystyle=\lim_{x\to\infty}\int_{0}^{1}\frac{\Pi(vx)}{\Pi(x)}\cdot\frac{I\circ\mu^{-1}(vx)}{vx}\cdot
v^{2}dv,$ (74)
$\displaystyle=\int_{0}^{1}v^{2}\lim_{x\to\infty}\left[\frac{\Pi(vx)}{\Pi(x)}\right]dv,$
(75)
$\displaystyle=\int_{0}^{1}v\lim_{x\to\infty}\left[\frac{\int_{0}^{vx}\mu^{-1}(u)du}{\int_{0}^{x}\mu^{-1}(u)du}\right]dv,$
(76)
$\displaystyle=\int_{0}^{1}v^{2}\lim_{x\to\infty}\left[\frac{\mu^{-1}(vx)}{\mu^{-1}(x)}\right]dv,$
(77)
$\displaystyle=\int_{0}^{1}v^{2}\lim_{x\to\infty}\left[\frac{I^{-1}(vx)}{I^{-1}(x)}\right]dv,$
(78)
$\displaystyle=\int_{0}^{1}v^{2}\lim_{x\to\infty}\left[\frac{1+\sqrt{x}}{1+\sqrt{vx}}\right]dv,$
(79) $\displaystyle=\int_{0}^{1}v^{3/2}dv,$ (80) $\displaystyle=2/5.$

Here (74) and (77) arise from the fact $I(\phi)\sim\mu(\phi)$ as $\phi\to 0$,
because the fluid shear stress approaches the slurry shear stress. (76) comes
from L’Hôpital’s rule. We conclude that the equations governing Newtonian flow
are not the same as those in the zero-proppant slurry flow limit.

##  Appendix C Matrix $(P)_{ij}$, when $N=4$

The matrix $(P)_{ij}$ for $N=4$, as provided in [9], is given in table 4.

|  |  |  |  | j |  |
---|---|---|---|---|---|---|---
|  |  | -1 | 0 | 1 | 2 | 3
| -1 |  | 1.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000
| 0 |  | 0.0000 | 0.9560 | 1.2730 | 0.4101 | 0.3145
i | 1 |  | 0.0000 | 0.0991 | -0.0185 | 0.4068 | 0.0610
| 2 |  | 0.0000 | 0.0018 | -0.0429 | -0.0244 | 0.2293
| 3 |  | 0.0000 | 0.0017 | 0.0039 | -0.0416 | -0.0141
| 4 |  | 0.0000 | 0.0005 | 0.0026 | -0.0032 | -0.0372
Table 4. Matrix $(P)_{ij}$, for N=4.

## References

  * [1] Wells, Bruce A., ed. (2007). ”Shooters”. The Petroleum Age. American Oil and Gas Historical Society. 4 (3): 8–9. ISSN 1930-5915
  * [2] Charlez, Philippe A. (1997). Rock Mechanics: Petroleum Applications. Paris: Editions Technip. p. 239. ISBN 9782710805861.
  * [3] National Earthquake Hazards Reduction Program (U.S.), Geological Survey (U.S.), Office of Earthquakes, Volcanoes, and Engineering, U.S. National Committee for Rock Mechanics (1983). Hydraulic Fracturing Stress Measurements. Volume 26 of International journal of rock mechanics and mining sciences and geomechanics abstracts.
  * [4] Pierce, Brenda (2010). Geothermal Energy Resources. National Association of Regulatory Utility Commissioners (NARUC).
  * [5] Miller, Bruce G. (2005). Coal Energy Systems. Sustainable World Series. Academic Press. p. 380. ISBN 9780124974517.
  * [6] E. Rivalta, B. Taisne, A.P. Bunger, R.F. Katz (2015). A review of mechanical models of dike propagation: Schools of thought, results and future directions. Tectonophysics. Volume 638,2015. Pages 1-42. ISSN 0040-1951.
  * [7] Petford, N., Koenders, M.A. (1998). Granular flow and viscous fluctuations in low Bagnold number granitic magmas. Journal of the Geological Society, 155 (5), pp. 873-881. 10.1144/gsjgs.155.5.0873
  * [8] Spence, D.A., Sharp, P.W. (1985). Self-similar solution for elastohydrodynamic cavity flow. Proc. Roy. Soc. London, Ser. A (400),289–313.
  * [9] A.A. Savitski, E. Detournay (2002). Propagation of a penny-shaped fluid-driven fracture in an impermeable rock: asymptotic solutions, International Journal of Solids and Structures, Volume 39, Issue 26, Pages 6311-6337.
  * [10] Einstein, A. (1906). A new determination of molecular dimensions. Ann. Phys. 4 (19), 289–306.
  * [11] Boyer F., Guazzelli É., Pouliquen O. (2011). Unifying suspension and granular rheology. Phys. Rev. Lett. 107 (18), 188301.
  * [12] Dontsov EV, Boronin SA, Osiptsov AA, Derbyshev DY. (2019). Lubrication model of suspension flow in a hydraulic fracture with frictional rheology for shear-induced migration and jamming. Proc. R. Soc. A 475: 20190039.
  * [13] Lecampion, Garagash (2014). Confined flow of suspensions modelled by a frictional rheology. J. Fluid Mech. (2014), vol. 759, pp. 197–235. Cambridge University Press 2014. doi:10.1017/jfm.2014.557
  * [14] Niall J. O’Keeffe, Herbert E. Huppert & P. F. Linden (2018). Experimental exploration of fluid-driven cracks in brittle hydrogels. J. Fluid Mech., vol. 844, pp. 435–458.
  * [15] Garagash, D.I., Detournay, E. (2000). The tip region of a fluid-driven fracture in an elastic medium. ASME J. Appl. Mech. 67, 183–192.
  * [16] Sneddon, I.N., (1951). Fourier Transforms. McGraw-Hill, New York, NY
  * [17] Rice, J.R., (1968). Mathematical analysis in the mechanics of fracture. In: Liebowitz, H. (Ed.), Fracture, an Advanced Treatise. Vol. II.Academic Press, New York, NY, pp. 191–311 (Chapter 3).
  * [18] Dagois-Bohy S., Hormozi S., Guazzelli É, Pouliquen O. (2015). Rheology of dense suspensions of non-colloidal spheres in yield-stress fluids. Journal of Fluid Mechanics, 776, R2. doi:10.1017/jfm.2015.329
  * [19] Garside, J., Al-Dibouni, M. R. (1977). Velocity-voidage relationships for fluidization and sedimentation in solid–liquid systems. Ind. Eng. Chem. Process Des. Dev. 16 (2), 206–214.
  * [20] Richardson, J., Zaki, W. (1954) Sedimentation and fluidization: Part I. Trans. Inst. Chem. Engrs 32, 35–47.
  * [21] Bacri, J.-C., Frenois, C., Hoyos, M., Perzynski, R., Rakotomalala, N. & Salin, D. (1986). Acoustic study of suspension sedimentation. Europhys. Lett. 2 (2), 123–128.
  * [22] Shiozawa, S., Mcclure, M. (2016). Simulation of proppant transport with gravitational settling and fracture closure in a three-dimensional hydraulic fracturing simulator. J. Petrol. Sci. Engng, 138, 298–314.
  * [23] Chen Zhixi, Chen Mian, Jin Yan, Huang Rongzun (1997). Determination of rock fracture toughness and its relationship with acoustic velocity, International Journal of Rock Mechanics and Mining Sciences, Volume 34, Issues 3–4, 1997, Pages 49.e1-49.e11, ISSN 1365-1609
  * [24] Feng Liang, Mohammed Sayed, Ghaithan A. Al-Muntasheri, Frank F. Chang, Leiming Li (2016). A comprehensive review on proppant technologies. Petroleum, Volume 2, Issue 1, March 2016, Pages 26-39.
  * [25] Dontsov EV, Boronin SA, Osiptsov AA, Derbyshev DY. (2019). Lubrication model of suspension flow in a hydraulic fracture with frictional rheology for shear-induced migration and jamming. Proc. R. Soc. A 475: 20190039.
  * [26] Abramowitz, M., Stegun, I.A. (Eds.), (1964). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.Applied Mathematics Series, 55. US Govt. Print. Off, Washington, DC.
  * [27] Lagarias, J. C., J. A. Reeds, M. H. Wright, & P. E. Wright (1998). Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions. SIAM Journal of Optimization. Vol. 9, Number 1, 1998, pp. 112–147.
  * [28] Dontsov, E. V., Peirce, A. P. (2014). Slurry flow, gravitational settling and a proppant transport model for hydraulic fractures. J. Fluid Mech. 760, 567–590.
  * [29] Ching-Yao Lai, Zhong Zheng, Emilie Dressaire, Guy Z. Ramon, Herbert E. Huppert, & Howard A. Stone (2016). Elastic Relaxation of Fluid-Driven Cracks and the Resulting Backflow. Physical Review Letters 117, 268001.
  * [30] A.M. Rubin. (1993). Tensile fracture of rock at high confining pressure: implications for dike propagation. J. Geophys. Res., 98 (B9) (1993), pp. 15,919-15,935.
  * [31] E. Detournay & D. Garagash (2003). The tip region of a fluid-driven fracture in a permeable elastic solid. J. Fluid Mech., 494, pp. 1-32.
  * [32] D. Garagash (2006). Propagation of a plane-strain hydraulic fracture with a fluid lag: Early-time solution, International Journal of Solids and Structures 43, 5811–5835.
  * [33] Jiehao Wang, Derek Elsworth & Martin K. Denison (2018). Propagation, proppant transport and the evolution of transport properties of hydraulic fractures. J. Fluid Mech., vol. 855, pp. 503–534.
  * [34] Wang, J. & Elsworth, D. (2018). Role of proppant distribution on the evolution of hydraulic fracture conductivity. J. Petrol. Sci. Engng 166, 249–262.