††thanks: Corresponding author<EMAIL_ADDRESS>
# Understanding the rheology of kaolinite clay suspensions using Bayesian
inference

Ranjiangshang Ran  Department of Mechanical Engineering and Applied Mechanics,
University of Pennsylvania, Philadelphia, PA 19104, USA Shravan Pradeep
Department of Earth and Environmental Science, University of Pennsylvania,
Philadelphia, PA 19104, USA Sébastien Kosgodagan Acharige  Department of
Mechanical Engineering and Applied Mechanics, University of Pennsylvania,
Philadelphia, PA 19104, USA Department of Earth and Environmental Science,
University of Pennsylvania, Philadelphia, PA 19104, USA Brendan C. Blackwell
Department of Mechanical Engineering and Applied Mechanics, University of
Pennsylvania, Philadelphia, PA 19104, USA Christoph Kammer  Department of
Mechanical Engineering and Applied Mechanics, University of Pennsylvania,
Philadelphia, PA 19104, USA Douglas J. Jerolmack  Department of Earth and
Environmental Science, University of Pennsylvania, Philadelphia, PA 19104, USA
Paulo E. Arratia  Department of Mechanical Engineering and Applied Mechanics,
University of Pennsylvania, Philadelphia, PA 19104, USA

###### Abstract

Mud is a suspension of fine-grained particles (sand, silt, and clay) in water.
The interaction of clay minerals in mud gives rise to complex rheological
behaviors, such as yield stress, thixotropy and viscoelasticity. Here, we
experimentally examine the flow behaviors of kaolinite clay suspensions, a
model mud, using steady shear rheometry. The flow curves exhibit both yield
stress and rheological hysteresis behaviors for various kaolinite volume
fractions ($\phi_{k}$). Further understanding of these behaviors requires
fitting to existing constitutive models, which is challenging due to numerous
fitting parameters. To this end, we employ a Bayesian inference method, Markov
chain Monte Carlo (MCMC), to fit the experimental flow curves to a
microstructural viscoelastic model. The method allows us to estimate the
rheological properties of the clay suspensions, such as viscosity, yield
stress, and relaxation time scales. The comparison of the inherent relaxation
time scales suggests that kaolinite clay suspensions are strongly viscoelastic
and weakly thixotropic at relatively low $\phi_{k}$, while being almost
inelastic and purely thixotropic at high $\phi_{k}$. Overall, our results
provide a framework for predictive model fitting to elucidate the rheological
behaviors of natural materials and other structured fluids.

##  I Introduction

Mud — a suspension of fine-grained particles (such as sand, silt, and clay)
and water — forms the foundation of Earth’s surface and plays important roles
in many natural (e.g., mudslides [1], soil erosion [2], river morphology [3,
4]) and industrial processes (e.g., civil construction [5], oil drilling [6]).
At low stresses, mud exhibits solid-like behavior, which gives the material
its load-bearing capacity. Above a certain threshold stress, however, these
suspensions can experience sudden mechanical failure and start to flow like a
liquid, as seen during the rapid and catastrophic mudslides [7]. This
threshold at which the transition from solid- to fluid-like behaviors occurs
is referred to as the “yield stress” [8, 9]. Ultimately, to understand
yielding behaviors in mud-like materials, we need models that can predict the
yield stress of mud suspensions based on its constituents and bulk rheology.
If successful, such rheological models could help improve failure predictions
for industrial processes and prevent hazardous natural events.

The rheology of mud is dominated by the attractive nature of clay minerals,
since these charged plate-shaped mud constituents are observed to aggregate
and form gel-like microstructure [10, 11, 12, 13, 14, 15]. This gives rise to
“thixotropy” [16, 17], a rheological behavior characterized by time-dependent
viscosity resulting from the evolution of the microstructure associated with
flow history [18]. Thixotropic behaviors often involve two processes: an
increase in viscosity due to the buildup of microstructure as the fluid is at
rest (or “aging”), and a decrease in viscosity due to the breakdown of
microstructure as the fluid is sheared (or “rejuvenation”) [19, 20].
Rheological models of thixotropy often use a structural parameter to relate
the fluid microstructure to its macroscopic properties (e.g., viscosity),
along with a kinetic equation to describe the evolution of the structural
parameter [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. A
similar structural parameter and accompanied kinetic equation are also used to
describe the rheological behaviors of “soft glassy” materials that are usually
not referred to as thixotropic [36, 37, 38, 39, 40, 41, 42, 43]. While the
structural parameter is often taken to be a single scalar quantity [23, 24,
25, 26, 27, 28, 29, 30, 31, 34, 35], it can also be represented as tensors
[22, 40, 43], scalar fields [41, 42], or a combination of multiple scalar
quantities [32, 33]. Despite its clear physical meaning, the estimation of the
structural parameter from experimental measurements is challenging.

Besides yield stress and thixotropy, mud and other clay-based suspensions can
also exhibit viscoelasticity [44, 45, 46, 47], since clay particles can form
gel-like microstructure that behaves like weak elastic solids [18, 48, 49].
Viscoelastic behaviors are characterized by the relaxation of elastic stress
over an observed time scale [50]. The time scales associated with stress
relaxation and microstructural evolution are usually different and
independent, but it is often difficult to distinguish them [18, 19, 20, 34].
Therefore, to capture the rheological behaviors of mud and other structured
fluids, such as food products [51, 52, 53, 54] and crude oils [55, 56],
constitutive models that incorporate both thixotropy and viscoelasticity are
necessary [24, 25, 26, 27, 28, 29, 30, 31, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 55, 56]. Such models commonly contain Maxwell- or Jeffrey-type
viscoelastic elements whose viscosity and elasticity depend on the structural
parameter. An exception is when the viscoelastic time scale (i.e., stress
relaxation) is much shorter than that of thixotropic time scale (i.e.,
microstructural evolution), in which case the material is called “ideally”
thixotropic and the effects of viscoelasticity can be neglected [32, 19, 20].

The interplay of thixotropy and viscoelasticity also lead to hysteresis
behavior [57, 58, 42, 59, 60]. A common procedure to identify such behavior is
by sweeping the shear rate (or shear stress) downward and upward over a
prescribed range [57, 42], which produces hysteresis loops in stress versus
shear rate flow curves. Large amplitude oscillatory shear (LAOS) rheology can
also identify hysteresis behavior [61, 62, 29, 30, 31, 63, 64]. Despite the
advances in rheological models and measurements, it remains challenging to fit
constitutive models to the experimental flow curves obtained from steady shear
or LAOS, due to the model complexity and numerous fitting parameters.
Therefore, common model fitting approaches often apply statistical methods
such as Bayesian inference [65] and neural networks [66, 67, 68, 69].

In this study, we experimentally examine the flow behavior of kaolinite clay
suspensions — a model mud — using shear rate sweep tests. The flow curves
exhibit both yield stress and hysteresis behaviors for different kaolinite
volume fractions ($\phi_{k}$). In a effort to understand these results, we fit
the experimental flow curves to a constitutive model that is both viscoelastic
and thixotropic, using Markov chain Monte Carlo (MCMC), a Bayesian inference
method. The model fitting results provide estimates of the rheological
properties of the clay suspensions, such as viscosity, yield stress, and
inherent relaxation time scales. A comparison of the relaxation time scales of
elastic stress and microstructure suggests that kaolinite clay suspensions are
strongly viscoelastic and weakly thixotropic at relatively low $\phi_{k}$,
while they are almost ideally thixotropic and inelastic at high $\phi_{k}$.
Our results also reveal that kaolinite clay suspensions possess multiple time
scales of stress and microstructural relaxation. Overall, our results provide
insights into the development of the rheological models of natural materials
such as mud and clay, as well as other structured fluids.

##  II Materials and Methods

###  II.1 Kaolinite clay suspensions

Kaolinite clay suspensions are prepared by dispersing natural kaolinite
particles (03584 Sigma-Aldrich) into deionized water at four different volume
fractions: $\phi=$ 0.05, 0.12, 0.19, and 0.26. Kaolinite is a layered aluminum
silicate ($\mathrm{Al}_{2}\mathrm{O}_{3}\cdot 2\mathrm{SiO}_{2}\cdot
2\mathrm{H}_{2}\mathrm{O}$), which consists of stacks of individual plate-
shaped particles (particle density, $\rho_{p}=2.61~{}\mathrm{g/cm}^{3}$)
joined together by hydrogen bonds. Kaolinite particles have a modal particle
size, $d\approx 10~{}\umu$m [70, 14], and an aspect ratio $d/h\approx 10$,
where $h$ is the average thickness of the particles. The zeta potential of
kaolinite particles at $\mathrm{pH}=7.0$ is estimated to be $-30$ mV [70],
which indicates that kaolinite particles are weakly attractive and tend to
aggregate in water.

###  II.2 Shear rate sweep test

A stress-controlled cone-plate rheometer (Bohlin, Malvern) is used to perform
shear rate sweep tests that are similar to the protocols developed by Divoux
et al. [57]. First, the sample is pre-sheared at a high shear rate
($\dot{\gamma}_{\mathrm{max}}$) for 180 s to eliminate previous flow history
and generate a reproducible initial microstructure. Next, two consecutive
shear rate sweeps are performed: a downward sweep from the highest shear rate
($\dot{\gamma}_{\mathrm{max}}=10^{3}~{}\mathrm{s}^{-1}$) to the lowest shear
rate ($\dot{\gamma}_{\mathrm{min}}=10^{-1}~{}\mathrm{s}^{-1}$) through $N=50$
logarithmically spaced steps of duration $\Delta t$ each, followed by an
upward sweep from $\dot{\gamma}_{\mathrm{min}}$ back to
$\dot{\gamma}_{\mathrm{max}}$ with the same $N$ logarithmic steps in a reverse
order. In this manuscript, we adopt the terms “upsweep” for the upward sweep
in shear rates
($\dot{\gamma}_{\mathrm{min}}\rightarrow\dot{\gamma}_{\mathrm{max}}$) and
“downsweep” for the downward sweep in shear rates
($\dot{\gamma}_{\mathrm{max}}\rightarrow\dot{\gamma}_{\mathrm{min}}$), as used
in previous literature [42]. Two different sweep time $\Delta t$ are used:
$\Delta t=10$ s and 20 s, leading to total experimental duration of $2N\Delta
t=2000$ s and 4000 s, respectively. Figure 1(a) shows the prescribed shear
rate as a function of time for $\Delta t=10~{}$s. The room temperature is
maintained at 21$~{}{}^{\circ}$C during all the tests.

Figure 1(b) shows the measured stress in a shear rate sweep test, for
kaolinite suspension of volume fraction $\phi=0.05$ and $\Delta t=10~{}$s. We
find that the suspension stress is asymmetric in time even though the applied
shear rate is symmetric in time 1(a); this asymmetry leads to rheological
hysteresis in the $\sigma$-$\dot{\gamma}$ flow curves. Figure 1(c) shows
examples of the rheological hysteresis behaviors for kaolinite suspensions of
different volume fractions: $\phi=$0.05, 0.12, and 0.19 with $\Delta
t=10~{}$s. At low $\dot{\gamma}$, the suspension stress in upsweep (red) lies
below the downsweep (blue) but beyond a critical shear rate
($\dot{\gamma_{c}}$) the flow curves reverse their behavior.

[width = 6.69 in]Fig1.png

Figure 1: (a) The prescribed shear rate ($\dot{\gamma}$) as a function of time
in the shear rate sweep test: a downsweep (blue) from
$\dot{\gamma}_{\mathrm{max}}=10^{3}~{}\mathrm{s}^{-1}$ to
$\dot{\gamma}_{\mathrm{min}}=10^{-1}~{}\mathrm{s}^{-1}$ in 50 logarithmic
steps of duration $\Delta t=10~{}$s is followed by an upsweep (red) from
$\dot{\gamma}_{\mathrm{min}}$ back to $\dot{\gamma}_{\mathrm{max}}$ with the
same steps in reverse order. (b) Experimental data for shear stress ($\sigma$)
as a function of time $t$, for kaolinite suspension of $\phi_{k}=0.05$ and
$\Delta t=10~{}$s, showing an asymmetry in time. (c) Flow curves $\sigma$ vs.
$\dot{\gamma}$ for kaolinite suspensions of different volume fractions:
$\phi_{k}=0.05$, 0.12, 0.19, and $\Delta t=10~{}$s, showing rheological
hysteresis behaviors.

###  II.3 Rheological model

A viscoelastic thixotropic model [29] is used here to capture the clay
suspensions’ rheological behavior described above. The model is represented by
a simple three-element Jeffreys canonical form in Fig. 2, where a constant
infinite shear viscosity ($\eta_{\infty}$) is modeled in parallel with a
viscoelastic Maxwell component that is dependent on a scalar structural
parameter, $\xi$. The viscous and elastic parts of the Maxwell component are
assumed to be linearly dependent on the structural parameter as $\xi\eta_{A}$
and $\xi G_{A}$, respectively, where $\eta_{A}$ is the structural viscosity
and $G_{A}$ is the elastic modulus. The time evolution of the shear stress
$\sigma$ under an applied shear rate $\dot{\gamma}$ in a Jeffreys-type fluid
can be described as:

$\sigma+\lambda\frac{d\sigma}{dt}=\left(\eta_{\infty}+\xi\eta_{A}\right)\left(\dot{\gamma}+\lambda_{r}\frac{d\dot{\gamma}}{dt}\right),$
(1)

where the fluid has a total viscosity $\eta=\eta_{\infty}+\xi\eta_{A}$, a
single relaxation time $\lambda=\eta_{A}/G_{A}$ independent of microstructure,
and a microstructural-dependent retardation time
$\lambda_{r}=\lambda\eta_{\infty}/(\eta_{\infty}+\xi\eta_{A})$. Equation (1)
can be further reduced to the following form [29, 30]:

$\frac{\sigma}{\eta_{A}}+\frac{1}{G_{A}}\frac{d\sigma}{dt}=\left(\xi+\frac{\eta_{0}}{\eta_{A}}\right)\dot{\gamma}+\frac{\eta_{\infty}}{G_{A}}\frac{d\dot{\gamma}}{dt}.$
(2)

We assume homogeneous flow conditions in our model, and thus Equations (1) and
(2) are given in the scalar form.

[width=3.37in]Fig2.png

Figure 2: A mechanical analog of the canonical three-element viscoelastic
thixotropic model.

The time-evolution of the suspension microstructure is described using a
kinetic equation of the form [21]:

$\frac{d\xi}{dt}=k_{A}(1-\xi)-k_{R}\left|\dot{\gamma}\right|\xi,$ (3)

where $\xi$ is a dimensionless structural parameter that ranges between 0 and
1 representing the unstructured and fully structured states, respectively.
Equation (3) describes two different processes associated with the evolution
of microstructure, namely “aging” and “rejuvenation”, represented by the
kinetic rate constants, $k_{A}$ and $k_{R}$, respectively. Aging is a
spontaneous buildup of microstructure when the fluid is undisturbed; and
rejuvenation refers to a flow-induced breakdown of microstructure when the
fluid is under shear [18]. To minimize the number of fitting parameters, we
use the combination of Equations (2) and (3) to provide a minimal model of
viscoelastic thixotropy. In this model, there is a total of 5 fitting
parameters: infinite-shear viscosity ($\eta_{\infty}$), structural viscosity
($\eta_{A}$), elastic modulus ($G_{A}$), aging rate constant ($k_{A}$), and
rejuvenation rate constant ($k_{R}$). Despite the simplicity of this model, we
still need a statistical inference method to fit it to the experimental flow
curves, as discussed in Sec. II.4.

###  II.4 Markov chain Monte Carlo (MCMC) method

Due to the relatively large number of fitting parameters, we need a
statistical inference method to fit the proposed rheological model to our
experimental data. To this end, we use the Markov chain Monte Carlo (MCMC)
method, which is a class of algorithms that can generate a sequence of random
samples from a probability distribution where direct sampling is difficult. To
setup the problem, we consider a model $\mathcal{M}$ determined by a set of
parameters $\vec{\theta}$ that can make predictions on data
$\vec{\mathcal{D}}$. Here, $\mathcal{M}$ is the five-parameter rheological
model proposed in the previous section with parameters
$\vec{\theta}=[\eta_{\infty},\eta_{A},G_{A},k_{A},k_{R}]$, and the data
$\vec{\mathcal{D}}=\\{\sigma,\dot{\gamma}\\}$l; where the latter contains the
experimentally measured shear stress and shear rate. The parameters
$\vec{\theta}$ have a probability density $\mathbb{P}(\vec{\theta})$
satisfying some prior distributions, or initial guess. The distribution of
interest is the posterior distribution of the density
$\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})$, or the probability that the
parameters $\vec{\theta}$ can provide the best model $\mathcal{M}$ given the
data $\vec{\mathcal{D}}$. We can compute this probability from Bayes’ theorem:

$\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})=\frac{\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})\mathbb{P}(\vec{\theta})}{\mathbb{P}(\vec{\mathcal{D}})},$
(4)

where
$\mathbb{P}(\vec{\mathcal{D}})=\int\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})\mathbb{P}(\vec{\theta})\,d\vec{\theta}$
is a normalizing constant [71]. The parameters that can best fit the
experimental data are the ones maximizing Equation (4), which is given as:

$\vec{\theta}_{\mathrm{best}}=\underset{\vec{\theta}}{\mathrm{argmax}}\left\\{\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})\mathbb{P}(\vec{\theta})\right\\},$
(5)

and also known as the maximum likelihood estimation (MLE).

For the five-parameter model used here, an explicit formula of the posterior
probability $\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})$ is difficult to
obtain, and direct sampling from $\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})$
is computationally prohibitive. Hence, we resort to stochastic sampling
methods (MCMC) to generate observations from the posterior distribution [65].
Specifically, we use the Metropolis-Hastings (M-H) algorithm that was
developed by Metropolis et al. [72] and later generalized by Hastings [73]. A
standard procedure of the M-H algorithm to generate a Markov chain consisting
of the following steps:

  * Step 0. Propose an initial guess $\vec{\theta}=\vec{\theta_{0}}$.

  * Step 1. Denote the current state as $\vec{\theta}$; propose a new state $\vec{\theta^{\prime}}$ according to a transition kernel $q(\vec{\theta}\to\vec{\theta^{\prime}})$.

  * Step 2. Calculate the acceptance ratio transitioning from $\vec{\theta}$ to $\vec{\theta^{\prime}}$:

$\alpha=\min\left({1,\frac{\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta^{\prime}})\mathbb{P}(\vec{\theta^{\prime}})}{\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})\mathbb{P}(\vec{\theta})}}\right).$
(6)

  * Step 3. Draw an arbitrary number $u$ from a uniform distribution $\mathcal{U}[0,1]$. If $u<\alpha$, accept the new state $\vec{\theta}=\vec{\theta^{\prime}}$; otherwise reject $\vec{\theta^{\prime}}$ and keep $\vec{\theta}=\vec{\theta}$. Return to Step 1.

Here, the transition kernel $q(\vec{\theta}\to\vec{\theta^{\prime}})$ is
chosen to be a random walk:

$q(\vec{\theta}\to\vec{\theta^{\prime}}):\vec{\theta^{\prime}}=\vec{\theta}+\mathcal{A}\vec{\epsilon},$
(7)

where $\vec{\epsilon}\sim\mathcal{N}(0,\epsilon_{\mathrm{rw}}^{2}\bm{I})$ is a
normally distributed dimensionless perturbation vector, with
$\epsilon_{\mathrm{rw}}$ being a measure of step size of the random walk and
$\bm{I}$ being the identity matrix; $\mathcal{A}$ is a diagonal matrix
adapting the step size in different dimensions [74, 75], since all five
parameters of $\vec{\theta}$ have different orders of magnitude and units.

The prior distribution of $\vec{\theta}$ is assumed to be normal:
$\vec{\theta}\sim\mathcal{N}(\vec{\theta_{0}},\epsilon_{\theta}^{2}\mathcal{A})$,
where $\epsilon_{\theta}$ is a measure of uncertainty in the prior
distribution, and $\mathcal{A}$ is diagonal so that the components of
$\vec{\theta}$ are independent. The posterior distribution of $\vec{\theta}$
relies on the definition of $\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})$. We
denote the experimental data as $\vec{\mathcal{D}}$, and the data generated by
model $\mathcal{M}$ with parameters $\vec{\theta}$ as
$\vec{\mathcal{D}^{\prime}}$. We can then calculate the “distance” between
$\vec{\mathcal{D}}$ and $\vec{\mathcal{D}^{\prime}}$ by a metric $\rho$,
defined as:

$\rho(\vec{\mathcal{D}},\vec{\mathcal{D}^{\prime}})=\left\|\log(\sigma)-\log(\sigma^{\prime})\right\|,$
(8)

where $\left\|\cdot\right\|$ denotes the $L^{2}$-norm, and $\sigma$ and
$\sigma^{\prime}$ are the experimental and model predicted shear stresses,
respectively. Here, we use the logarithmic stress $\log(\sigma)$ to put an
equal weight to the stress data of low shear rates, which is order of
magnitude smaller than that of high shear rates. Note that the structural
parameter $\xi$ is not included in Equation (8) despite being predicted by the
model $\mathcal{M}$, because $\xi$ is unobservable in the experiments.
Alternative methods quantifying the prediction of unobserved quantities can be
found in [76]. Ideally, the metric
$\rho(\vec{\mathcal{D}},\vec{\mathcal{D}^{\prime}})$ is zero for a “perfect”
model. We assume the metric is normally distributed with a zero mean:
$\rho(\vec{\mathcal{D}},\vec{\mathcal{D}^{\prime}})\sim\mathcal{N}(0,\epsilon_{\rho}^{2})$,
therefore we have the conditional probability:

$\mathbb{P}(\vec{\mathcal{D}}|\vec{\theta})=\left(1/\sqrt{2\pi}\epsilon_{\rho}\right)\exp\left(-\rho^{2}(\vec{\mathcal{D}},\vec{\mathcal{D}^{\prime}})/2\epsilon_{\rho}^{2}\right).$
(9)

[width=3.37in]Fig3.png

Figure 3: A graphical representation of the MCMC method, M-H algorithm.

To visualize the aforementioned MCMC method in a simple manner, a graphical
representation of M-H algorithm is shown in Fig. 3. Starting with the initial
guess $\vec{\theta_{0}}$ (gray dot), the current sample $\vec{\theta}$
performs a random walk to a new location $\vec{\theta^{\prime}}$. Samples with
low $\alpha$ are rejected (red dots marked with X), while samples with high
$\alpha$ are accepted and recorded (blue dots). The Markov chain composed of
all accepted samples will slowly drift towards the MLE that maximizes the
posterior probability $\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})$, or global
minimum in terms of the loss function of the predicted variables (green dot).
The Markov chain will then fluctuate around the MLE and reach a stationary
distribution. One can use the recorded samples from a stationary Markov chain
to approximate the posterior distribution
$\mathbb{P}(\vec{\theta}|\vec{\mathcal{D}})$. Figure 4(a) shows an example of
a stationary Markov chain in a three-dimensional subspace of $\eta_{\infty}$,
$\eta_{A}$, and $G_{A}$, for the model fitting of a kaolinite suspension at
$\phi=0.05$. The recorded samples (blue dots) are used to generate marginal
probability densities (colormap). Figure 4(b) shows the evolution of
parameters $\eta_{\infty}$, $\eta_{A}$, and $G_{A}$ as a function of the
number of iterations for the same kaolinite suspension at $\phi=0.05$,
starting with an initial guess of an order of magnitude estimation given as
$\vec{\theta_{0}}=[\eta_{\infty},\eta_{0},G_{A}]=[10^{-3}~{}\mathrm{Pa.s},10^{-1}~{}\mathrm{Pa.s},10^{-3}~{}\mathrm{Pa}]$,
the Markov chain becomes stationary in approximately 2000 iterations. Samples
taken before the Markov chain is stationary are discarded. Only the samples in
stationary Markov chain are used for the model prediction shown here.

[width = 6.69 in]Fig4.png

Figure 4: (a) Accepted samples (blue dots) from a Markov chain for the fitting
of a kaolinite suspension of $\phi=0.05$, in a three-dimensional subspace of
$\eta_{\infty}$, $\eta_{A}$, and $G_{A}$ of the five-dimensional parametric
space. Projections (colormap) are marginal density functions. (b) Evolution of
the parameters $\eta_{\infty}$, $\eta_{A}$, and $G_{A}$ for a kaolinite
suspension of $\phi=0.05$ as a function of iteration number. The Markov chain
becomes stationary after around 2000 iterations.

##  III Results and Discussions

###  III.1 Model parameter fitting results

We first discuss the fitting results of the five model parameters that provide
estimates of the rheological properties of kaolinite clay suspensions. Table 1
summarizes the values of the five model parameters: $\eta_{\infty}$,
$\eta_{A}$, $G_{A}$, $k_{A}$, $k_{R}$ (see Sec. II.3), for different kaolinite
volume fractions $\phi_{k}$ and sweep time $\Delta t$. Also listed are the
critical shear rate $\dot{\gamma}_{c}$ and yield stress $\sigma_{y}$ values
which are discussed in Sec. III.2. Since MCMC method provides the probability
distributions of model parameters, we report both the mean and the standard
deviation in Table 1. Figures 5(a) and 5(b) show the uncertainty bounds for
the model fittings of the stress ($\sigma$) and the structural parameter
($\xi$), respectively.

Table 1: Fitting results of the five model parameters: $\eta_{\infty}$, $\eta_{A}$, $G_{A}$, $k_{A}$, $k_{R}$; also listed are critical shear rate $\dot{\gamma}_{c}$ and yield stress $\sigma_{y}$. All results are reported with three 3 significant figures. $\phi_{k}$ | $\eta_{\infty}\times 10^{3}$ | $\eta_{A}$ | $G_{A}\times 10^{3}$ | $k_{A}\times 10^{3}$ | $k_{R}\times 10^{3}$ | $\dot{\gamma}_{c}=k_{A}/k_{R}$ | $\sigma_{y}$
---|---|---|---|---|---|---|---
| [Pa.s] | [Pa.s] | [Pa] | [s-1] | [–] | [s-1] | [Pa]
$\Delta t=10~{}$s |  |  |  |  |  |  |
0.05 | $2.33\pm 0.05$ | $0.114\pm 0.002$ | $1.53\pm 0.08$ | $17.1\pm 1.1$ | $4.03\pm 0.28$ | $4.31\pm 0.09$ | $0.481\pm 0.007$
0.12 | $4.26\pm 0.32$ | $0.451\pm 0.006$ | $11.6\pm 0.7$ | $18.7\pm 5.6$ | $2.14\pm 0.19$ | $8.74\pm 0.24$ | $3.90\pm 0.08$
0.19 | $11.4\pm 0.4$ | $1.82\pm 0.02$ | $36.8\pm 0.9$ | $11.2\pm 0.9$ | $3.30\pm 0.25$ | $3.36\pm 0.05$ | $6.07\pm 0.05$
0.26 | $131\pm 11$ | $2.48\pm 0.02$ | $1650\pm 260$ | $11.7\pm 0.8$ | $0.629\pm 0.041$ | $18.6\pm 1.8$ | $43.7\pm 0.3$
$\Delta t=20~{}$s |  |  |  |  |  |  |
0.05 | $2.27\pm 0.06$ | $0.091\pm 0.001$ | $0.84\pm 0.02$ | $17.0\pm 0.8$ | $2.80\pm 0.14$ | $6.07\pm 0.12$ | $0.538\pm 0.008$
0.12 | $3.78\pm 0.29$ | $0.405\pm 0.004$ | $6.72\pm 0.13$ | $8.02\pm 0.63$ | $0.751\pm 0.082$ | $10.7\pm 0.3$ | $4.31\pm 0.08$
0.19 | $9.11\pm 0.07$ | $1.63\pm 0.02$ | $19.6\pm 0.8$ | $7.19\pm 0.31$ | $1.83\pm 0.07$ | $3.93\pm 0.04$ | $6.39\pm 0.02$
0.26 | $122\pm 10$ | $2.26\pm 0.03$ | $1260\pm 360$ | $6.42\pm 0.42$ | $0.310\pm 0.022$ | $20.7\pm 2.5$ | $44.2\pm 0.5$

For all the $\phi_{k}$ values, the fitting parameter, the infinite-shear
viscosity ($\eta_{\infty}$) for kaolinite clay suspensions, is much larger
than the viscosity of water ($\eta_{\mathrm{w}}=10^{-3}$ Pa.s). For example,
the kaolinite suspension of the lowest volume fraction ($\phi=0.05$) has a
infinite-shear viscosity of $\eta_{\infty}\approx 2\eta_{w}$. The structural
viscosity ($\eta_{A}$) of the clay suspensions is approximately 50 to 200
times the $\eta_{\infty}$. Both $\eta_{\infty}$ and $\eta_{A}$ increase
nonlinearly with $\phi_{k}$. As the sweep time increases from $\Delta t=10$ s
to 20 s, we find that the values of both the viscosities $\eta_{\infty}$ and
$\eta_{A}$ decrease slightly by $\sim 5\%$ for all $\phi_{k}$ values. This is
the result of the suspension exhibiting thixotropy [18, 19, 20], where the
viscosity decreases as the suspension is sheared for a longer time.

[width = 3.37 in]Fig5.png

Figure 5: Experimental data (dots) and model fittings (solid curves) of (a)
stress ($\sigma$) and (b) structural parameter ($\xi$) as a function of time
for a kaolinite suspension at $\phi=0.19$. Gray shades are uncertainty bounds
of the fitting results.

Next, we find that the kinetic rate constants $k_{R}$ and $k_{A}$ are (on
average) reduced by half (except for $\phi=$0.05 for $k_{R}$) for all
kaolinite clay suspensions as $\Delta t$ increases from 10 s to 20 s.
Nonetheless, the ratio between these two rate constants,
$\dot{\gamma}_{c}=k_{A}/k_{R}$, increases by approximately 15%. The physical
meaning of the ratio $\dot{\gamma}_{c}$ is discussed in Sec. III.2. The
changes in $k_{R}$ and $k_{A}$ with different $\Delta t$ reveal that: (i)
kaolinite clay suspensions have multiple time scales corresponding to the
relaxation of microstructure; and (ii) the model captures only the dominant
time scales corresponding to aging ($k_{A}^{-1}$) and rejuvenation
($|k_{R}\dot{\gamma}|^{-1}$). As $\Delta t$ increase from 10 s to 20 s, the
dominant time scales increase and the values of $k_{A}$ and $k_{R}$ decrease.
Moreover, we find that the values of elastic moduli ($G_{A}$) are halved, on
average, as $\Delta t$ is doubled. This indicates that kaolinite clay
suspensions have more than a single stress relaxation time scale,
$\lambda=\eta_{A}/G_{A}$; and the value of $\lambda$ for $\Delta t=20~{}$s is
almost twice as long as that for $\Delta t=10~{}$s since $\eta_{A}$ remains
roughly a constant as $\Delta t$ increases (see Table 2). The multi-timescale
nature of kaolinite clay suspensions, as well as the changes in these time
scales with $\phi_{k}$ are discussed in more detail in Sec. III.3.

###  III.2 Yield stress and rheological hysteresis behaviors

[width = 6.69 in]Fig6.png

Figure 6: (a-b) Stress $\sigma$ vs. shear rate $\dot{\gamma}$ flow curves in
linear-linear scale, for kaolinite suspensions of different volume fractions:
(a) $\phi=0.05$ and (b) $\phi=0.19$. Arrows indicate the direction of the
shear rate sweep: first a downsweep (blue, “1”), followed by an upsweep (red,
“2”). Symbols are experimental data; solid curves are fittings of the
rheological model. Yield stress $\sigma_{y}$ is indicated by red dashed lines.
(c-d) The same flow curves as in (a) and (b), but in log-log scale, for (c)
$\phi=0.05$ and (d) $\phi=0.19$. Insets: apparent viscosity,
$\eta_{\mathrm{app}}=\sigma/\dot{\gamma}$, of the same flow curves. (e-f)
structural parameter $\xi$ as a function of $\dot{\gamma}$, obtained from the
model fittings, for (e) $\phi=0.05$ and (f) $\phi=0.19$. Insets: actual
viscosity of the same data, $\eta=\eta_{\infty}+\xi\eta_{A}$.

In this section, we focus on the yield stress and rheological hysteresis
behaviors exhibited in the stress ($\sigma$) vs. shear rate ($\dot{\gamma}$)
flow curves, and how well our model captures these behaviors. Figures 6(a) and
6(b) show the experimental data and model fitting results of
$\sigma$-$\dot{\gamma}$ flow curves in linear-linear scale, for kaolinite
suspensions of $\phi=0.05$ and 0.19, respectively. We find that the stress
increases rapidly with the shear rate until it reaches a critical stress value
$\sigma_{y}$; above $\sigma_{y}$, the stress exhibits a power-law-like
relation with the shear rate. These are typical behaviors of a yield stress
fluid, which are captured by the model fitting. Unlike equations of the
Herschel-Bulkley form [77]: $\sigma=\sigma_{0}+K\dot{\gamma}^{n}$, our model
does not include an explicit yield stress $\sigma_{0}$ in the equation of
stress [ Equation (2)]. Rather, the model captures the yield stress behaviors
through a dramatic decrease in viscosity from $(\eta_{A}+\eta_{\infty})$ to
$\eta_{\infty}$ with increasing shear rate. Nevertheless, one can estimate a
yield stress from the model as:
$\sigma_{y}=(\eta_{A}-\eta_{\infty})\dot{\gamma}_{c}$, where
$\dot{\gamma}_{c}=k_{A}/k_{R}$ is a critical shear rate at which the material
starts to fluidize [29]. We can see that the model predicted yield stress
$\sigma_{y}$ is located relatively accurately at the yield transition of the
flow curves for both $\phi_{k}$ [red dashed lines in Fig. 6(a) and 6(b)].
Above yield, the stress during upsweep (red) overshoots and stays higher than
the stress during downsweep (blue). This rheological hysteresis behavior is
also captured by the model. Table 1 lists the values of the yield stress
$\sigma_{y}$ for different volume fraction $\phi_{k}$ and sweep time $\Delta
t$. We find that $\sigma_{y}$ increases nonlinearly with the volume fraction
$\phi_{k}$, since the additional kaolinite particles increase the viscosities
$\eta_{A}$ and $\eta_{\infty}$. As $\Delta t$ is doubled, $\sigma_{y}$ stays
roughly the same (Table 1), which suggests the yield stress of kaolinite clay
suspensions is approximately time-independent.

Figures 6(c) and 6(d) highlights the rheological behaviors of kaolinite
suspensions below yield stress, by showing the same flow curves as in Fig.
6(a) and 6(b) but in log-log scale. At low shear rate, we find that the stress
during upsweep (red) undershoots and stays below that during downsweep (blue),
until it reaches the yield stress. This creates a hysteresis loop below yield
that is inverted compared to the flow curves above yield, in which case the
stress overshoots during upsweep. Our model captures this double-loop
hysteresis behaviors quite well. Similar double-loop hysteresis behaviors have
been observed for other structured fluids in previous studies [54, 42]. And it
has been proposed that the looping behavior below yield is the signature of
viscoelasticity [42]: the elastic stress initially accumulated at the high
shear rates continues to relax during upsweep, which leads to lower stress
values during upsweep than downsweep. This mechanism requires the time scale
of stress relaxation to be comparable to the time scale of microstructural
relaxation; and we shall see that this is the case for kaolinite suspensions
at these volume fractions in Sec. III.3. Overall, these results suggest that
the behaviors of kaolinite clay suspensions are primarily viscoelastic below
yield transition, which is consistent with previous modelling of structured
fluids [25, 26].

Insets of Figs. 6(c) and 6(d) show the apparent viscosity,
$\eta_{\mathrm{app}}=\sigma/\dot{\gamma}$, of the same flow curves. We find
that $\eta_{\mathrm{app}}$ increases non-monotonically during downsweep
(blue), which is likely due to the coexistence of multiple stress relaxation
times. This non-monotonic behavior is not captured by our model, since our
model only has a single stress relaxation time $\lambda$. During upsweep
(red), the model fits the flow curves relatively better than downsweep. This
can be understood as follows. Since elastic stress is mainly accumulated at
the start of downsweep, faster stress relaxation modes have already fully
relaxed at the end of downsweep. And our model captures the longest relaxation
time $\lambda$ associated with the slowest stress relaxation mode during
upsweep. Another limitation of the model is that it only has a single
microstructural relaxation time $k_{A}^{-1}$, thus it fails to predict non-
monotonic microstructural relaxation behaviors [35]. Both limitations of the
model can be responsible for the discrepancy between the experimental data and
the model fittings.

Figures 6(e) and 6(f) show the structural parameter $\xi$ as a function of
shear rate $\dot{\gamma}$, for kaolinite suspensions of $\phi_{k}=0.05$ and
0.19. Since $\xi$ is unobservable, there is no experimental data shown in the
figures. We find that for all shear rates, the structural parameter during
upsweep (“2”) is always larger than that during downsweep (“1”), leading to
hysteresis loops in the $\xi$-$\dot{\gamma}$ flow curves. Since the fluid
viscosity is directly related to the structural parameter by
$\eta=\eta_{\infty}+\xi\eta_{A}$, the viscosity during the upsweep is also
always higher than that during downsweep, as shown in the insets of Fig. 6(e)
and 6(f). This is expected since downsweep has a much stronger shear history
than upsweep: the high $\dot{\gamma}$ values at the start of downsweep destroy
most microstructure in the fluid, while the low $\dot{\gamma}$ values at the
beginning of upsweep allow more microstructure to build up. These results
provide further evidence that the hysteresis behaviors in Figs. 6(c) and 6(d)
are caused by viscoelasticity. The fluid has a lower stress value despite
having a higher viscosity during upsweep (than downsweep), which indicates the
looping behaviors below yield stress are the results of viscoelasticity,
instead of the change in viscosity due to thixotropy.

###  III.3 Stress vs. microstructural relaxation time scales

Kaolinite suspensions exhibit thixotropy since kaolinite clay particles form
microstructure in the suspensions that physically break down and build up over
time [78, 13, 14, 15]. In this section, we focus on the comparison of two
relaxation time scales: $\lambda$ associated with the relaxation of elastic
stress, and $k_{A}^{-1}$ associated with the relaxation of microstructure in
the suspension. The ratio between the two time scales, $\lambda k_{A}$,
provides a possible classification of thixotropy into “ideal” and
“viscoelastic” [18, 19, 20]. “Ideal” thixotropic materials have a much shorter
relaxation time of elastic stress than that of microstructure, and $\lambda
k_{A}\to 0$; “viscoelastic” thixotropic materials have a stress relaxation
time that is comparable to the microstructural relaxation time, and $\lambda
k_{A}\sim\mathcal{O}(1)$.

Table 2: Stress relaxation time scale $\lambda$, microstructural relaxation time scale $k_{A}^{-1}$, and the ratio between the 2 time scales $\lambda k_{A}$. $\phi_{k}$ |  $\lambda$ [s] |  $k_{A}^{-1}$ [s] |  $\lambda k_{A}$ [-]
---|---|---|---
$\Delta t=10~{}$s |  |  |
0.05 | $77.0\pm 5.4$ | $58.7\pm 3.7$ | $1.32\pm 0.07$
0.12 | $38.8\pm 1.4$ | $54.0\pm 5.6$ | $0.73\pm 0.06$
0.19 | $49.5\pm 3.4$ | $93.1\pm 4.5$ | $0.56\pm 0.02$
0.26 | $1.5\pm 0.4$ | $85.4\pm 6.4$ | $0.02\pm 0.00$
$\Delta t=20~{}$s |  |  |
0.05 | $107.3\pm 1.7$ | $59.1\pm 2.8$ | $1.82\pm 0.08$
0.12 | $60.8\pm 2.7$ | $125.8\pm 6.5$ | $0.49\pm 0.05$
0.19 | $83.1\pm 3.7$ | $139.4\pm 5.8$ | $0.60\pm 0.01$
0.26 | $1.8\pm 0.6$ | $156.2\pm 7.2$ | $0.01\pm 0.00$

Table 2 lists the fitting results of the stress relaxation time $\lambda$, the
microstructural relaxation time $k_{A}^{-1}$, and the ratio between the two
time scales $\lambda k_{A}$, for different volume fractions $\phi_{k}$ and
sweep time $\Delta t$. We first discuss the changes in these time scales with
$\phi_{k}$. We find that both time scales are at the order of 30 seconds to
several minutes (except for $\phi=0.26$, which is discussed later). The stress
relaxation time $\lambda$ decreases with increasing $\phi_{k}$, while the
microstructural relaxation time $k_{A}^{-1}$ increases with $\phi_{k}$.
Nonetheless, the time scale ratio $\lambda k_{A}$ decreases from $\lambda
k_{A}\sim\mathcal{O}(10^{0})$ at $\phi=0.05$ to $\lambda
k_{A}\sim\mathcal{O}(10^{-2})$ at $\phi=0.26$ indicating that kaolinite
suspensions becomes more thixotropic and less viscoelastic. At relatively low
volume fractions ($\phi<0.26$), $\lambda k_{A}\sim\mathcal{O}(1)$. At the
highest volume fraction ($\phi=0.26$), however, $\lambda k_{A}$ decreases
dramatically to $\sim\mathcal{O}(0.01)$, which indicates the kaolinite
suspension is almost ideally thixotropic with negligible viscoelasticity. The
changes in the rheological behaviors of kaolinite suspensions can be explained
as follows. Kaolinite particles are charged platelets; the faces and edges of
kaolinite particles carry opposite charges. As $\phi_{k}$ increases, the
dominant interaction between kaolinite particles transitions from an edge-face
attraction to a face-face repulsion, due to steric effects [14]. This
transition from attractive to repulsive interaction will reduce the
viscoelastic behaviors of kaolinite suspensions [79], which is responsible for
the dramatic decreases in $\lambda$ (from $\sim 30$ s to $\sim 1$ s) and
$\lambda k_{A}$ (from $\sim 0.5$ to $\sim 0.01$) at the highest $\phi_{k}$.

Next, we focus on the changes in stress and microstructural relaxation times
with different sweep time $\Delta t$. We find that both time scales $\lambda$
and $k_{A}^{-1}$ increase as $\Delta t$ is doubled from 10 s to 20 s (by an
average of 53% and 68%, respectively). As discussed in Sec. III.1, this is a
result of the multi-timescale nature of kaolinite clay suspensions, as they
possess multiple stress and microstructural relaxation times. Since the model
used here only has a single stress relaxation time $\lambda$ and a single
microstructural relaxation time $k_{A}^{-1}$, it captures only the dominant
time scales of the relaxation processes. As $\Delta t$ increases, the dominant
time scales captured by the model will increase accordingly, which is shown by
the increase in stress and microstructural relaxation times with $\Delta t$
(Table 2). Since $\lambda$ and $k_{A}^{-1}$ increase by roughly the same
percentage, the ratio between them $\lambda k_{A}^{-1}$ remains similar as
$\Delta t$ is increased.

[width=3.37in]Fig8.png

Figure 7: (a) Examples of kaolinite suspensions possessing linearly correlated
spectra of stress relaxation time scale $\lambda$ and microstructural
relaxation time scale $k_{A}^{-1}$, for $\phi_{k}=0.05$ and $\Delta t=10~{}$s.
(b) Examples of kaolinite suspensions possessing a single dominant stress
relaxation time scale $\lambda$, and multiple dominant microstructural
relaxation time scales $k_{A}^{-1}$, for $\phi_{k}=0.12$ and $\Delta
t=10~{}$s. Here, joint distributions are illustrated by scatter plots of the
accepted samples from the Markov chain (blue dots), and marginal distributions
are shown by the histograms on the side (bar plots); red dashed lines are
linear regressions of the data.

Another piece of evidence of the multi-timescale nature of kaolinite
suspensions can be shown by the probability distribution phase of the
relaxation time scales. Figure 7(a) shows the joint probability distribution
of $\lambda$ and $k_{A}^{-1}$ for $\phi=0.05$ and $\Delta t=10~{}$s. The
marginal distributions (bar plots) of $\lambda$ and $k_{A}^{-1}$ show that
both time scales vary by approximately $\pm 30$ s, and there is no evident
modal time scale (peak) in the distributions, indicating the existence of a
spectrum of time scales. Interestingly, the joint distribution exhibits a
linear correlation between $\lambda$ and $k_{A}^{-1}$, which suggests we can
have a relatively accurate estimate of the intrinsic time scale ratio $\lambda
k_{A}$, despite not having accurate measures of $\lambda$ and $k_{A}^{-1}$
separately. Figure 7(b) shows the joint probability distribution of $\lambda$
and $k_{A}^{-1}$ for $\phi=0.12$ and $\Delta t=10~{}$s. The marginal
distribution of $k_{A}^{-1}$ shows two modal time scales, peaked at 54 s and
62 s, respectively. The marginal distribution of $\lambda$ shows a single
modal time scale at approximately $38.5~{}$s, with narrow tails of a standard
deviation of $\sim 1~{}$s. The joint distribution of $\lambda$ and
$k_{A}^{-1}$ still demonstrates a weak linear correlation, but not as evident
as in Fig. 7(a). Both figures show that kaolinite clay suspensions possess
multiple structural and/or stress relaxation times. These results suggest that
models with multiple or a spectrum of relaxation times, such as, “SGR model”
[36, 37] and “multimode structural kinetics model” [32, 33], may better
capture the rheological behavior of clay suspensions, especially the non-
monotonic relaxation behaviors in Fig. 6(c) and 6(d).

###  III.4 Gravitational aging and microstructural relaxation

We now discuss the origin of “aging”, i.e., the recovery of microstructure in
kaolinite clay suspensions. In many colloidal suspensions, the recovery of
microstructure is considered to be induced by Brownian motion [80, 18, 47, 33,
20, 81]. We shall see that this may not be the case for kaolinite clay
suspensions, since kaolinite particles are essentially non-Brownian. The
extent of non-thermal nature of the constituent particles is estimated using
the gravitational Péclet number, which is the ratio between the time scale of
Brownian motion and the sedimentation time scale [82, 83]:

$Pe_{g}=\frac{4\pi r^{4}\Delta\rho g}{3k_{B}T},$ (10)

where $r$ is the average particle radius, $\Delta\rho$ is the density
difference between kaolinite particles and water, $g$ is the gravitational
acceleration, $k_{B}$ is the Boltzmann constant, and $T$ is the Kelvin
temperature. Kaolinite particles have a modal nominal diameter of $d\approx
10~{}\umu$m [70, 14], or a particle radius of $r\approx 5~{}\umu$m. This leads
to a gravitational Péclet number of $Pe_{g}\sim\mathcal{O}(10^{4})$, which
indicates that kaolinite particles are non-Brownian in nature and aging in
kaolinite suspensions is mainly driven by gravitational settling, rather than
Brownian motion. The aging process in kaolinite suspensions can be understood
as follows. As kaolinite particles sediment, they pile up onto each other, and
aggregate due to electrostatic forces, which constantly generates new
microstructures in the suspensions. We should notice that the system does not
restore its original state by forming new microstructure. Therefore, we prefer
the terminology “restructuring”, instead of “recovery” of microstructure.

[width=3.37in]Fig7.png

Figure 8: Normalized stress $\sigma/\sigma_{y}$ vs. the shear Péclet number
$Pe_{\dot{\gamma}}=9\eta_{\mathrm{w}}\dot{\gamma}/(2\Delta\rho gr)$. We find
that the yield stress $\sigma_{y}$ is attained approximately at
$Pe_{\dot{\gamma}}\sim 1$.

This hypothesis allows us to non-dimensionalize our $\sigma$-$\dot{\gamma}$
flow curves. To re-scale our shear rate, we estimate the shear Péclet number
$Pe_{\dot{\gamma}}$, which is the ratio of the sedimentation time scale and
the time scale of the fluid shear (or experimental time scale):

$Pe_{\dot{\gamma}}=\frac{9\eta_{\mathrm{w}}\dot{\gamma}}{2\Delta\rho gr}.$
(11)

Figure 8 shows the normalized stress, scaled using inherent yield stress
scale, $\sigma/\sigma_{y}$ as a function of the shear Péclet number
$Pe_{\dot{\gamma}}$. We find that in the re-scaled flow curves, the yield
stress $\sigma_{y}$ matches approximately at $Pe_{\dot{\gamma}}=1$. Below
yield stress, $Pe_{\dot{\gamma}}<1$, the fluid shear is too weak to overrun
the gravitational aging. Above yield stress, $Pe_{\dot{\gamma}}>1$, the shear
flow is strong enough to rejuvenate the sample from aging, and the suspension
starts to fluidize. The yield stress of kaolinite suspensions is a result of
the competition between two time scales: the time scale of fluid shear
associated with the breakdown of microstructure, and the time scale of
gravitational settling associated with the buildup of microstructure.
Apparently, the simple scaling of $Pe_{\dot{\gamma}}$ considers neither the
electrostatic force nor the interaction between kaolinite clay particles.
Nonetheless, this limitation does not prevent the scaling from qualitatively
capturing the yield transition of kaolinite clay suspensions.

##  IV Conclusions

To summarize, we experimentally examine the rheology of kaolinite clay
suspensions using a shear rate sweep protocol consisted of consecutive
downsweep and upsweep; the flow curves show both yield stress and rheological
hysteresis behaviors for various kaolinite volume fractions. These behaviors
are captured by a microstructural viscoelastic model using Markov chain Monte
Carlo (MCMC) fitting method. The method also allow us to estimate the
rheological properties of the clay suspensions, such as viscosities, stress
and microstructural relaxation time scales, and the yield stress. The
comparison of stress relaxation time $\lambda$ and microstructural relaxation
time $k_{A}^{-1}$ suggests that kaolinite clay suspensions are strongly
viscoelastic and weakly thixotropic at relatively low volume fractions
($\phi_{k}<0.26$), while being almost inelastic and purely thixotropic at the
highest volume fractions ($\phi_{k}=0.26$). We have also shown that kaolinite
suspensions possess multiple or a spectrum of stress and microstructural
relaxation time scales. Overall, these results can provide insights for
developing predictive models of natural materials like mud and clay, as well
as other yield stress materials and known structured fluids.

###### Acknowledgements.

We thank Yogesh Joshi, Paris Perdikaris, Georgios Kissas, Yibo Yang, and Larry
Galloway for helpful discussions. This work was supported by NSF Grant
DMR-1709763, and partially supported by ARO Grant W911-NF-16-1-0290.

## References

  * Hungr  _et al._ [2014] O. Hungr, S. Leroueil, and L. Picarelli, The varnes classification of landslide types, an update, Landslides 11, 167 (2014).
  * Pimentel [2006] D. Pimentel, Soil erosion: A food and environmental threat, Environ. Dev. Sustain. 8, 119 (2006).
  * Dunne and Jerolmack [2020] K. B. Dunne and D. J. Jerolmack, What sets river width?, Sci. Adv. 6 (2020).
  * Phillips  _et al._ [2022] C. B. Phillips, C. C. Masteller, L. J. Slater, K. B. Dunne, S. Francalanci, S. Lanzoni, D. J. Merritts, E. Lajeunesse, and D. J. Jerolmack, Threshold constraints on the size, shape and stability of alluvial rivers, Nat. Rev. Earth Environ. 3, 406 (2022).
  * Barnes [2010] G. E. Barnes,  _Soil Mechanics: Principles and Practice_ (Palgrave Macmillan, 2010).
  * Van Dyke [1998] K. Van Dyke,  _Drilling Fluids, Mud Pumps, and Conditioning Equipment_ (University of Texas at Austin, Petroleum Extension Service, 1998).
  * Jerolmack and Daniels [2019] D. J. Jerolmack and K. E. Daniels, Viewing earth’s surface as a soft-matter landscape, Nat. Rev. Phys. 1, 716 (2019).
  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89, 035005 (2017).
  * Frigaard [2019] I. Frigaard, Simple yield stress fluids, Curr. Opin. Colloid Interface Sci. 43, 80 (2019).
  * Jogun and Zukoski [1999] S. M. Jogun and C. F. Zukoski, Rheology and microstructure of dense suspensions of plate-shaped colloidal particles, J. Rheol. 43, 847 (1999).
  * Brown  _et al._ [2000] A. B. D. Brown, S. M. Clarke, P. Convert, and A. R. Rennie, Orientational order in concentrated dispersions of plate-like kaolinite particles under shear, J. Rheol. 44, 221 (2000).
  * Zbik  _et al._ [2008] M. S. Zbik, R. S. Smart, and G. E. Morris, Kaolinite flocculation structure, J. Colloid Interface Sci. 328, 73 (2008).
  * Gupta  _et al._ [2011] V. Gupta, M. A. Hampton, J. R. Stokes, A. V. Nguyen, and J. D. Miller, Particle interactions in kaolinite suspensions and corresponding aggregate structures, J. Colloid Interface Sci. 359, 95 (2011).
  * Lin  _et al._ [2015] Y. Lin, N. Phan-Thien, J. B. P. Lee, and B. C. Khoo, Concentration dependence of yield stress and dynamic moduli of kaolinite suspensions, Langmuir 31, 4791 (2015).
  * Neelakantan  _et al._ [2018] R. Neelakantan, F. Vaezi G., and R. S. Sanders, Effect of shear on the yield stress and aggregate structure of flocculant-dosed, concentrated kaolinite suspensions, Miner. Eng. 123, 95 (2018).
  * Mewis [1979] J. Mewis, Thixotropy – a general review, J. Nonnewton. Fluid Mech. 6, 1 (1979).
  * Barnes [1997] H. A. Barnes, Thixotropy – a review, J. Nonnewton. Fluid Mech. 70, 1 (1997).
  * Mewis and Wagner [2009] J. Mewis and N. J. Wagner, Thixotropy, Adv. Colloid Interface Sci. 147-148, 214 (2009).
  * Larson [2015] R. G. Larson, Constitutive equations for thixotropic fluids, J. Rheol. 59, 595 (2015).
  * Larson and Wei [2019] R. G. Larson and Y. Wei, A review of thixotropy and its rheological modeling, J. Rheol. 63, 477 (2019).
  * Moore [1959] F. Moore, The rheology of ceramic slips and bodies, Trans. Br. Ceram. Soc. 58, 470 (1959).
  * Stickel  _et al._ [2006] J. J. Stickel, R. J. Phillips, and R. L. Powell, A constitutive model for microstructure and total stress in particulate suspensions, J. Rheol. 50, 379 (2006).
  * Grillet  _et al._ [2009] A. M. Grillet, R. R. Rao, D. B. Adolf, S. Kawaguchi, and L. A. Mondy, Practical application of thixotropic suspension models, J. Rheol. 53, 169 (2009).
  * Dullaert and Mewis [2006] K. Dullaert and J. Mewis, A structural kinetics model for thixotropy, J. Nonnewton. Fluid Mech. 139, 21 (2006).
  * de Souza Mendes [2009] P. R. de Souza Mendes, Modeling the thixotropic behavior of structured fluids, J. Nonnewton. Fluid Mech. 164, 66 (2009).
  * de Souza Mendes [2011] P. R. de Souza Mendes, Thixotropic elasto-viscoplastic model for structured fluids, Soft Matter 7, 2471 (2011).
  * de Souza Mendes and Thompson [2012] P. R. de Souza Mendes and R. L. Thompson, A critical overview of elasto-viscoplastic thixotropic modeling, J. Nonnewton. Fluid Mech. 187-188, 8 (2012).
  * de Souza Mendes and Thompson [2013] P. R. de Souza Mendes and R. L. Thompson, A unified approach to model elasto-viscoplastic thixotropic yield-stress materials and apparent yield-stress fluids, Rheol. Acta 52, 673 (2013).
  * Blackwell and Ewoldt [2014] B. C. Blackwell and R. H. Ewoldt, A simple thixotropic–viscoelastic constitutive model produces unique signatures in large-amplitude oscillatory shear (LAOS), J. Nonnewton. Fluid Mech. 208-209, 27 (2014).
  * Blackwell and Ewoldt [2016] B. C. Blackwell and R. H. Ewoldt, Non-integer asymptotic scaling of a thixotropic-viscoelastic model in large-amplitude oscillatory shear, J. Nonnewton. Fluid Mech. 227, 80 (2016).
  * Armstrong  _et al._ [2016] M. J. Armstrong, A. N. Beris, S. A. Rogers, and N. J. Wagner, Dynamic shear rheology of a thixotropic suspension: Comparison of an improved structure-based model with large amplitude oscillatory shear experiments, J. Rheol. 60, 433 (2016).
  * Wei  _et al._ [2016] Y. Wei, M. J. Solomon, and R. G. Larson, Quantitative nonlinear thixotropic model with stretched exponential response in transient shear flows, J. Rheol. 60, 1301 (2016).
  * Wei  _et al._ [2018] Y. Wei, M. J. Solomon, and R. G. Larson, A multimode structural kinetics constitutive equation for the transient rheology of thixotropic elasto-viscoplastic fluids, J. Rheol. 62, 321 (2018).
  * Agarwal  _et al._ [2021] M. Agarwal, S. Sharma, V. Shankar, and Y. M. Joshi, Distinguishing thixotropy from viscoelasticity, J. Rheol. 65, 663 (2021).
  * Joshi [2022] Y. M. Joshi, Thixotropy, nonmonotonic stress relaxation, and the second law of thermodynamics, J. Rheol. 66, 111 (2022).
  * Sollich  _et al._ [1997] P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Rheology of soft glassy materials, Phys. Rev. Lett. 78, 2020 (1997).
  * Sollich [1998] P. Sollich, Rheological constitutive equation for a model of soft glassy materials, Phys. Rev. E 58, 738 (1998).
  * Bocquet  _et al._ [2009] L. Bocquet, A. Colin, and A. Ajdari, Kinetic theory of plastic flow in soft glassy materials, Phys. Rev. Lett. 103, 036001 (2009).
  * Mansard  _et al._ [2011] V. Mansard, A. Colin, P. Chauduri, and L. Bocquet, A kinetic elasto-plastic model exhibiting viscosity bifurcation in soft glassy materials, Soft Matter 7, 5524 (2011).
  * Jeyaseelan and Giacomin [2008] R. S. Jeyaseelan and A. J. Giacomin, Network theory for polymer solutions in large amplitude oscillatory shear, J. Nonnewton. Fluid Mech. 148, 24 (2008).
  * Moorcroft  _et al._ [2011] R. L. Moorcroft, M. E. Cates, and S. M. Fielding, Age-dependent transient shear banding in soft glasses, Phys. Rev. Lett. 106, 055502 (2011).
  * Radhakrishnan  _et al._ [2017] R. Radhakrishnan, T. Divoux, S. Manneville, and S. M. Fielding, Understanding rheological hysteresis in soft glassy materials, Soft Matter 13, 1834 (2017).
  * Moorcroft and Fielding [2013] R. L. Moorcroft and S. M. Fielding, Criteria for shear banding in time-dependent flows of complex fluids, Phys. Rev. Lett. 110, 086001 (2013).
  * Moan  _et al._ [2003] M. Moan, T. Aubry, and F. Bossard, Nonlinear behavior of very concentrated suspensions of plate-like kaolin particles in shear flow, J. Rheol. 47, 1493 (2003).
  * Bossard  _et al._ [2007] F. Bossard, M. Moan, and T. Aubry, Linear and nonlinear viscoelastic behavior of very concentrated plate-like kaolin suspensions, J. Rheol. 51, 1253 (2007).
  * Joshi  _et al._ [2008] Y. M. Joshi, G. K. Reddy, A. L. Kulkarni, N. Kumar, and R. P. Chhabra, Rheological behaviour of aqueous suspensions of laponite: new insights into the ageing phenomena, Proc. R. Soc. London A 464, 469 (2008).
  * Joshi and Petekidis [2018] Y. M. Joshi and G. Petekidis, Yield stress fluids and ageing, Rheol. Acta 57, 521 (2018).
  * Rueb and Zukoski [1997] C. J. Rueb and C. F. Zukoski, Viscoelastic properties of colloidal gels, J. Rheol. 41, 197 (1997).
  * Abbasi Moud  _et al._ [2021] A. Abbasi Moud, J. Poisson, Z. M. Hudson, and S. G. Hatzikiriakos, Yield stress and wall slip of kaolinite networks, Phys. Fluids 33, 053105 (2021).
  * Larson [1999] R. Larson,  _The rheology and structure of complex fluids_ (Oxford University Press, New York, 1999).
  * Augusto  _et al._ [2012] P. E. D. Augusto, V. Falguera, M. Cristianini, and A. Ibarz, Rheological behavior of tomato juice: Steady-state shear and time-dependent modeling, Food Bioproc. Technol. 5, 1715 (2012).
  * Glicerina  _et al._ [2015] V. Glicerina, F. Balestra, M. Dalla Rosa, and S. Romani, Effect of manufacturing process on the microstructural and rheological properties of milk chocolate, J. Food Eng. 145, 45 (2015).
  * Glicerina  _et al._ [2016] V. Glicerina, F. Balestra, M. Dalla Rosa, and S. Romani, Microstructural and rheological characteristics of dark, milk and white chocolate: A comparative study, J. Food Eng. 169, 165 (2016).
  * Krystyjan  _et al._ [2016] M. Krystyjan, M. Sikora, G. Adamczyk, A. Dobosz, P. Tomasik, W. Berski, M. Łukasiewicz, and P. Izak, Thixotropic properties of waxy potato starch depending on the degree of the granules pasting, Carbohydr. Polym. 141, 126 (2016).
  * Dimitriou and McKinley [2014] C. J. Dimitriou and G. H. McKinley, A comprehensive constitutive law for waxy crude oil: a thixotropic yield stress fluid, Soft Matter 10, 6619 (2014).
  * Geri  _et al._ [2017] M. Geri, R. Venkatesan, K. Sambath, and G. H. McKinley, Thermokinematic memory and the thixotropic elasto-viscoplasticity of waxy crude oils, J. Rheol. 61, 427 (2017).
  * Divoux  _et al._ [2013] T. Divoux, V. Grenard, and S. Manneville, Rheological hysteresis in soft glassy materials, Phys. Rev. Lett. 110, 018304 (2013).
  * Puisto  _et al._ [2015] A. Puisto, M. Mohtaschemi, M. J. Alava, and X. Illa, Dynamic hysteresis in the rheology of complex fluids, Phys. Rev. E 91, 042314 (2015).
  * Fielding [2020] S. M. Fielding, Elastoviscoplastic rheology and aging in a simplified soft glassy constitutive model, J. Rheol. 64, 723 (2020).
  * Jamali  _et al._ [2019] S. Jamali, R. C. Armstrong, and G. H. McKinley, Multiscale nature of thixotropy and rheological hysteresis in attractive colloidal suspensions under shear, Phys. Rev. Lett. 123, 248003 (2019).
  * Ewoldt  _et al._ [2008] R. H. Ewoldt, A. E. Hosoi, and G. H. McKinley, New measures for characterizing nonlinear viscoelasticity in large amplitude oscillatory shear, J. Rheol. 52, 1427 (2008).
  * Dimitriou  _et al._ [2013] C. J. Dimitriou, R. H. Ewoldt, and G. H. McKinley, Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress), J. Rheol. 57, 27 (2013).
  * Lee  _et al._ [2019] J. C.-W. Lee, K. M. Weigandt, E. G. Kelley, and S. A. Rogers, Structure-property relationships via recovery rheology in viscoelastic materials, Phys. Rev. Lett. 122, 248003 (2019).
  * Donley  _et al._ [2020] G. J. Donley, P. K. Singh, A. Shetty, and S. A. Rogers, Elucidating the G” overshoot in soft materials with a yield transition via a time-resolved experimental strain decomposition, Proc. Natl. Acad. Sci. U.S.A. 117, 21945 (2020).
  * Freund and Ewoldt [2015] J. B. Freund and R. H. Ewoldt, Quantitative rheological model selection: Good fits versus credible models using Bayesian inference, J. Rheol. 59, 667 (2015).
  * Mahmoudabadbozchelou  _et al._ [2021] M. Mahmoudabadbozchelou, M. Caggioni, S. Shahsavari, W. H. Hartt, G. Em Karniadakis, and S. Jamali, Data-driven physics-informed constitutive metamodeling of complex fluids: A multifidelity neural network (MFNN) framework, J. Rheol. 65, 179 (2021).
  * Mahmoudabadbozchelou and Jamali [2021] M. Mahmoudabadbozchelou and S. Jamali, Rheology-informed neural networks (RhINNs) for forward and inverse metamodelling of complex fluids, Sci. Rep. 11, 12015 (2021).
  * Mahmoudabadbozchelou  _et al._ [2022a] M. Mahmoudabadbozchelou, G. E. Karniadakis, and S. Jamali, nn-PINNs: Non-newtonian physics-informed neural networks for complex fluid modeling, Soft Matter 18, 172 (2022a).
  * Mahmoudabadbozchelou  _et al._ [2022b] M. Mahmoudabadbozchelou, K. M. Kamani, S. A. Rogers, and S. Jamali, Digital rheometer twins: Learning the hidden rheology of complex fluids through rheology-informed graph neural networks, Proc. Natl. Acad. Sci. 119, e2202234119 (2022b).
  * Seiphoori  _et al._ [2021] A. Seiphoori, A. Gunn, S. Kosgodagan Acharige, P. E. Arratia, and D. J. Jerolmack, Tuning sedimentation through surface charge and particle shape, Geophys. Res. Lett. 48, e2020GL091251 (2021).
  * Marjoram  _et al._ [2003] P. Marjoram, J. Molitor, V. Plagnol, and S. Tavaré, Markov chain Monte Carlo without likelihoods, Proc. Natl. Acad. Sci. U.S.A. 100, 15324 (2003).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953).
  * Hastings [1970] W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57, 97 (1970).
  * Haario  _et al._ [2001] H. Haario, E. Saksman, and J. Tamminen, An adaptive Metropolis algorithm, Bernoulli 7, 223 (2001).
  * Andrieu and Thoms [2008] C. Andrieu and J. Thoms, A tutorial on adaptive MCMC, Stat. Comput. 18, 343 (2008).
  * Oliver  _et al._ [2015] T. A. Oliver, G. Terejanu, C. S. Simmons, and R. D. Moser, Validating predictions of unobserved quantities, Comput. Methods Appl. Mech. Eng. 283, 1310 (2015).
  * Herschel and Bulkley [1926] W. H. Herschel and R. Bulkley, Konsistenzmessungen von Gummi-Benzollösungen, Kolloid-Z. 39, 291 (1926).
  * Rand and Melton [1977] B. Rand and I. E. Melton, Particle interactions in aqueous kaolinite suspensions: I. effect of ph and electrolyte upon the mode of particle interaction in homoionic sodium kaolinite suspensions, J. Colloid Interface Sci. 60, 308 (1977).
  * Chattopadhyay  _et al._ [2022] S. Chattopadhyay, S. Nagaraja, and S. Majumdar, Effect of adhesive interaction on strain stiffening and dissipation in granular gels undergoing yielding, Commun. Phys. 5, 126 (2022).
  * Dullaert and Mewis [2005] K. Dullaert and J. Mewis, Thixotropy: Build-up and breakdown curves during flow, J. Rheol. 49, 1213 (2005).
  * Wei  _et al._ [2019] Y. Wei, M. J. Solomon, and R. G. Larson, Time-dependent shear rate inhomogeneities and shear bands in a thixotropic yield-stress fluid under transient shear, Soft Matter 15, 7956 (2019).
  * Hunter and Weeks [2012] G. L. Hunter and E. R. Weeks, The physics of the colloidal glass transition, Rep. Prog. Phys. 75, 066501 (2012).
  * Pradeep and Hsiao [2020] S. Pradeep and L. C. Hsiao, Contact criterion for suspensions of smooth and rough colloids, Soft Matter 16, 4980 (2020).