<figure><img src="image_1.png" /></figure>

HAL
open science

# Confined yield stress lubrication flows for cement paste backfill in underground stopes

Paul Vigneaux, Yajian Shao, Ian Frigaard

► To cite this version:

Paul Vigneaux, Yajian Shao, Ian Frigaard. Confined yield stress lubrication flows for cement paste backfill in underground stopes. Cement and Concrete Research, 2023, 164, pp.107038. 10.1016/j.cemconres.2022.107038. hal-03892951

HAL Id: hal-03892951

https://hal.science/hal-03892951

Submitted on 29 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Confined yield stress lubrication flows for cement paste backfill in underground stopes

Paul Vigneaux*

UMPA, CNRS UMR 5669 ENS de Lyon, 46 allée d'Italie, F-69364, Lyon Cedex 07, France & Pacific Institute for the Mathematical Sciences, University of British Columbia, 2207 Main Mall, Vancouver, BC Canada V6T 1Z4 & LAMFA CNRS UMR 7352 Université de Picardie Jules Verne, 33 Rue Saint Leu, F-80000, Amiens Cedex 1, France.

e-mail: paul.vigneaux@math.cnrs.fr

Yajian Shao

School of Civil and Resource Engineering, University of Science and Technology Beijing, No. 30 Xueyuan Road Haidian District, Beijing, China 341000. & Department of Mechanical Engineering, University of British Columbia, 2054-6250 Applied Science Lane, Vancouver, BC Canada V6T 1Z4

Ian A. Frigaard

Departments of Mathematics and Mechanical Engineering, University of British Columbia, 1984 Mathematics Road, Vancouver, BC Canada V6T 1Z2.

Tel: 1-604-822-1316; e-mail: frigaard@math.ubc.ca

## Abstract

Along with the increasing excavation of underground spaces, the mining industry has developed various techniques for refilling these voids. One of them is cemented paste backfill (CBP), coming from the mixing of tailings with water and a binder. These pastes are then transported to underground cavities through pipes and into voids, which they fill in a layering pattern. In this paper we explore the capabilities of a simplified thin-film model to describe the evolution of the paste shape within the stope. The mathematical difficulties are to include the yield stress behaviour of the paste and to deal with the walls in a consistent manner. We show the possible behaviours of the associated solutions to these problems, governed by 2 dimensionless parameters. Two regimes of filling are exhibited, when studying physical parameters typically used in CBP. Finally, we compare the simulations with physical experiments performed in a scaled lab setting using realistic tailings pastes from lead-zinc and nickel mines in North-West China. A good agreement is found between simulations and experiments.

**Keywords:** paste backfill; cement; yield stress fluid; Herschel-Bulkley; lubrication approximation; closed domains; imposed flux; laboratory experiments; numerical simulations

## 1. Introduction

The mining industry is a globally important industry with demand driven largely by increasing consumption and populations. Massive demand means a large amount of mineral consumption, mining waste handling, emissions and other consequent environmental issues (Glotov et al., 2018; Achterberg et al., 1999). As an example, the total volume of solid waste from the mines

of China has been more than 25 billion tons. Approximately 12,000 tailings ponds exist in China, with a total tailings stack of more than 10 billion tons, which increases at a rate of 600 million tons per year (Sun et al., 2018). At the same time, excavation of economic minerals produces a variety of underground spaces named as goaf, stope, cave or gob. According to preliminary statistics, by the end of 2015, the national metal and non-metallic underground mines had a total of 1.28 billion m³ of goaf, distributed in 28 provinces (cities, districts) across the country.

Therefore, during the historical development of min-

<footnote>*Corresponding author</footnote>
ing technology, several methods of refilling or backfilling these voids have been developed. Amongst these a relatively new filling technique is cemented paste backfill (CPB). CPB is a homogeneous material in which tailings are bonded together by a hydraulic binder and mixing water with maximum concentration value that allows transportability. Its components are mixed in a plant, usually located at surface, and then, these mixtures are transported to underground cavities through pipes by pumping or gravity (Fang and Fall, 2018). CPB is potentially one of the best approaches for the management of process tailings since it offers significant environmental, technical and economic benefits (Edraki et al., 2014; Yilmaz and Fall, 2017).

In recent years, there have been some difficulties in the application of paste technology. Free flow and accumulation processes in the stope is one of them. During filling of the stope, the slurry can continue to accumulate from the discharge port and spread to the surrounding area, forming a domed surface. Since the stope is a closed underground space, the discharge port is fixed on the top wall and cannot easily be moved during the filling process. Thus, questions arise as to how to predict the filling process and hence optimise it. Because the slurry requires a certain hardening time (usually 10 to 16 hours), workers can enter the mine to make adjustments, e.g. the stope may require frequent pipe movement.

One common difficulty is to effectively link the paste surface with the top wall. Gaps that are left empty reduce the consumption efficiency of the tailings and compromises the ability of the fill to support any load above. It is thus important to understand the key factors governing the filling process such as: the position of discharge inlet, number of discharge inlets working at the same time, and the optimal/maximum roof-contact rate, etc. At the same time, the position of the discharge inlet is determined empirically with no strict theoretical basis. In this way, multiple pipe transfer operations and random discharge positioning greatly reduce the overall efficiency of the filling operation.

As a consequence we would like to understand better (and predict) the evolution of the paste shape as it enters and fills a stope. In this paper we develop a simple model for the paste flow, based on simplifying the Navier-Stokes equations via a thin-film (or lubrication) scaling approach. We have selected this methodology as it is versatile in being able to incorporate at leading order the shear rheology of pastes and is able to model varying base topologies, as are common in an underground space. It is also very economical computationally, which is important in consider-

ing the relatively long process timescales that would make e.g. three dimensional computational fluid dynamics (3D CFD) simulations very tiresome. Here we start with a simple rectangular (cuboid) geometry.

As a further complication, it is well known that paste slurries have significant non-Newtonian properties, in particular a yield stress. This adds a level of complexity to the models. In the context of paste and thickened tailings we can mention the works of Simms and colleagues, that have built on previously derived Bingham and Herschel-Bulkley (HB) asymptotic models, such as we develop below. Modelling free surface flows of a non-Newtonian material over a varying topography is important to many applications such as avalanches, lava or debris flows. Consequently, these flows have received considerable attention, see e.g. (Ancey, 2007; Balmforth et al., 2014) for a review. (Henriquez and Simms, 2009) compared deposits of gold paste tailings (both in the lab and in the field) with the Bingham fluid lubrication theory of (Liu and Mei, 1989) (essentially in 1D). They found good agreement in terms of shape of the final deposit, both for one layer or multiple layers of successive deposits. In particular, they draw attention to the fact that the equilibrium beach profile is not characterised by a unique angle, and that the overall angle of the deposit is a function of the flow rate. Mizani et al. (2013) extend this study to 3-dimensional deposits and to potential change of the yield stress during the flow. In particular, it is found that shape of the final deposit can sometimes exhibit a change of convexity of the curvature after several layers are deposited. Gao and Fourie (2018), using the bi-viscosity model approximation of Bingham fluids implemented in Ansys Fluent, performed 2D (vertical) and 3D simulations in flume test configurations. They studied the influence of involved parameters on the slope of final profiles, before the slurry touches the opposite wall of the flume.

In this paper, we start with development of a thin-film model for the paste flow within the stope. While much of the derivation of this type of model is standard, there is some novelty and difficulty in dealing with the end wall condition, the correct handling of which drives the layering in the stope. This end wall condition also means that the semi-analytical solutions that are present in the literature are not applicable here, i.e. solutions for unbounded regions. The model is described below in §2 and results are presented in §3 covering the range of dimensionless parameters common in CPB processes. Later in the paper we present validation results from a series of laboratory experiments performed in a dimensionally scaled model using realistic tailings pastes (§4). The agreement is generally good. The paper concludes
with a forward-looking discussion of problem areas that we hope to study using our model in the future.

## 2. Thin film framework for modelling CBP

As discussed in §1, thin film (lubrication) models provide an economical and versatile method for predicting process features in CBP operations. They can be readily adapted to varying filling geometries, underlying topography and pumping procedures. We discuss some of these adaptations in §5, but here focus on the simplest model only. Below we will refer to the filled geometry as a stope, although the terminology used varies.

### 2.1. Rectangular inclined stope

In the context of the modelling of stope CBP filling, we consider at first a 3-dimensional rectangular cuboid of size $\hat{L} \times \hat{H} \times \hat{W}$, inclined by an angle $\varphi$ with respect to the horizontal (cf. Fig 1, (a)). Assuming that $\hat{W}$ is large relative to the typical height variation of a layer (or taking a vertical slice along the middle of the stope), we may simplify the geometry from 3D to 2D, leading to the configuration of Fig 1, (b).

We assume that the paste can be described by a Herschel-Bulkley constitutive law, which allows for both yield stress behaviour and shear-dependency. The rheological properties are the yield stress $\hat{\tau}_y$ (kg.m$^{-1}$.s$^{-2}$ or Pa), the consistency $\hat{\kappa}$ (kg.m$^{-1}$.s$^{n-2}$ or Pa.s$^n$) and the power law index $n$ (dimensionless). It is assumed that these properties are determined independently from paste samples using usual rheometric methods. We denote by $\hat{\rho}$ (kg.m$^{-3}$) the paste density and by $\hat{g} = 9.81$ (m$^2$.s$^{-1}$) the gravitational acceleration. It is assumed that the injected flow rate of paste $\hat{Q}_0$ (m$^3$.s$^{-1}$) is known. All the aforementioned can be considered as input parameters for the model and the aim is to predict the placement of the paste.

We adopt directly the lubrication model of (Balmforth et al., 2006) derived for a Herschel-Bulkley fluid. The premise of such models is that the height variation of the paste layer is small relative to the length $\hat{L}$ and that the underlying balance in the flow is between viscous stresses and gravitational spreading. If $\hat{H}_0$ is a layer thickness scale and $\hat{U}$ a representative lateral velocity, then evidently: $\hat{Q}_0 = \hat{W}\hat{H}_0\hat{U}$.

The balance of gravitational spreading with viscous resistance leads to:

$$ \hat{U} = \frac{\hat{g}\hat{H}_0^3 \cos \varphi}{\hat{v}\hat{L}} = \frac{\hat{Q}_0}{\hat{W}\hat{H}_0}, \quad (1) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Geometry of the stope. (a) 3D view, variables with dimensions. Note that the stope may be inclined with an angle $\varphi$ from the horizontal. (b): a slurry in the stope – 1D model (obtained by taking a vertical slice in the middle of the stope), non dimensionalised variables. The blue curve is the height $h$ of the material at two successive times (the flow is from left to right).</figcaption></figure>
where the characteristic kinematic viscosity is given by

$$ \hat{v} = \frac{\hat{k}}{\hat{\rho}} \left(\frac{\hat{U}}{\hat{H}_0}\right)^{n-1}. \qquad (2) $$

From this we may eliminate $\hat{U}$ and find $\hat{H}_0$:

$$ \hat{H}_0 = \left( \frac{\hat{k}\hat{L}}{\hat{\rho}\hat{g}\cos\varphi} \left( \frac{\hat{Q}_0}{\hat{W}} \right)^n \right)^{\frac{1}{2(1+n)}}. \qquad (3) $$

To obtain the thin film model we assume that the aspect ratio $\epsilon = \hat{H}_0/\hat{L}$ is asymptotically small (i.e. a thin film). The derivation is then quite standard and we refer to (Balmforth et al., 2006) for details. The model is expressed dimensionlessly, with heights scaled with $\hat{H}_0$ and distances along the stope scaled with $\hat{L}$. Note that typically $\hat{H}_0 \ll \hat{H}$: the stope height is not relevant to the dynamics of the layering, but simply determines a maximal height for the backfill.

The dimensionless height of the top of the paste layer is: $z = h(x, t)$ (cf. Fig 1, (b)), at time $t$ and distance $x$. Here time has been scaled with $\hat{L}/\hat{U}$. The evolution of $h$ is given by the following partial differential equation:

$$ \frac{\partial}{\partial t}h(x, t) + \frac{\partial}{\partial x}q(x, t) = 0, \forall x \in [0, 1[, \forall t > 0, \quad (4) $$

where $q(x, t)$ is the flux function (areal flow rate):

$$ q(x, t) = \int_{0}^{h(x,t)} u(x, z, t) dz. \qquad (5) $$

The simplicity of the thin film approach lies in that the lateral velocity component $u(x, z, t)$, can be calculated directly. It is driven by both the inclination of the stope and the gradients of the interface. We find:

$$ u(x, z, t) = \frac{n \left| S - \frac{\partial h}{\partial x} \right|^{1/n}}{n + 1} \operatorname{sgn} \left( S - \frac{\partial h}{\partial x} \right) \times \begin{cases} [Y^{1+1/n} - (Y-z)^{1+1/n}], & 0 \le z \le Y \\ Y^{1+1/n}, & Y < z \le h. \end{cases} \qquad (6) $$

Here $S$ is the slope parameter:

$$ S = \frac{\tan \varphi}{\epsilon}, \qquad (7) $$

and $z = Y(x, t)$ denotes the position of the yield surface. We have seen in (6) that the velocity is plug-like in the range $Y < z \le h$, which is characteristic of yield stress fluids such as the Herschel-Bulkley fluid. For $z \le Y$ the fluid adjacent to the wall is yielded (i.e. a sheared zone

<figure><img src="image_4.png" /><figcaption>Fig. 2. Detailed sketch for the inner variables of the lubrication model (adapted from (Balmforth et al., 2006)), illustrating the yield surface defined by $Y$ and the $z$ profile of the velocity $u$.</figcaption></figure>

where the deformation is non-zero); cf. Fig 2. The yield surface $Y$ is given by

$$ Y(x, t) = \max \left( h(x, t) - \frac{B}{|S - \frac{\partial h}{\partial x}|}, 0 \right), \qquad (8) $$

where $B$ is the Bingham number:

$$ B = \frac{\hat{H}_0 \hat{\tau}_y}{\hat{\rho} \hat{\nu} \hat{U}}. \qquad (9) $$

To understand the physical meaning of $B$, note that we can rewrite as

$$ B = \frac{\hat{\tau}_y}{\hat{\rho}\hat{g}\hat{H}_0\epsilon\cos\varphi}. \qquad (10) $$

This represents the balance between the yield stress and the driving static pressure acting along the slope. We can see that $Y(x, t)$ depends on the constant parameters $(B, S, n)$, as well as on $h$ and $\frac{\partial h}{\partial x}$. The pivotal role of $B$ is seen in (8). Where $B$ exceeds the combined effects of the slope and interface slope $(B, h|S - \frac{\partial h}{\partial x}|)$, we find $Y = 0$. Where $Y = 0$ (see Fig 2) the paste layer does not yield and the flow stops.

We also see that the $(x, t)$ dependency in both $Y$ and $u$, enters only via both $h$ and $\frac{\partial h}{\partial x}$. Consequently, as we integrate $u(x, z, t)$ over $z$ to obtain $q(x, t)$, it follows that $q$ also depends only on $h$ and $\frac{\partial h}{\partial x}$, parametrically also on $(B, S, n)$. Thus formally at least, we can rewrite (4) as:

$$ \frac{\partial h}{\partial t} + \frac{\partial q}{\partial h} \frac{\partial h}{\partial x} = - \frac{\partial q}{\partial (\frac{\partial h}{\partial x})} \frac{\partial^2 h}{\partial x^2}. \qquad (11) $$

We see that in (4), $h(x, t)$ changes in time due to two mechanisms: (i) advection and dispersion, from the 2nd term in (11); (ii) diffusion, from the right-hand side of (11). Both mechanisms are nonlinear. It can be shown that the diffusivity remains positive and this is the mechanism of gravitational spreading. For the purposes of computation, we can integrate the expression
for $u(x, z, t)$ and find $q$ directly:

$$ q(x, t) = \operatorname{sgn}\left(S - \frac{\partial h}{\partial x}\right) \left|S - \frac{\partial h}{\partial x}\right|^{1/n} \times \frac{nY^{1+1/n}}{(n+1)(2n+1)} ((2n+1)h - nY), \quad (12) $$

where sgn(.) is the usual sign function.

To close the model mathematically, we need to have both initial and boundary conditions for (4). The flows studied here are different from those of Balmforth et al. (2006), although we note that the same authors have treated elsewhere similar problems, such as the flow of lava through a vent. For the stope, we assume that the paste is pumped and deposited at the start, $x = 0$, whereas a wall at the end of the stope at $x = 1$ prevents flow (cf. Fig 1, (b)). We assume therefore, the following initial condition:

$$ \forall x \in [0, 1], \quad h(x, t = 0) = 0 \quad (\text{initial condition}) \quad (13) $$

and the boundary conditions

$$ \forall t > 0, \quad q(x = 0, t) = 1 \quad \text{and} \quad (14) $$

$$ \forall t > 0, \quad q(x = 1, t) = 0. \quad (15) $$

## 2.2. Physical parameter ranges

There are many goafs, stopes and similar underground spaces that might be filled using cemented paste backfill (CBP). Broad parameter ranges may be estimated from (Yilmaz and Fall, 2017; Henriquez and Simms, 2009; Mizani et al., 2013) and other similar studies. These give typical ranges: $B = 0.5 - 250$, $S = 1 - 120$. Concerning the power law index, although both shear-thinning and shear-thickening pastes may be found rheologically, generally $n \approx 1$.

We note that the range taken for $(B, S)$ also includes extreme values and more commonly we might find CBP parameters around $(B, S) \sim (150, 12)$. Also note that we have assumed that the angle $\varphi = 0.5 - 5$ (degrees), in defining $S$. In particular, this means that the stope is always inclined in the direction of flow. A perfectly horizontal stope, $S = 0$, is possible but not considered and similarly there might be situations in which $S < 0$. For the latter, one could consider an equivalent problem with the inflow specified at $x = 1$.

Lastly, note that for typical values the Reynolds number associated with the thin-film approach, $Re = (\hat{U}\hat{H}_0^2)/(\hat{L}\hat{v})$ is generally very small, which is required for the continuing validity of the method.

## 2.3. Computational overview

The evolution equation (4), together with initial condition (13) and boundary conditions (14)-(15), is solved numerically using an upwind (space) discretisation for the flux term (LeVeque, 2002), implemented within Matlab. The time discretisation is done with an explicit Euler scheme on the discrete unknown $h$. For given spatial discretisation ($\Delta x$) we adapt the time step ($\Delta t$) dynamically at each time iteration. This is done by taking the minimum of the two classical numerical stability constraints on $\Delta t$ associated to (4) (see also 11). Namely, a constraint due to the advection component of (4) (CFL condition) and a constraint due to the explicit treatment of the diffusion component ($\Delta t$ constrained by $\Delta x^2$) (LeVeque, 2002). In the simulations presented in this article, we have used 150 points to discretise $x \in [0, 1]$. We checked that the results are converged under a mesh refinement of the discretisation.

Due to the inflow boundary condition ($q = 1$) at $x = 0$ and the upwind nature of the scheme, the time evolution of $h(x = 0, t)$ is obtained by solving the associated discrete non-linear problem $q(h(x = 0)) = 1$ at $x = 0$. This is done with the Matlab `fzero` function, at each time step.

# 3. Computational results and parameter variations

In this section we present the various typical behaviours found in the $(B, S)$ parameter space, fixing $n = 1.0$ to begin with. Figure 3 shows the evolution of $h$ for 4 extremes of the $(B, S)$ parameter space considered.

We observe that the model shows 2 types of qualitative behaviour. In the first case, we have a relatively rapid steady propagation of a near-constant height wavefront towards the end wall $x = 1$, followed by increase of the layer height as a backwards propagating wave. This tends to happen for larger $S$ (e.g. $S = 120$ here). Secondly, where the $h$ profile is always increasing in time at every $x$, (seen e.g. in the case $B = 250$, $S = 1$). The case $B = 0.5$, $S = 1$ is somehow at the transition between these two observed behaviours: showing slightly the backwards wave behaviour.

In these plots, the two thicker black curves depict the interface when it first reaches the wall at $x = 1$ and at a later final time after the front climbs this wall. The backwards propagating interfacial wave does not imply a negative fluid velocity, as we shall see later. To understand this transition we need to dig into the boundary condition at $x = 1$. Setting the volume flux $q = 0$ is of course natural at the wall, but has multiple solutions for $h$ and $\frac{\partial h}{\partial x}$. During the initial propagation of the wave to
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Evolution of *h* profiles (propagation is from left to right as time increases). The profile at the contact time with the wall is marked with a broken black line. The final profile (thick black curve) is computed at the time given in the header of each sub-figure. Here, *n* = 1.</figcaption></figure>
the wall, we have $h(1, t) = 0$, which means $Y(1, t) = 0$ and $q = 0$ is satisfied.

On attaining the end wall, since $q > 0$ within the propagating layer and $q = 0$ at $x = 1$, we find $\frac{dq}{dx} < 0$ at $x = 1$ and from (4) we see that $h$ must increase. This increase in $h$ continues in time, so how is $q = 0$ realised when $h > 0$? The variable $Y$ gives the distance above the bottom of the stope where the paste becomes unyielded. Initially, $\frac{\partial h}{\partial x} < 0$ along the propagating layer and hence the slope of the paste layer reinforces $S$. If $|S - \frac{\partial h}{\partial x}| < B/h$, the combined effects of stope and interface slope are insufficient to yield the fluid at the wall. For fixed yield stress (effectively $B$) there is a range of $h$ and $\frac{\partial h}{\partial x}$ for which this is satisfied. Thus, we may continue with increasing $h$ while the fluid at the end wall remains unyielded.

The switch to a backwards propagating wave can occur when $q \approx 1$ everywhere away from the end wall. This happens naturally when we have these near-constant propagating layers, i.e. with $\frac{\partial h}{\partial x}$ very small except at the sharper front. In this case, the decrease in $q$ becomes most rapid near the end wall and hence the increase in $h$ is also largest there. We switch from negative to positive $\frac{\partial h}{\partial x}$ and the wave propagate backwards. For the backwards wave the slope of the interface now opposes that of the base of the stope. Note, large $S$ arises either from a larger angle of inclination or a smaller aspect ratio. Therefore, although in the dimensionless variables we see a significant positive gradient in $h$, we do not see a change of sign in $q$, which would require that $\frac{\partial h}{\partial x} > S + B/h$.

This latter point becomes apparent if we look at the velocity profiles during the interface propagation. Figure 4 shows the velocity profiles at two $x$ positions, shortly after the front has crossed $x \sim 0.8$, for the same parameter pairs as in Fig. 3. The layer thickness decreases monotonically towards the front and the velocity too. For the two cases with $B = 0.5$ there is barely any yield surface visible, just close to the interface. For the two cases with $B = 250$ there is only a thin yielded layer near the lower wall and most of the paste moves as a plug.

Considering now a time just after contact with the wall at $x = 1$, Fig. 5 shows that all the velocities are still positive, vanishing at $x = 1$. For the two cases with $S = 120$ the velocity under the sharply inclined part of the interface is greatly diminished, i.e. due to the competition between the slope of the interface and the base of the stope.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Vertical velocity profiles at $x = 0.07$ and $x = 0.77$ for a surface $h(x, t)$ just as the front of the paste layer crosses $x \sim 0.8$. Parameters (B, S, t) are in the sub-figure headers and the maximum velocity $umax$ for each profile is also given.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 5. Vertical velocity profiles at x = 0.07 and x = 0.91 for a surface h(x, t) at some time t after contact with the wall. Parameters (B, S, t) are in the sub-figure headers and the maximum velocity umax for each profile is also given.</figcaption></figure>

### 3.1. Front shape and propagation speed

As we have seen in Fig. 3, the initial propagation of the front to the end wall is strongly influenced by (B, S). We compute the front position $x_f(t)$ during this initial propagation, as shown in Fig. 6. As might be expected the front speed decreases with B and increases with S. It is of interest to note that the contact time, (i.e. that at which the front reaches the end wall), $x = 1$ is of O(1) over the range of practical (B, S) studied, confirming that the scaling chosen was appropriate for the initial propagation.

The contact time variations are plotted in Fig. 7 for clarity. In order to recover dimensional times from these results, recall that the timescale is $\hat{L}/\hat{U}$, which may be expressed in terms of the rheology and process variables as:

$$ 
\begin{aligned}
\text{Timescale} &= \frac{\hat{L}}{\hat{U}} = \frac{\hat{W}\hat{L}\hat{H}_0}{\hat{Q}_0}, \\
&= \left[ \frac{\hat{k}}{\hat{\rho}\hat{g}\hat{L}\cos\varphi} \left( \frac{\hat{W}\hat{L}^2}{\hat{Q}_0} \right)^{2+n} \right]^{\frac{1}{2(1+n)}} = \left[ t_v^n t_f^{n+2} \right]^{\frac{1}{2(1+n)}}.
\end{aligned}
\quad (16) $$

In (16) we see that the timescale results from a combination of a viscous timescale ($\hat{t}_v = [\hat{k}/(\hat{\rho}\hat{g}\hat{L}\cos\varphi)]^{1/n}$), representing a balance of density driven slumping with viscous resistance, and a filling timescale $\hat{t}_f = \hat{W}\hat{L}^2/\hat{Q}_0$. Note that $\hat{Q}_0/(\hat{W}\hat{L})$ is the mean vertical filling velocity.

Although we attach no deeper significance to the timescales above, it should not be surprising to see the density vs rheology balance appear. The other effects of rheology of course come in the yield stress, captured in B, which we have already seen has a profound effect on the shape of the propagating paste layer. To illustrate this in more detail Fig. 8 shows the variation in the layer shape at the initial contact with the end wall.

### 3.2. Variations with n

Recall that all the previous results were computed with $n = 1$. We also performed the same sets of (B, S) simulations by taking first $n = 0.8$ and $n = 1.2$. The obtained equivalent figures of Figures 3 to 8 are virtually the same in terms of both qualitative and quantitative behaviours for most of (B, S) values. Larger $n$ is not common and for shear thinning effects to be visible, significantly smaller $n$ is required. However, this is not encountered in the CPB context.

The case for which the difference is the biggest is for (B, S) = (0.5, 120), and for brevity we only show this case to illustrate how close the results are. Figures 9 and 10 can be compared with Figs. 3 and 4, respectively. We
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 6. Evolution of $x_f(t)$ for $n = 1$: a) $B = 0.5$, $S$ varying; b) $B = 250$, $S$ varying; c) $S = 1$, $B$ varying; d) $S = 120$, $B$ varying.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 7. Contact time as a function of (B, S) for n = 1.</figcaption></figure>

observe that increasing *n* (while keeping (B, S) fixed) induces an increase of the maximal height *h* and a slight decrease of velocity, hence an increase of the typical times. The qualitative behaviour is however the same for these 3 values of *n*. This is an interesting observation from the perspective of comparisons with the experiments of the next section.

**4. Experimental verification**

We have also carried out a series of experiments to verify that the simple model developed above is able to capture the key features of filling dynamics in large aspect ratio slightly inclined stopes, using realistic pastes. These pastes were made from field materials, intended for filling of two mines: including all materials and their composition. Both mines have adopted the full tailings paste filling process, which consists of only thickened tailings slurry and cement. The paste description and cement dosage are shown in Table 2.

**4.1. Materials**

The paste slurry must satisfy characteristics of non-segregation, non-bleeding and non-layering. For a given composition material, the variable range of slurry concentration is thus limited. The driving force in the paste flow process in the stope mainly comes from gravity; therefore, the bulk density of the paste slurry plays an important role. In order to obtain paste materials with different rheological parameters and bulk density, two tailings were used as experimental materials in this study.

Tailings samples were obtained from a lead-zinc mine (mine-A) and a nickel mine (Mine-B) located in Northwest China, named Tailings-A and Tailings-B. The grading curve of the tailings was detected using the Top-Sizer Laser Particle Size Analyzer Instrument (OMEC

Company, UK), as shown in Fig. 11. The fines (<20 µm) content of the tailings were determined to be 26.13% and 54.71%, and the tailings were classified as fine and coarse tailings respectively (ASTM D2487-17). The binder is Portland cement (P.C. 42.5) formulated according to the Chinese standard (GB175-2007)), and the mixing water is urban tap water in this study.

Other geotechnical properties of tailings, such as specific gravity, porosity, liquid/plastic limit, etc were obtained according to the test standards (ASTM D4318-17e1). The chemical composition of the tailings samples (Table 1) was determined by X-Ray Fluorescence Spectrometer (EDX 8000, SHIMADZU Company, JPN) and wet chemical methods using ICP-OES Spectrometer (ICPOES730, Agilent Company, USA) after acid digestion according to literature (Ercikdi et al., 2013). The test water was deionised and deoxidised before use.

The tailings were tested in an Ultima IV diffractometer with monochromatic radiation of Cu Kα (40 kV; 40 mA). A scanning speed of 0.02 every 0.5s and scanning angles between 10 and 80 in units of 2-theta, were used. The primary minerals in the tailings sample are quartz, a small amount of calcite and mayenite, as shown in Fig. 12.

**4.2. Rheology test**

The tailings and cement were thoroughly mixed under dry conditions. Subsequently, the dry mixture was combined with water in a planetary mortar mixer (Model No. JJ-5), as described in the ASTM C305 standard (ASTM C305-14), and according to our specified concentrations. Two types of test were performed: a slump test and rotational rheometry.

For the slump test, a truncated cone mould was used downscaled from the ASTM-slump-cone geometry (ASTM C143/C143M-20): 36 mm top diameter, 60 mm bottom diameter, 60 mm height, as described in the methods for test uniformity of concrete admixture (GB/T 8077-2012). The paste was dropped onto a calibrated plane marked with grid squares of 1 × 1 mm². Freshly mixed paste was immediately cast into the cone, and then the cone was lifted as slowly as possible (<0.005 m/s) to minimise inertial effects (Tan et al., 2017). The columnar structure slowly transforms into a diffusive spreading/slumping sample. After 60 s, the final diameter (workability) was recorded; see Table 2.

For the rheometry, the paste specimen were held in cylindrical glass beakers. The height of the thickened tailings was about 100 mm. After casting the paste slurry into the beakers, the specimens were tamped 25 times using a glass rod (5 mm diameter) (Wu et al.,
<figure><img src="image_16.png" /><figcaption>**Fig. 8.** Variation of *h* profiles at the time of wall contact with the (*B*,*S*) given in the legend and *n* = 1. Values are grouped to ease the visualisation.
(a) *B* = 0.5, (b) *B* = 250, (c) *S* = 1, (d) *S* = 120.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 9.** Evolution of *h* profiles (propagation is from left to right as time increases) for (*B*, *S*) = (0.5, 120), to be compared with Figure 3: (a) *n* = 0.8; (b) *n* = 1.2. The profile at the contact time with the wall is marked with a broken black line. The final profile (thick black curve) is computed at the time given in the header of each sub-figure.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 10.** Vertical velocity profiles at *x* = 0.07 and *x* = 0.77 for a surface *h*(*x*, *t*) just as the front of the paste layer crosses *x* ~ 0.8. Parameters (*B*, *S*, *t*) are in the sub-figure headers and the maximum velocity *umax* for each profile is also given: (a) *n* = 0.8; (b) *n* = 1.2, to be compared with Fig. 4.</figcaption></figure>

<table><caption>Table 1: Chemical and physical properties of tailings and cement.</caption><thead><tr><th>Chemical composition(%)</th><th>Tailings-A</th><th>Tailings-B</th><th>Cement</th><th>Physical properties</th><th>Tailings-A</th><th>Tailings-B</th><th>Cement</th></tr></thead><tbody><tr><td>Si</td><td>17.906</td><td>19.58</td><td>13.454</td><td>Specific gravity</td><td>2.74</td><td>2.966</td><td>3.092</td></tr><tr><td>Fe</td><td>0.7</td><td>0.39</td><td>0.89</td><td>MDD(kg m<sup>-3</sup>)</td><td>1.551</td><td>1.617</td><td>1.384</td></tr><tr><td>Ca</td><td>25.107</td><td>2.596</td><td>35.971</td><td>Porosity(%)</td><td>43.39</td><td>45.48</td><td>55.24</td></tr><tr><td>Mg</td><td>2.052</td><td>18.979</td><td>0.08</td><td>Liquid limit(%)</td><td>15.31</td><td>19.45</td><td>--</td></tr><tr><td>Al</td><td>3.626</td><td>2.021</td><td>4.749</td><td>Plastic limit (%)</td><td>13.52</td><td>15.65</td><td>--</td></tr><tr><td>K</td><td>0.04</td><td>0.363</td><td>0.889</td><td>OWC(%)</td><td>13.69</td><td>16.69</td><td>--</td></tr><tr><td>S</td><td>0.076</td><td>1.6012</td><td>1.376</td><td>SSA(cm<sup>2</sup>g<sup>-1</sup>)</td><td>3257.7</td><td>3581.5</td><td>4175.4</td></tr><tr><td>Na</td><td>0.059</td><td>0.009</td><td>0.007</td><td>D<sub>10</sub></td><td>7.274</td><td>2.185</td><td>3.014</td></tr><tr><td>Mn</td><td>1.836</td><td>0.129</td><td>0.124</td><td>D<sub>30</sub></td><td>24.574</td><td>7.962</td><td>14.582</td></tr><tr><td>Ti</td><td>0.186</td><td>0.204</td><td>0.366</td><td>D<sub>50</sub></td><td>52.751</td><td>16.972</td><td>26.104</td></tr><tr><td>Cr</td><td>0.017</td><td>0.312</td><td>--</td><td>D<sub>60</sub></td><td>76.152</td><td>24.031</td><td>31.902</td></tr><tr><td>Ni</td><td>0.011</td><td>0.001</td><td>--</td><td>D<sub>90</sub></td><td>262.121</td><td>74.963</td><td>82.452</td></tr><tr><td>Zn</td><td>0.026</td><td>0.005</td><td>0.012</td><td>C<sub>u</sub></td><td>10.47</td><td>11.00</td><td>10.58</td></tr><tr><td>Cu</td><td>--</td><td>0.241</td><td>0.015</td><td>C<sub>c</sub></td><td>1.09</td><td>1.21</td><td>2.211</td></tr><tr><td>C</td><td>--</td><td>0.022</td><td>--</td><td>Fine(&lt;20 µm)(%)</td><td>26.13</td><td>54.71</td><td>43.56</td></tr><tr><td>LOI</td><td>15.351</td><td>10.254</td><td>14.823</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Notes: MDD = Max dry density; OWC = Optimum water content; SSA = Specific surface area; LOI = Loss on ignition; C<sub>u</sub> = coefficient of uniformity (D<sub>60</sub>/D<sub>10</sub>); C<sub>c</sub> = coefficient of curvature ((D<sub>30</sub>)<sup>2</sup>/(D<sub>10</sub>×D<sub>60</sub>)); -- indicates that the content is lower than the detection threshold. The chemical composition data of the test materials did not include oxygen content. If converted to oxide form, the elemental content totals approximately 100%.
<figure><img src="image_20.png" /><figcaption>Fig. 11. Particle size distributions of tailings A, B and cement.</figcaption></figure>

2015). A R/S four-paddle rotational rheometer (Brookfield RST, MA, USA) was used to test the rheological behaviour of the slurry samples. During the experiment, a four-paddle rotor/vane (V40–20–3tol) was immersed into the slurry, rotating at controlled shear rate (CSR). The temperature of the water bath equipped on this equipment remained at 20 ± 2 °C, which is approximately that of the lab in which experiments were later performed.

During pre-testing, it was found that there is a significant stress overshoot during the initial part of the CSR ramp (Panchal et al., 2018). To eliminate this, following Han and Ferron (2016), a low creep shear preconditioning process was applied prior to conducting the flow curve test, (cf. Fig. 13). For the CSR ramp we raise from 1 to 160 s⁻¹ in increments of 1 s⁻¹ via a step-up approach. Under shear, the constitutive relationship between shear rate and shear stress appears to be approximately represented by the Bingham model or Herschel-Bulkley (HB) model. The rheology parameters are obtained by a least squares regression method based on the test results (see Table 2). In order to avoid the stress overshoot, we fit the rheological models to data in the range 10 to 120 s⁻¹. Evidently the degree of shear-thinning is limited at higher shear rates, i.e. n ≈ 1. Due to the observed stress overshoot, a second estimate of the yield stress of the paste was made at the initial moment of deformation using a “creep shear” test (or stress growth test) (Roussel and Coussot, 2005). During the creep shear test, the optimum speed 2.0 rpm (0.5 s⁻¹) for 30 s was adopted as shear speed, based on the research report ICAR-105-3F (Koehler and Fowler, 2004). The peak stress (τ̂<sub>peak</sub>) was recorded from the stress creep test.

### 4.3. Box backfill experiment

We conduct our experiments in a laboratory setup that is approximately 1/40 of the size of an onsite stope. The lab model has height of 185 mm and total length and width of 778 mm and 130 mm. Recall that the aforementioned thin film model is based not only on geometric similarity but also on the scaling assumptions to reduce the Navier-Stokes equations to the thin-film approximation. This requires that the typical layer thickness variation is significantly smaller than the length of the stope. In our experiments this is verified *a posteri* ori, with typical aspect ratios < 1/20. The model is constructed from acrylic board, allowing side visualisation. The sidewalls were lubricated with hydrophobic grease, to reduce velocity gradients across the width. The bottom is covered with emery cloth to simulate the rough base of the stope floor. The laboratory model stope is shown in Fig. 14.

During the backfill process, the paste slurry is transported by long pipelines (~ km), so that it enters the stope well sheared and with approximately steady flow rate. In the laboratory model the paste slurry is loaded in an upper reservoir maintained at approximately constant height throughout the experiment. In the laboratory model we can control the flow rate by either changing the height in the upper reservoir or by varying the discharge port diameter. Our experiments used flow rates in the range 0.4 - 4.0 litres per minute (LPM). A camera is used to record evolution of the slurry interface. A typical flow at different times during an experiment is shown in Fig. 15.

### 4.4. Results

Figure 16 shows the height variation of the paste layer for Pastes 1-4, each at different representative times during experiments. The experiments are compared with numerical simulation results for two values of yield stress: the one fitted to the Bingham model (B) and the one calculated from the creep shear experiment (*B*<sub>creep</sub>). For the numerical simulations, here are the associated values in terms of (B, S): Paste 1 (S = 0.22, B = 43, *B*<sub>creep</sub> = 58), Paste 2 (S = 0.28, B = 102, *B*<sub>creep</sub> = 122), Paste 3 (S = 0.28, B = 74, *B*<sub>creep</sub> = 98), Paste 4 (S = 0.26, B = 68, *B*<sub>creep</sub> = 93). Note that with these values of (B, S), simulations exhibit a *h* which is always increasing with the time *t*, for all Pastes (see the discussion on the two types of regimes mentioned in section 3).

We observe that the paste surface profile in the physical experiments is consistent with the numerical simulation results. Different rheological parameters of the
<figure><img src="image_21.png" /><figcaption>Fig. 12. Mineral composition of tailings (a) Tailings-A; (b) Tailings-B.</figcaption></figure>

<table><caption>Table 2: Rheology parameters of paste slurry.</caption><thead><tr><th rowspan="2">No.</th><th>Concentration</th><th rowspan="2">Cement dosage/%</th><th>Creep shear</th><th colspan="2">Bingham model</th><th colspan="2">Herschel-Bulkley</th><th colspan="2">Workability</th></tr><tr><th>%</th><th>τ<sub>peak</sub>/Pa</th><th>τ<sub>0</sub>/Pa</th><th>μ<sub>p</sub>/Pa.s</th><th>τ<sub>0</sub>/Pa</th><th>κ̂/Pa.s<sup>n</sup></th><th>n</th><th>/mm</th></tr></thead><tbody><tr><td>Paste-1</td><td>62.6</td><td>10</td><td>63.00</td><td>47.10</td><td>0.167</td><td>48.51</td><td>0.169</td><td>0.995</td><td>97 96</td></tr><tr><td>Paste-2</td><td>63.6</td><td>10</td><td>80.00</td><td>66.84</td><td>0.132</td><td>64.98</td><td>0.086</td><td>1.129</td><td>87 85</td></tr><tr><td>Paste-3</td><td>63.1</td><td>10</td><td>63.00</td><td>47.52</td><td>0.196</td><td>50.52</td><td>0.078</td><td>1.137</td><td>92 96</td></tr><tr><td>Paste-4</td><td>76.8</td><td>25</td><td>85.30</td><td>62.76</td><td>0.318</td><td>64.52</td><td>0.239</td><td>1.041</td><td>112 114</td></tr></tbody></table>

Notes: Concentration: the solid mass concentration of the paste slurry; Cement dosage: the percentage of cement in the solid material of the paste slurry; τ₀: Yield stress, μₚ: Plastic viscosity; n: Power law index; κ: Consistency; τₚₖₐₗ: Yield stress at creep shear or stress growth test stage; Workability: The diffusion radius of the paste slurry originates from the miniature slump cylinder (Cone diameter: φₐₗₐₚ=36 mm, φₐₗₐ₉₆=60 mm; height=60mm). Paste-1, 2 and 3 are from Tailings-A and Paste-4 is from Tailings-B.

<figure><img src="image_22.png" /><figcaption>Fig. 13. Rheology test program.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 14. Stope backfill simulation system (Length 77.8 cm, width 13 cm, height 18.5 cm) (a) outlet with diameter 8 ~ 12 mm; (b) Bottom covered with emery cloth.</figcaption></figure>

slurry do have observable effects on the contour of the paste surface, but for the range of B and S in the experiments, effects are minor. The layer thickness obtained by the simulation based on the Bingham model data fit is smaller than that computed with the yield stress from the peak stress value in the creep test. The creep test value Bₐₗₐₚ produces simulation results closer to the experimental values.

We can see that during the initial stage (before the slurry touches the right wall), the agreement is very good in all cases. However, as the slurry touches the
<figure><img src="image_24.png" /><figcaption>Fig. 15. Paste flow example for Paste 1: (a) initial stage; (b) touching the wall; (c) rising of the level; (d) final profile.</figcaption></figure>

right wall and gradually rises, the difference between simulation and experiments becomes more apparent, especially away from the right wall, close to the injection point. Part of this is due to the mound formed directly under the discharge port. This part of the flow is not 2D and it is not surprising that our model cannot capture this feature. Fig 17 shows an example of the mound flow for Paste 2.

**5. Summary and discussion**

In this paper, we have adapted the viscoplastic lubrication model of (Balmforth et al., 2006) to the geometric configuration of a closed underground stope, in order to study the cemented paste backfill (CBP) process. The study was conducted by using values of the dimensionless parameters ($B, S$) stemming from typical CBP applications. Starting from an initial empty stope and filling it from one side with a yield stress fluid, we showed through numerical simulations that the model exhibits two types of evolution of the free surface. First, the case of a near-constant wave front that rapidly reaches the end wall, followed by a backwards propagating wave of increasing layer height. The second type consists of a height profile which is increasing for all space and at all times. The first type of filling tends to happen for large slope parameter $S$, while the second type occurs for smaller $S$. We also showed that the results vary only marginally with the power law index of the paste $n$, over the usually observed range of $n$ from 0.8 to 1.2.

We also performed laboratory experiments using realistic pastes coming from mines in NW China. We compared these laboratory experiments with the numerical simulations using aforementioned model. It appears that the comparison of results gives very good agreement, even after the paste touches the wall. One source of discrepancy was the fact that 3D mounding can happen close to the injection point, with such thick yield stress pastes, i.e. the surface only becomes 2D at a distance from the injection pipe. This 3D effect is obviously not taken into account by the 2D model. When the flow height becomes high (i.e. close to the stope's roof), 3D mounding can be pronounced and the surface can touch the injection point, leading to injection obstruction and potentially stopping of the flow, before the cavity is totally filled. However, for most of the filling process, the 2D model is a good tool for the simulation and prediction of these flows.

A key point here is to note that the model validity rests with the *slope* of the interface remaining small, since this is the main influence on the streamlines. In our derivation, we have used the parameter $\epsilon = \hat{H}_0/\hat{L} \ll 1$, but $\hat{H}_0/\hat{L}$ only needs to be interpreted as a scale for the slope of the interface and streamlines. Similarly, an underlying slowly varying topography can be included in such models without affecting validity, and there are many examples of lubrication/thin film models in 2 lateral directions. Thus, extension to stopes with sidewalls is feasible, and interesting in general.
<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 16.** Paste flow profile, (a) Paste 1, (b) Paste 2, (c) Paste 3, (d) Paste 4. See text for the associated (B,S) values.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 17.** Typical 3D mounding effect. Here observed on Paste 2: (left) soon after injection starts and (right) when the paste has significantly filled the stope. Mounding effect can be significant in the vicinity of the injector</figcaption></figure>

One direction we have not pursued is to make the rheology more complex. Paste tailings are varied and certainly many will exhibit thixotropic effects. However, the tailings are strongly sheared in the piping leading to the placement. Thus, aging starts from the injection time. However, thixotropic (aging) timescales are typically minutes-to-hours, which overlap with those of other processes relevant to the final quality of the filled material. In particular, the cement binder is likely to hydrate and any fluid loss to the surroundings will occur. In essence, a more complicated model is needed to incorporate chemical and transport (compositional) changes as well as rheological. Here we have simplified to the hydrodynamics of the filling process.

Future directions naturally come out of this study. In practice, it is frequent that various discharge inlet positions are used to solve the above obstruction effect. Hence good filling is often obtained by making several injections on already existing deposited layers of cemented paste. We plan to extend our model to be able to simulate different injection points as well as a thorough study of the effect of non constant bottom of the stope. The latter may result either from stope topography or from the setting of pre-existing cement layers. The latter requires a model of the changing cement rheology with time.

Full 3D simulations are another route to explore. They are needed for better modelling of the 3D mounding effect at injection but also the sidewall dynamics may have significant effects in some cases. However, note these 3D simulations can be computationally very expensive. We believe that for operational design aspects in the field, extensions of the 2D lubrication model as presented in this paper is probably the most useful tool due to computational affordability.
Acknowledgements

This research was initiated at the University of British Columbia during research visits of YS and PV YS gratefully acknowledges the China Scholarship Council Foundation for the financial support during PhD studies (Grant number: 201806460043), that made this possible; as well as the National Nature Science Foundation of China (52104124). PV gratefully acknowledges funding from the Pacific Institute for the Mathematical Sciences (CNRS IRL 3069) during his CNRS' délégation in 2018-2019, as well as the UMPA (CNRS UMR 5669) and CNRS Interdisciplinary program "Défi InFliti" (2018). PV also acknowledges funding from ANR (French National Research Agency) via the AAPG 2020 programme (Grant number: ANR-20-CE46-0006). IF acknowledges funding from NSERC via the discovery grants programme.

References

Achterberg, E.P., Braungardt, C., Morley, N.H., Elbaz-Poulichet, F., Leblanc, M., 1999. Impact of Los Frailes mine spill on riverine, estuarine and coastal waters in southern Spain. Water Research 33, 3387-3394. doi:https://doi.org/10.1016/S0043-1354(99)00282-1.

Ancey, C., 2007. Plasticity and geophysical flows: a review. Journal of non-Newtonian fluid mechanics 142, 4-35. doi:https://doi.org/10.1016/j.jnnfm.2006.05.005.

ASTM C143/C143M-20, 2020. Standard test method for slump of hydraulic-cement concrete. URL: https://www.astm.org/c0143_c0143m-20.

ASTM C305-14, 2014. Standard practice for mechanical mixing of hydraulic cement pastes and mortars of plastic consistency. doi:https://doi.org/10.1520/C0305-14.

ASTM D2487-17, 2017. Standard practice for classification of soils for engineering purposes (unified soil classification system). doi:https://doi.org/10.1520/D2487-17.

ASTM D4318-17e1, 2017. Standard test methods for liquid limit, plastic limit, and plasticity index of soils. doi:https://doi.org/10.1520/D4318-17E01.

Balmforth, N.J., Craster, R.V., Rust, A.C., Sassi, R., 2006. Viscoplastic flow over an inclined surface. Journal of Non-Newtonian Fluid Mechanics 139, 103-127. doi:https://doi.org/10.1016/j.jnnfm.2006.07.013.

Balmforth, N.J., Frigaard, I.A., Ovarlez, G., 2014. Yielding to Stress: Recent Developments in Viscoplastic Fluid Mechanics. Annual Review of Fluid Mechanics 46, 121-146. doi:https://doi.org/10.1146/annurev-fluid-010313-141424.

Edraki, M., Baumgartl, T., Manlapig, E., Bradshaw, D., Franks, D.M., Moran, C.J., 2014. Designing mine tailings for better environmental, social and economic outcomes: a review of alternative approaches. Journal of Cleaner Production 84, 411-420. doi:https://doi.org/10.1016/j.jclepro.2014.04.079.

Ercikdi, B., Baki, H., İzki, M., 2013. Effect of desliming of sulphide-rich mill tailings on the long-term strength of cemented paste backfill. Journal of Environmental Management 115, 5-13. doi: http://dx.doi.org/10.1016/j.jenvman.2012.11.014.

Fang, K., Fall, M., 2018. Effects of curing temperature on shear behaviour of cemented paste backfill-rock interface. International

Journal of Rock Mechanics and Mining Sciences 112, 184-192. doi:https://doi.org/10.1016/j.ijrmms.2018.10.024.

Gao, J., Fourie, A., 2018. Studies on thickened tailings deposition in flume tests using the computational fluid dynamics (CFD) method. Canadian Geotechnical Journal 56, 249-262. doi:10.1139/cgj-2017-0228. publisher: NRC Research Press.

GB175-2007, 2017. Common portland cement. URL: https://www.spc.org.cn/online/a5fff5d1e2da6c2491b1ee1e973ed8c2.html.

GB/T 8077-2012, 2012. Methods for testing uniformity of concrete admixture. URL: https://www.spc.org.cn/online/9d0723d4ac3db9bcfd91e10f1930c85e.html.

Glotov, V.E., Chlachula, J., Glotova, L.P., Little, E., 2018. Causes and environmental impact of the gold-tailings dam failure at Karamken, the Russian Far East. Engineering geology 245, 236-247. doi:https://doi.org/10.1016/j.enggeo.2018.08.012.

Han, D., Ferron, R.D., 2016. Influence of high mixing intensity on rheology, hydration, and microstructure of fresh state cement paste. Cement and Concrete Research 84, 95-106. doi:10.1016/j.cemconres.2016.03.004.

Henriquez, J., Simms, P., 2009. Dynamic imaging and modelling of multilayer deposition of gold paste tailings. Minerals Engineering 22, 128-139. doi:https://doi.org/10.1016/j.mineng.2008.05.010.

Koehler, E.P., Fowler, D.W., 2004. Development of a Portable rheometer for fresh Portland cement concrete. Report. International Center for Aggregates Research. URL: https://repositories.lib.utexas.edu/handle/2152/35338.

LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. doi:https://doi.org/10.1017/CBO9780511791253.

Liu, K.F., Mei, C.C., 1989. Slow spreading of a sheet of Bingham fluid on an inclined plane. Journal of Fluid Mechanics 207, 505-529. doi:https://doi.org/10.1017/S0022112089002685.

Mizani, S., He, X., Simms, P., 2013. Application of lubrication theory to modeling stack geometry of high density mine tailings. Journal of Non-Newtonian Fluid Mechanics 198, 59-70. doi:https://doi.org/10.1016/j.jnnfm.2013.03.002.

Panchal, S., Deb, D., Sreenivas, T., 2018. Variability in rheology of cemented paste backfill with hydration age, binder and superplasticizer dosages. Advanced Powder Technology 29, 2211-2220. doi:10.1016/j.apt.2018.06.005.

Roussel, N., Coussot, P., 2005. "Fifty-cent rheometer" for yield stress measurements: From slump to spreading flow. Journal of Rheology 49, 705-718. URL: https://sor.scitation.org/doi/abs/10.1122/1.1879041, doi:10.1122/1.1879041.

Sun, W., Wang, H., Hou, K., 2018. Control of waste rock-tailings paste backfill for active mining subsidence areas. Journal of Cleaner Production 171, 567-579. doi:https://doi.org/10.1016/j.jclepro.2017.09.253.

Tan, Z., Bernal, S.A., Provis, J.L., 2017. Reproducible mini-slump test procedure for measuring the yield stress of cementitious pastes. Materials and Structures 50, 235-246. doi: http://dx.doi.org/10.1617/s11527-017-1103-x.

Wu, A., Wang, Y., Wang, H., 2015. Estimation model for yield stress of fresh uncemented thickened tailings: Coupled effects of true solid density, bulk density, and solid concentration. International Journal of Mineral Processing 143, 117-124. doi:https://doi.org/10.1016/j.minpro.2015.09.010.

Yilmaz, E., Fall, M., 2017. Paste tailings management. Springer. doi:https://doi.org/10.1007/978-3-319-39682-8.