Hostname: page-component-54dcc4c588-tfzs5 Total loading time: 0 Render date: 2025-09-12T02:08:15.777Z Has data issue: false hasContentIssue false

Viscoplastic slumps supported by a barrier

Published online by Cambridge University Press:  12 August 2025

Nitay Ben-Shachar
Affiliation:
School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Australia Department of Physics, California Institute of Technology, Pasadena, CA 91125, USA
Douglas R. Brumley
Affiliation:
School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Australia
Andrew J. Hogg
Affiliation:
School of Mathematics, University of Bristol, Woodland Road, Bristol BS8 1UG, UK
Edward M. Hinton*
Affiliation:
School of Mathematics and Statistics, The University of Melbourne, Parkville, VIC 3010, Australia
*
Corresponding author: Edward M. Hinton, edward.hinton@unimelb.edu.au

Abstract

The shape of a free-surface slump of viscoplastic material supported by an oblique barrier on an inclined plane is investigated theoretically and experimentally. The barrier is sufficiently tall that it is not surmounted by the viscoplastic fluid, and a focus of this study is the largest volume of rigid viscoplastic fluid that can be supported upstream of it. A lubrication model is integrated numerically to determine the transient flow as the maximal rigid shape is approached. Away from the region supported by the barrier, the viscoplastic layer attains a uniform thickness in which the gravitational stresses are in balance with the yield stress of the material. However, closer to the barrier, the layer thickens and the barrier bears the additional gravitational loading. An exact solution for the rigid shape of the viscoplastic material is constructed from the steady force balance and computed by integrating Charpit’s equations along characteristics that emanate from the barrier wall. The characteristics represent the late-time streamlines of the flow as it approaches the rigid shape. The exact solution depends on a single dimensionless group, which incorporates the slope inclination, the barrier width and the fluid’s yield stress. It is shown that the shape is insensitive to the transient flow from which it originates. The force exerted by the slump is calculated for different barrier shapes. The results of new laboratory experiments are reported; these show that although convergence to the final rigid state is slow, there is good agreement with the experimental measurements at long times.

Information

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2025. Published by Cambridge University Press

1. Introduction

Obstructions that hold up viscoplastic material occur in environmental and industrial contexts, including mine tailing dams, porous networks such as forests, and defences against hazardous flows such as lava and debris (e.g. Barberi et al. Reference Barberi, Brondi, Carapezza, Cavarra and Murgia2003; Sigtryggsdottir et al. Reference Sigtryggsdóttir, Hrafnsdóttir, Steingrímsson and Gudmundsson2025). Key engineering questions include: what is the maximum height or volume of material that is held up by a barrier, what is the force on the barrier, and can properties of the fluid be inferred from the shape of the supported material (Pedersen et al. Reference Pedersen, Pfeffer, Barsotti, Tarquini, de’Michieli Vitturi, Óladóttir and Prastarson2023)? Indeed, it has been shown that markers of the maximum lava flow height on the surface of trees can be used to effectively recreate aspects of the flow history (Chevrel et al. Reference Chevrel, Harris, Ajas, Biren, Gurioli and Calabrò2019).

Figure 1. Schematics of viscoplastic flow around an obstruction. (a) Steady flow. (b) Shape of the stationary slump following termination of the upstream source. (c) Cross-section of the steady flow along the centreline. (d) Photo of laboratory experiment 1; see § 5 (obstruction has extent 200 mm).

Barriers to lava flows have been used or proposed at a variety of sites over the past 70 years with mixed success, including in Hawaii and at Etna, and most recently, to protect the town of Grindavik on the Reykjanes peninsula in Iceland (Moore Reference Moore1982; Barberi et al. Reference Barberi, Brondi, Carapezza, Cavarra and Murgia2003; Bjarnason & Bubola Reference Bjarnason and Bubola2024; Sigtryggsdottir et al. Reference Sigtryggsdóttir, Hrafnsdóttir, Steingrímsson and Gudmundsson2025). Improvements in the effectiveness of this technology require a comprehensive understanding of the physical mechanisms that control the interaction of lava flows with barriers (Dietterich et al. Reference Dietterich, Cashman, Rust and Lev2015).

In this paper, we use an analogue theoretical and laboratory model to investigate how lava is held up by a barrier (see figure 1). Lava is represented by an isothermal viscoplastic fluid (e.g. Griffiths Reference Griffiths2000). We calculate the shape of the remnant material supported by the barrier, and furthermore, we demonstrate that the model may be used to approximate the yield stress from observations of the supported fluid (cf. Hinton Reference Hinton2022).

In contrast to viscous fluids, viscoplastics have a yield stress that must be exceeded for motion to occur, and this bounds their spreading. Even when subjected to a gravitational body force, they come to rest at some finite slumped shape, whose thickness is related to the yield stress. In the absence of an obstruction, the final slumped shape of a viscoplastic has been widely studied on both horizontal and inclined planes (Balmforth et al. Reference Balmforth, Craster and Sassi2002, Reference Balmforth, Craster, Perona, Rust and Sassi2007; Dubash et al. Reference Dubash, Balmforth, Slim and Cochard2009). These studies demonstrated that the slump shape can be used as a simple rheometer to infer the yield stress (Roussel & Coussot Reference Roussel and Coussot2005).

The majority of these investigations exploit a shallow approximation to simplify the governing equations. For viscoplastic fluids, some care is required in applying the shallow model to free-surface flows (Balmforth & Craster Reference Balmforth and Craster1999). Generally, the flow is partitioned into two regions: a yielded region near the substrate, where the deviatoric stress is greater than the yield stress, and a region adjacent to the free surface, where the stress does not exceed the yield stress and the material moves as a pseudo-plug; see Balmforth & Craster (Reference Balmforth and Craster1999). The boundary between the two regions is called the ‘yield surface’; as the material thins and approaches its final slumped shape, the yield surface approaches the substrate, and the driving gravitational stresses are balanced by the yield stress.

Shallow free-surface models of viscoplastics have been used to investigate steady and transient flows around obstructions and topography (Ionescu Reference Ionescu2013; Hinton & Hogg Reference Hinton and Hogg2022; Hinton, Hewitt & Hogg Reference Hinton, Hewitt and Hogg2023). It was shown that ‘dead’ zones – in which the entire thickness of the layer is rigid – occur in the vicinity of the upstream stagnation point. In general, the size of the dead zone is an increasing function of the yield stress and a decreasing function of the upstream flow rate. The present study focuses on the geometry of the arrested state; it is effectively the limit of a zero upstream flow rate, and the dead zone comes to encompass the entirety of the upstream face of the obstruction.

Viscoplastic slumps supported by barriers have some similarities with their granular counterparts, and both can form arrested wedge-shaped deposits (Cui & Gray Reference Cui and Gray2013; Tregaskis et al. Reference Tregaskis, Johnson, Cui and Gray2022). Both materials exhibit a yield criterion such that they deform only when the deviatoric shear stress exceeds a threshold. For viscoplastic materials, this threshold, the yield stress, is an intrinsic material property, whereas for granular material, it is proportional to the normal stress. Thus the yielding condition is modelled differently and manifests in rather different forms of governing equations. However, it is shown in Appendix C that when the yield stress of the viscoplastic material is much greater than the gravitational weight of the layer per unit area (a condition that is made precise in Appendix C), the leading-order governing equation for the arrested viscoplastic becomes identical to the granular flow equation apart from at the edges of the arrested material.

The paper is structured as follows. The model of Hinton et al. (Reference Hinton, Hewitt and Hogg2023) for steady flow around a barrier is summarised in § 1.1 to define governing parameters and equations. In § 2, we analyse the evolution of this flow following termination of the upstream source. In § 3, an exact expression for the final rigid shape is given. Features of the rigid shape and the force it exerts on the barrier are discussed in § 4. Results from new laboratory experiments are presented and compared to the theory in § 5. Concluding remarks are given in § 6. Some ancillary details are included in the appendices. The numerical method is described in Appendix A, the integration of Charpit’s equations in Appendix B, the asymptotic behaviour at low and high yield stress in Appendix C, and the generalisation of the analysis to non-rectilinear obstructions in Appendix D.

1.1. Steady viscoplastic flow around an obstruction

Hinton et al. (Reference Hinton, Hewitt and Hogg2023) analysed the steady, shallow and inertialess free-surface flow of a Bingham fluid down an inclined plane around an obstruction; see figure 1(a). Here, we extend this formalism to a Herschel–Bulkley fluid of power-law index $n$ , yield stress $\tau _Y$ , and consistency $k$ . The plane is inclined at an angle $\beta$ with respect to horizontal. The $\hat {x}$ axis is directed downslope, the $\hat {y}$ axis is directed cross-slope, and the $\hat {z}$ axis is perpendicular to the plane. The fluid has thickness $\hat {h}(\hat {x},\hat {y})$ and density $\rho$ . The obstruction has square cross-section and lies in $(\hat {x},\hat {y})\in [-L,L]\times [-L,L]$ . The flow is supplied by an upstream source of strength $\hat {Q}_\infty$ per unit length in the $\hat {y}$ direction.

The variables are non-dimensionalised via

(1.1) \begin{equation} (x,y) =(\hat {x},\hat {y})/L, \quad {(z,H) =(\hat {z},\hat {h})/\hat {H}_\infty , } \end{equation}

where $\hat {H}_\infty$ is the critical rigid height of viscoplastic material on an inclined plane,

(1.2) \begin{align} \hat {H}_\infty = \frac {\tau _Y}{\rho g \sin \beta } . \end{align}

The flux (the depth-integrated velocity) in the $x$ and $y$ directions was derived following Balmforth, Craster & Sassi (Reference Balmforth, Craster and Sassi2002) as

(1.3) \begin{align} \boldsymbol{q}=\frac {3nY(s Y)^{1/n}\left [(2n+1)H-nY)\right ]}{s(n+1)(2n+1)} \left (1-b\frac {\partial {\textit{H}}}{\partial x},-b \frac {\partial {\textit{H}}}{\partial y}\right )\!, \end{align}

where

(1.4) \begin{align} s=\sqrt {\left (1-b\frac {\partial {\textit{H}}}{\partial x}\right )^2+\left (b\frac {\partial {\textit{H}}}{\partial y}\right )^2} , \end{align}

and the height of the yield surface is

(1.5) \begin{align} Y(x,y)=\text{max}\left (0,H-\frac {1}{\sqrt {\left (1-b\dfrac {\partial {\textit{H}}}{\partial x}\right )^2+\left (b\dfrac {\partial {\textit{H}}}{\partial y}\right )^2}}\right )\! . \end{align}

The flowing solution is governed by two dimensionless parameters: the modified Bingham number

(1.6) \begin{align} b=\frac {\tau _Y}{\rho g L \sin \beta \tan \beta }=\frac {\hat {H}_\infty }{L\tan \beta } , \end{align}

and the dimensionless upstream flux $q_\infty$ . The dimensionless group $b$ is sometimes referred to as the plastic number (e.g. Thompson & Soares Reference Thompson and Soares2016), but here we retain the term ‘Bingham number’ for consistency with the literature on viscoplastic gravity currents (Balmforth & Craster Reference Balmforth and Craster1999; Balmforth et al. Reference Balmforth, Craster, Perona, Rust and Sassi2007; Hogg & Matson Reference Hogg and Matson2009).

Our choice of scaling for the flow thickness in (1.1) differs from the scales used by Hinton et al. (Reference Hinton, Hewitt and Hogg2023), who non-dimensionalised the flow thickness with the thickness of Newtonian flow (this is equivalent to scaling the upstream flux by its far-field value). The scales employed here are chosen for ease of analysis of the stationary slump obtained once the upstream source is terminated.

The flowing solution of Hinton et al. (Reference Hinton, Hewitt and Hogg2023) is controlled by two parameters, $B$ and $\mathcal{L}$ (defined in their (2.5a,b)), and here $b=B/\mathcal{L}$ . In contrast to the flowing solution developed by Hinton et al. (Reference Hinton, Hewitt and Hogg2023), the model for the arrested state developed here features only the single parameter $b$ .

The steady flow is governed by

(1.7) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} = 0, \end{align}

with no-flux applied on the boundaries of the obstruction,

(1.8) \begin{align} \boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{n} = 0, \end{align}

where $\boldsymbol{n}$ is the boundary-normal vector. Far upstream (and far downstream), the dimensionless flux satisfies $\boldsymbol{q} \to (q_\infty ,0)$ , whilst $H \to H_\infty \geqslant 1$ , which reflects a constant flux in the far field. These are related via

(1.9) \begin{align} q_\infty = \frac {3n(H_\infty -1)^{1+1/n}\left ((n+1)H_\infty +n\right )}{(n+1)(2n+1)}, \end{align}

which is obtained directly from (1.3). Solutions for the steady free surface from Hinton et al. (Reference Hinton, Hewitt and Hogg2023) are shown in figures 2(a) and 2(b) for two sets of parameter values.

2. Unsteady model and switching off the source

We analyse the unsteady flow that arises when the upstream source is terminated. The dimensionless variables are given by (1.1), and dimensionless time is given by

(2.1) \begin{equation} { t = \frac {\rho g \hat {H}_\infty ^2 \sin \beta }{3\mu L}\hat {t}}, \end{equation}

where

(2.2) \begin{align} {\mu =k^{1/n} \left (\frac {\rho g \hat {H}_\infty \sin \beta }{3}\right )^{1-1/n}} \end{align}

is a viscosity scale, and $t=0$ corresponds to the switch-off time. The evolution of the flow thickness $H(x,y,t)$ is governed by

(2.3) \begin{equation} \frac {\partial {\textit{H}}}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} = 0. \end{equation}

The initial condition is given by a steady flowing solution computed using the method described by Hinton et al. (Reference Hinton, Hewitt and Hogg2023), with examples plotted in figure 2(a–c).

Figure 2. Steady flowing solutions for $b=1$ and: (a) Bingham fluid, $H_\infty =2.53$ (corresponding to $q_\infty =7.1$ ); (b) Bingham fluid, $H_\infty =1.68$ ( $q_\infty =1$ ); (c) Herschel–Bulkley fluid, $n=0.5$ , $H_\infty =1.30$ ( $q_\infty =0.03$ ); see (1.3). When the flux source is switched off, all of these converge to the same final stationary slump, shown in (d). Note the different scale bars used across the panels.

Equation (2.3) is numerically integrated in time; for details of the numerical method, see Appendix A. In the computation, the source is located relatively far upstream at $x=-6$ , so at $t=0$ , the boundary condition at the left-hand end of the domain is adjusted to

(2.4) \begin{equation} H =1 \quad \text{at} \ x =-6. \end{equation}

Note that a different location for the upstream source simply leads to a time offset in the ensuing dynamics after the switch off (the reduction in $H$ simply takes longer to propagate downstream towards the barrier). The no-flux boundary condition on the obstruction (1.8) is unchanged.

Fluid drains downstream, and the flow thickness converges to a rigid stationary shape everywhere; see figures 2 and 3.The final rigid shape is approached asymptotically with the following late-time scalings when $t \gg 1$ (Matson & Hogg Reference Matson and Hogg2007; Hogg & Matson Reference Hogg and Matson2009):

(2.5) \begin{align} {H (x,y,t) = H_{\textit{final}}(x,y) + t^{-n}\, H_1(x,y) +\cdots , \quad Y(x,y,t)= t^{-n}\, Y_1(x,y) +\cdots ,} \end{align}

which are verified numerically for our problem in figure 4. Here, $H_{\textit{final}}(x,y)$ is the final height of the free surface, while $H_1(x,y)$ and $Y_1(x,y)$ dictate the late-time approach of the height of the free surface and yield surface, respectively, to their final static shapes. Henceforth, we drop the subscript ‘final’ for brevity. Figure 4 shows that fluid near the upstream corners of the barrier takes longer to drain away and approach the arrested height than elsewhere. The yield surface exhibits a maximum at the upstream corners, and is also significant on the sides of the barrier and at the downstream corners. This is because fluid draining from upstream is diverted around the obstruction and thus flows past the upstream corners, giving rise to a large flux at these locations (see figure 5). The increased flux at the corners requires a higher yield surface. A similar phenomenon occurs at the downstream corners, where fluid flowing around the sharp corner induces a large flux, and therefore a relatively high yield surface, but to a lesser extent.

Figure 3. Late-time slumping of a Bingham fluid. (a,b) The height of the free surface (dashed black lines) at the upstream wall of the obstruction ( $x=-1$ ) at times $t\in \{\textrm{e}^1,\textrm{e}^2,\textrm{e}^3,\textrm{e}^4,\textrm{e}^5,\textrm{e}^6\}$ after the upstream flux (at $x=-6$ ) is terminated. (c,d) Corresponding results along the centreline ( $y=0$ ). In (a) and (c), parameter values are $b=1,\ q_\infty =1$ , whilst for (b) and (d), $b=1,\ q_{\infty }=1/2$ . The final shape is shown with a blue line (from § 3).

Figure 4. Late-time decay of (a,b,c) a Bingham fluid to the final state, and (d,e,f) a Herschel–Bulkley fluid with $n=1/4$ . (a,d) Transient evolution of the relative flow thickness at four points along the centreline ( $y=0$ ); these show the late-time decay rate $\sim 1/t^n$ . (b,e) The relative difference between the fluid thickness at $t=20$ and the final rigid thickness ( $t \to \infty$ ). (c, f) The relative height of the yield surface at $t=20$ . Parameter values are $b=1$ and $H_\infty =1.68$ .

Figure 5. Streamlines of the flow (white lines) and magnitude of the flux (colour bar) for a Bingham fluid at $t=20$ . Parameter values are $b=1$ and $H_\infty =1.68$ .

3. Static shape

Figure 2(a–c) show the thickness of the steady flow for different values of $q_\infty$ and $n$ . When the source is terminated, these three cases all converge to the same rigid shape, suggesting that it depends only on $b$ , and specifically, is independent of the power-law exponent $n$ and the upstream source flux $q_\infty$ . In this section, we derive an exact expression for the rigid shape.

In the far field, the flow thickness converges to the critical rigid height of viscoplastic material on an inclined plane, which furnishes

(3.1) \begin{align} H \to 1 \quad \text{as} \ x^2 +y^2 \to \infty . \end{align}

The rigid shape is found by evaluating (2.3) at late times with $Y \to 0$ . Thus the maximal rigid shape is found to satisfy (Balmforth et al. Reference Balmforth, Craster and Sassi2002)

(3.2) \begin{align} \left (1-b\frac {\partial {\textit{H}}}{\partial x}\right )^2 + \left (b \frac {\partial {\textit{H}}}{\partial y}\right )^2 = \frac {1}{H^2} , \end{align}

where $b$ is the single dimensionless parameter that characterises the shape of the rigid slump. The no-flux boundary condition on the obstruction walls, which applies at all times, imposes

(3.3) \begin{align} \left (1-b\frac {\partial {\textit{H}}}{\partial x}, -b\frac {\partial {\textit{H}}}{\partial y}\right ) \boldsymbol{\cdot }\boldsymbol{n} = 0, \end{align}

which is applied at $x=\pm 1,\ {-1}\leqslant y\leqslant 1$ and $y=\pm 1,\ {-1}\leqslant x\leqslant 1$ . The slump is determined maximal when it is ‘on the verge’ of yielding; any further addition of material to the slump will induce a flow that carries the added fluid downhill. The governing partial differential equation (3.2) and boundary conditions (3.1), (3.3) feature only $b$ , and this is corroborated by the different flows in figure 2(a,b,c) all adjusting to the same arrested state.

3.1. Shape of the free surface at the upstream wall

We parametrise the upstream wall of the obstruction by $y=s$ with $-1 \leqslant s \leqslant 1$ , $x=-1$ and $H(-1,s) = H_0(s)$ , which is to be determined. Combining (3.2) and (3.3) on the upstream wall furnishes

(3.4) \begin{align} {\textit{bH}}_0 \frac {\mathrm{d} H_0}{\mathrm{d}s} = -\text{sgn}(s). \end{align}

The boundary condition in the far field, (3.1), imposes the following condition at the corners of the obstruction:

(3.5) \begin{align} H_0(\pm 1) = 1 . \end{align}

This condition is required for a self-consistent solution. Specifically, one cannot satisfy (3.1) if the height at the obstruction corners exceeds the far-field height. Equation (3.4) is integrated to obtain

(3.6) \begin{align} H_0(s)&=\sqrt {2b^{-1}(1-|s|)+1}. \end{align}

This shape is shown as a blue line in figures 3(a) and 3(b), where it is compared to the numerical results for the unsteady flow. The gradient of $H_0$ is discontinuous at $s=0$ .

3.2. Charpit’s method

Equation (3.6) provides the boundary conditions to solve (3.2) in the region directly upstream of the obstruction where the rigid material is supported. The solution is obtained using Charpit’s method for first-order nonlinear partial differential equations. Defining

(3.7) \begin{align} p=\frac {\partial {\textit{H}}}{\partial x}, \quad q=\frac {\partial {\textit{H}}}{\partial y}, \end{align}

Charpit’s equations transform (3.2) into five ordinary differential equations along the characteristics (Ockendon et al. Reference Ockendon, Howison, Lacey and Movchan2003):

(3.8a) \begin{align} \dot {x} &= -2b(1-bp) , \end{align}
(3.8b) \begin{align} \dot {y} &= 2b^2 q , \end{align}
(3.8c) \begin{align} \dot {p} &= -2p H^{-3} , \end{align}
(3.8d) \begin{align} \dot {q} &= -2q H^{-3} , \end{align}
(3.8e) \begin{align} \dot {H} &= 2b^2q^2-2bp(1-bp) , \end{align}

where a dot indicates differentiation with respect to $\tau$ , which parametrises the characteristics. The characteristics (3.8a )–(3.8e ) emanate from the upstream wall of the obstruction where the boundary data from (3.6) are specified. The forms of the characteristics are qualitatively different depending on whether they emanate from the face of the upstream wall ( $-1\lt y\lt 1$ ) or from the corners ( $y=\pm\, 1$ ), so the analysis is separated accordingly.

3.2.1. Characteristics emanating from the upstream wall ( $-1\lt y\lt 1$ )

By combining (3.6) and (3.2), we obtain the following boundary data on the upstream wall for the characteristics:

(3.9a) \begin{align} x_0(s) &=-1 , \end{align}
(3.9b) \begin{align} y_0(s) &= s , \end{align}
(3.9c) \begin{align} p_0(s) &= \frac {1}{b} , \end{align}
(3.9d) \begin{align} q_0(s) &= -\frac {\text{sgn}(s)}{b \sqrt {2b^{-1}(1-|s|)+1}} , \end{align}
(3.9e) \begin{align} H_0(s) &= \sqrt {2b^{-1}(1-|s|)+1} , \end{align}

where subscript ‘0’ indicates evaluation at $\tau =0$ . Equations (3.8a )–(3.8e ) are integrated subject to the boundary data (3.9a )–(3.9e ). This furnishes part of the exact solution, given by (full details are reported in Appendix B)

(3.10a) \begin{align} x+1 &= b \left [H-\sqrt {2b^{-1}(1-|s|)+1}\right ] \nonumber \\&\quad -\, b \log \left [\frac {\pm b\left (\sqrt {2b^{-1}(1-|s|)+1}-1\right )-y+s}{(H-1)b}\right ]\! , \end{align}
(3.10b) \begin{align} s &= y\pm b\sqrt {2b^{-1}(1-|y|)+1-H^2} . \end{align}

This solution covers the region of the $(x,y)$ plane penetrated by characteristics emanating from the upstream wall, which is given by

(3.11) \begin{equation} -1\lt \, y\lt 1 \quad \text{and} \quad b \left (A(y)-1 \right ) -b \log \left (\frac {1-|y|}{b(A(y)-1)}\right ) \leqslant \, x+1\leqslant 0, \end{equation}

with

(3.12) \begin{align} A(y) \equiv H(y,s=1)= \sqrt {1+2b^{-1}(1-|y|) - b^{-2}(1-|y|)^2}. \end{align}

The characteristics and the region penetrated are shown by black lines in figure 6. We note that $A(0) \geqslant 1$ if $b\gt 1/2$ . Thus for $b \leqslant 1/2$ , the characteristics encompass the entirety of the centreline, $y=0$ , $x \in (-\infty ,-1]$ . The solution beyond this region is obtained by considering characteristics that emanate from the upstream corners.

Figure 6. (a) Characteristics for the thickness of the rigid shape supported by the obstruction for $b=1/2$ . (b) Characteristic projections in the $(x,y)$ plane corresponding to (a). (c) Characteristic projections for $b=1$ . (d) Characteristic projections for $b=2$ . The solid black lines represent characteristics emanating from the upstream wall, whilst the red dashed lines represent those emanating from the corners. See also figure 7.

3.2.2. Characteristics from the upstream corners

At the upstream corners, the no-flux condition (3.3) is applied at all incoming angles to the corner. That is, we parametrise the boundary normal vector at the upstream corners by

(3.13) \begin{align} \boldsymbol{n}=(\sin \theta , \cos \theta ), \end{align}

with $\theta \in [-\unicode{x03C0} /2,0]$ and $\theta \in [0, \unicode{x03C0} /2]$ for the corners at $y=1$ and $-1$ , respectively, where $\theta =0$ corresponds to the positive $x$ axis. Combining (3.3) and (3.2) gives the following boundary data for characteristics emanating from the corners:

(3.14a) \begin{align} x_0(\theta ) &= -1 , \end{align}
(3.14b) \begin{align} y_0(\theta ) &= \pm 1 , \end{align}
(3.14c) \begin{align} p_0(\theta ) &= \frac {1+\sin \theta }{b} , \end{align}
(3.14d) \begin{align} q_0(\theta ) &= \mp \frac {\cos \theta }{b} , \end{align}
(3.14e) \begin{align} H_0(\theta ) &= 1 , \end{align}

where the plus and minus signs are used to distinguish between the two upstream corners. The solution for the shape covered by these characteristics is (see Appendix B)

(3.15) \begin{equation} x+1 = b (H-1) - b \log \left [ \frac {H^2 -1 + b^{-2}(1-|y|)^2}{2(H-1)} \right ]\! , \end{equation}

valid in the domain

(3.16) \begin{equation} -1 \leqslant y \leqslant 1 \quad \text{and} \quad x+1 \lt b \left (A(y)-1 \right )-b \log \left (\frac {1-|y|}{b(A(y)-1)}\right )\!. \end{equation}

The characteristics and the region penetrated are shown by red dashed lines in figure 6.

Outside the domains dictated by (3.11) and (3.16), i.e. for $|y|\gt 1$ , the solution for the fluid thickness is $H=1$ .

Figure 7. The exact solution for the stationary slump ((3.10a ), (3.10b ) and (3.15)) for (a) $b=1/4$ , (b) $b=1/2$ , (c) $b=1$ , (d) $b=2$ .

Figure 8. (a) The dimensionless maximum slump height as function of $b$ . (b) The dimensional maximum height as a function of the yield stress $\tau _Y$ for $\rho g = 10^4$ $\mathrm{kg\,m^{-2}\,s^{-2}}$ , $L=0.1$ $\mathrm{m}$ and $\beta =5 ^{\circ }$ . (c) Dimensionless excess volume of the slump ((4.3), solid blue line), compared with the small and large $b$ asymptotic solutions (see (4.4a )–(4.4b )). (d) Dimensional volume as a function of $\tau _Y$ for the above-mentioned parameter values. The scaled maximum height and volume obtained in experiments 1 and 3 (see § 5) are shown in (a) and (c), where the error bars are smaller than the extent of the markers.

4. Discussion

Equations (3.10a ), (3.10b ) and (3.15) constitute complete and exact expressions for the final slump shape in terms of $b$ . Examples of the shape for four values of $b$ are shown in figure 7. It should be noted that a steeper slope corresponds to smaller $b$ , which is associated with a thicker and more localised slump. From the exact solution, many physical aspects of the slump can be computed. The maximum thickness of the slump, attained at $(x,y)=(-1,0)$ , is given by

(4.1) \begin{equation} H_{\textit{max}} = \sqrt {1+2/b}, \end{equation}

which is shown in figure 8(a). In dimensional terms, the maximum thickness is

(4.2) \begin{equation} \hat {H}_{\textit{max}} = \sqrt {\frac {\tau _Y^2}{(\rho g \sin \beta )^2}+\frac {2 \tau _Y L}{\rho g \cos \beta }}, \end{equation}

which is plotted in figure 8(b) for an example case with $\rho g = 10^4$ $\mathrm{kg\,m^{-2}\,s^{-2}}$ , $L=0.1$ m and $\beta =5 ^{\circ }$ , which is relevant to our laboratory experiments; see § 5. The dimensionless maximum thickness is a monotonically decreasing function of $b$ , whilst the dimensional maximum thickness is a monotonically increasing function of the yield stress $\tau _Y$ . This difference arises because the dimensionless thickness has been scaled by the thickness in the far field, which is linearly proportional to the yield stress. Thus although the slump thickness is larger for a fluid with a higher yield stress, the ratio of the slump thickness to the far-field thickness is smaller. A similar feature occurs with the excess volume, discussed below.

The excess volume of the slump, $V$ , is defined as the additional volume that is supported by the obstruction as compared to the case with no obstruction for which $H\equiv 1$ , i.e.

(4.3) \begin{align} V = \int _{-\infty }^{-1} \int _{-1}^1 \left ( H(x,y) -1 \right )\, \mathrm{d}y\, \mathrm{d}x. \end{align}

The excess volume is plotted as a function of $b$ in figure 8(c). The dimensional excess volume is shown in figure 8(d) (with the same parameter values as in figure 8 b) .

The parameter $b$ , which is proportional to the yield stress, represents the ratio of the far-field rigid thickness to the obstruction width and the slope gradient (see (1.6)). The regime of small $b$ corresponds to a relatively small far-field thickness (or wide obstruction), and large $b$ corresponds to a relatively large far-field thickness. For $b \gg 1$ , the governing equations away from the margins of the arrested slump become identical to those of a static granular pile on an empty plane (see § 7.1 of Tregaskis et al. Reference Tregaskis, Johnson, Cui and Gray2022); see also the calculations in § C.2 (though it should be noted that the static granular pile on an empty plane has $H \to 0$ away from the barrier, whereas here we have $H\to 1$ ). In this regime of a relatively large Bingham number, the static solution calculated here is also a good approximation for the height of the free surface of steady flowing material from a constant flux source. This steady flow driven by a constant-flux line source was analysed by Hinton et al. (Reference Hinton, Hewitt and Hogg2023), and in their notation this regime corresponds to $B\gg 1$ and $\mathcal{L}=O(1)$ .

For small and large $b$ , (4.3) admits the respective asymptotic limits

(4.4a) \begin{align} V &\sim 1 -\left (\frac {8-2\unicode{x03C0} }{3\sqrt {2}}\right ) \sqrt {b} , & b \ll 1 , \end{align}
(4.4b) \begin{align} V &\sim \frac {1}{18b} \left (5 + 6 \log (2b) \right )\!, & b \gg 1 ; \end{align}

for details, see Appendix C. In dimensional form, this is given by

(4.5a) \begin{align} \hat {V} &\sim \frac {L^2 \tau _Y}{\rho g \sin \beta } -\left (\frac {8-2\unicode{x03C0} }{3\sqrt {2}}\right ) \left (\frac {L \tau _Y}{\rho g \sin \beta } \right )^{3/2} (\tan \beta )^{-1/2} ,& b \ll 1 , \end{align}
(4.5b) \begin{align} \hat {V} &\sim \frac {L^3 \tan \beta }{18} \left (5 + 6 \log (2b) \right )\!, & b \gg 1. \end{align}

Interestingly, the dimensional excess volume exhibits different scaling laws with respect to the yield stress $\tau _Y$ in the regimes of small and large $b$ . In the relatively high yield stress regime ( $b\gg 1$ ), the excess volume is only weakly dependent on the yield stress through the $\log (b)$ term, and instead the volume depends primarily on the geometry of the slope and the barrier. Similar dependencies arise for a cylindrical obstruction; see Appendix D.

One other interesting observation concerning the rigid solution is that the characteristic projections (figure 6) are exactly the streamlines of the late-time fluid flow as it approaches the rigid state; this is shown in figure 9 for $b=1$ . The analogy arises because the characteristics follow $(\dot {x},\dot {y}) \propto (1-b\, \partial \textit{H}/\partial x, - b\, \partial \textit{H}/\partial y)$ , which is the flow direction. The characteristics thus emanate tangentially to the obstacle, and data for the rigid solution are carried in the opposite direction to the flow.

Figure 9. The streamlines of the numerically calculated late-time flow (solid blue lines) compared with the characteristic projections from Charpit’s method (dashed black lines), for $b=1$ . The obstruction is shown in grey.

The stationary free-surface thickness exhibits a seam along the centreline, $y=0$ . Here, the characteristics of the rigid solution collide and terminate. This corresponds to the late-time flow of the slumping fluid bifurcating at the centreline.

4.1. Force on the barrier

The dimensional hydrostatic force $(\hat {F},0)$ exerted on the barrier by the arrested slump is given by the integral of the pressure over the upstream face,

(4.6) \begin{equation} \hat {F}=\int _{-L}^L \int _0^{\hat {h}(-1,\hat {y})} \rho g \big(\hat {h}(-1,\hat {y})-\hat {z} \big) \cos \beta \, \mathrm{d} \hat {z}\, \mathrm{d} \hat {y}, \end{equation}

and using the solution for the shape at this face, (3.6), we obtain

(4.7) \begin{equation} \hat {F}= \left (\frac {\tau _Y}{\rho g \sin \beta } + L \tan \beta \right ) \tau _Y L \cot \beta . \end{equation}

4.2. Insensitivity to the upstream source

The stationary slump shape calculated in § 3 is insensitive to the upstream flux source provided that the flow exceeded the slump height. Here, we consider the flow generated by an off-centre vent source located at $-2\leqslant y\leqslant -1/2$ , $x=-4$ supplying constant flux (figure 10 a). Using the techniques described in § 2 and initialising the fluid at zero height i.e. $H(x,y,t=0)=0$ , the steady flow and the flow following the switch off are calculated; see figure 10. The steady flow exhibits an asymmetric profile. However, when the source is switched off, the stationary shape upstream of the obstruction matches the analytical calculations of § 3.

Figure 10. (a,c) Flowing and (b,d) stationary solutions for Bingham fluid ( $n=1$ ) with $b=1/4,\ q_\infty =1$ , for (a,b) off-centre vent source, (c,d) infinite line source.

4.3. Downstream behaviour

The shape of the free surface of rigid material downstream of the obstruction depends on the initial conditions of the flow prior to switch off as well as the upstream flux source and the viscous rheology of the fluid. This is in contrast to the upstream shape, which is independent of these parameters (see discussion above). The dependence of the downstream shape arises because parts of the downstream region are inaccessible to the flow. Thus if the fluid is initialised in this region to exhibit a static slump, then this shape will be maintained independently of the fluid flow impinging from the flux source; e.g. see unwetted region in figures 10(b) and 10(d). Moreover, the region inaccessible to the upstream flow depends on the strength of the flux source. The stationary downstream shape is altered at higher values of the upstream flux because more fluid invades the region directly downstream, and subsequently cannot flow away; see also Hinton & Hogg (Reference Hinton and Hogg2022), who observed similar behaviour upstream and downstream of a topographic mound. The non-universality of the downstream shape is exemplified by the different static shapes obtained for an infinite line-flux source and a vent source that are initialised to zero height prior to the the flux source being turned on; compare figures 10(b) and 10(d). We illustrate this explicitly by plotting the shape of the free surface at $x=1$ (i.e. at the downstream face of the obstruction) in figure 11, at different times after the flux source is terminated for both an infinite line flux source and a vent. This shows that for the vent source (figure 11 b), less fluid invades the downstream region near $y=1$ than near $y=-1$ , as compared with the line source (figure 11 a), where the downstream region invaded by the fluid is symmetric in $y$ . This exemplifies the non-universality of the downstream rigid shape.

Figure 11. Fluid thickness at $x=1$ (i.e. the downstream face of the obstruction), at various times after the source has been switched off: (a) infinite line source, (b) vent source, as described in figure 10. The dashed black lines indicate the edges of the obstruction’s downstream face.

5. Laboratory experiments

We report experiments on an inclined plane of width 40 cm and length 50 cm held at a shallow angle to the horizontal. A cuboid obstruction of total width 20 cm, or a cylindrical obstruction with diameter 17 cm, was held fixed on the inclined plane. Hair gel was poured onto the plane far upstream from the obstacle, with sufficient flux and quantity to ensure that the plane was entirely coated by fluid, and that the thickness exceeded the rigid shape. Excess fluid was allowed to drain off the end of the plane. The pouring was then stopped, and the free surface approached its arrested shape as excess material drained around the obstruction.

To quantify the thickness of fluid, we projected an array of 76 parallel blue lines (e.g. figure 1 d) using a video projector (Epson EH-TW5200) mounted perpendicularly to the sloping surface. Similar methods have been employed elsewhere for free-surface flow of complex fluids (Hewitt & Balmforth Reference Hewitt and Balmforth2013), free-surface flow with an obstruction (Hinton, Hogg & Huppert Reference Hinton, Hogg and Huppert2020), and elastic membranes (Pihler-Puzović et al. Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015). However, we extended the method here to determine the full profile of the surface in two dimensions, $\hat {h}(\hat {x},\hat {y})$ . This technique was verified by employing it to measure the thickness of a 4 mm calibration plate, which exhibited an error of less than 2 %.

Hair gel was used as a representative viscoplastic material, owing to its established rheological properties (Taylor-West & Hogg Reference Taylor-West and Hogg2024), and the fact that it does not separate over the time scales required for these experiments. We used one brand of hair gel throughout the study (Woolworths Essentials ‘Firm Hold’). A small quantity of white acrylic paint was added (approximately 2 % v/v) to introduce opacity, resulting in strong diffuse scattering of light from the surface of the hair gel. Although this addition modified the yield stress, it preserved the viscoplastic nature of the material (this was verified using a cup and bob test with the Anton Paar MCR 702 rheometer; see § 5.2). The difference in the yield stress across the three experiments is a result of small variations in the proportion of paint.

The hair gel had density approximately $1000\,\mathrm{kg\,m^{-3}}$ , and the surface tension is approximately given by that of water, $\sigma = 0.07\, \mathrm{N\,m^{-1}}$ (Taylor-West & Hogg Reference Taylor-West and Hogg2024), which furnishes a Bond number

(5.1) \begin{equation} {{Bo} = \frac {\rho g L^2 }{\sigma }=1400}, \end{equation}

suggesting that surface tension effects are negligible in the experiments. The Reynolds number is given by

(5.2) \begin{equation} {{Re} = \frac {\rho U L}{\mu }=\frac {\rho L^2}{\mu \left ( \dfrac {3 \mu L}{\rho g \hat {H}_\infty ^2 \sin \beta }\right )}=0.015,} \end{equation}

where we have used the time scale from (2.1) and viscosity scale from (2.2), and the consistency is determined with a rheometer; see § 5.2. The role of inertia is negligible.

Figure 12. (a) The fluid thickness at $\hat {t}=60$ minutes from experiment 1 (square obstruction located in $(\hat {x},\hat {y}) \in [-100 \text{ mm},100\text{ mm}]^2$ ), which is compared with the lubrication prediction for the arrested state in (b). (c,d) The fluid thickness for experiment 2 (circular obstruction) and the arrested state prediction, respectively. Parameter values are given in table 1. White in (c) denotes a region with no data due to the circle obstructing the camera’s view.

In previous experimental studies that use hair gels (or similar carbomer fluids), the elastic storage modulus $G'$ at low strain is in the range $200 {-} 800 \text{ Pa}$ (Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; van der Kolk, Tieman & Jalaal Reference van der Kolk, Tieman and Jalaal2023; Taylor-West & Hogg Reference Taylor-West and Hogg2024). For our experiments, this furnishes an elastic time scale $\lambda =(k/G')^{1/n}$ in the range $ 10^{-5}{-}10^{-3} \text{ s}$ (the consistency is determined with a rheometer; see § 5.2). Two dimensionless groups quantify the relative importance of elasticity: the Weissenberg number and the Deborah number. The former is $\textit{Wi}=2 \lambda \dot {\gamma }$ , which vanishes as the slump becomes arrested, and is always below $10^{-3}$ during the motion. The Deborah number is the ratio of the elastic time scale to the flow time scale (which is of order hours); ${De}=\lambda /T_{\textit{flo}w}$ , so ${De}=10^{-8}{-}10^{-6}$ . These values suggest that elastic effects are negligible in our experiments.

Images were recorded without and with the fluid present, using a DSLR camera (Canon 5D Mark III) equipped with a telephoto lens (Canon 70–200 mm f/4 L IS USM) at focal length 89 mm. An aperture f/22 was chosen with a large depth of field to ensure that the full fluid surface was in focus. For each of the photographs, the projected lines were identified using cubic spline fitting of points identified as exceeding two standard deviations above the mean pixel intensity. The dimensional thickness of the material was then determined using a spatially dependent pixel calibration, which considers the variable distance to the camera, and accounts for the angle of the camera relative to the inclined plane. The 76 one-dimensional slices were interpolated to determine the dimensional thickness of the entire two-dimensional fluid upstream of the obstacle (e.g. see figures 12 a,c).

5.1. Experiments 1 and 2

In experiments 1 and 2, after stopping the upstream supply of fluid, we waited 60 minutes for fluid to drain and the slump to converge towards the arrested state. Experiment 1 was conducted with the square obstruction, and experiment 2 with the cylindrical obstruction; results for the fluid thickness at $\hat {t}=60$ minutes are shown in figure 12.

Figure 13. Fluid thickness along the centreline measured in (a) experiment 1 and (b) experiment 2 after 60 minutes, and the final arrested state predicted by (3.10a )–(3.15).

Table 1. Parameters for the three experiments conducted. For experiments 1 and 2, the upstream height $\hat {H}_\infty$ was measured at $\hat {x}=-370$ mm, with $b$ calculated using (1.6), and the yield stress using (1.2). For experiment 3, the yield stress was measured using a rheometer, and $b$ and $\hat {H}_\infty$ were calculated using (1.6) and (1.2), respectively.

The measured shape of the free surface along the centreline is shown in figure 13 along with the predictions of (3.10a ), (3.10b ) and (3.15). Parameters for these two experiments are reported in table 1. The value of $b$ is obtained by measuring the thickness of the layer far upstream and then using (1.6); this does not require direct knowledge of the yield stress. The yield stress is inferred from the height in the far field (at $\hat {x}=-370$ mm) using (1.2).

Figure 13 indicates that the free surface along the centreline after 60 minutes is substantially larger than the rigid shape predicted by our lubrication model. The measured free surface height is also higher than the predicted rigid shape away from the centreline, as seen in figure 12. This is further seen in the maximum height and excess volume of the fluid slump for experiment 1 (plotted as a diamond in figure 8) exceeding the predictions of the theory. This volume prediction is calculated for $-300\;\text{mm}\leqslant \hat {x}\leqslant -100\;\text{mm}$ to avoid edge effects. We also note that the measured shape does not exhibit a sharp cusp at the centreline ( $\hat {y}=0$ ).

The transient evolution of nine points along the centreline with $\hat {t} \in [0,60]$ minutes for experiment 1 is shown in figure 14. This indicates that the convergence to the rigid shape is algebraic in time, as suggested by our theoretical results, and that the decay is very slow (cf. Ancey & Cochard Reference Ancey and Cochard2009). Results from the Anton Paar MCR 702 rheometer found that the hair gel mixture’s rheology was accurately fitted by a Herschel–Bulkley model with $n=0.3$ , which suggests convergence with $\hat {t}^{-0.3}$ . Our results suggest that in the first 60 minutes, the decay rate is even slower than this, especially near the barrier. The slower time variation could be due to a residual slip velocity, or because the late-time behaviour has not yet been reached at $\hat {t}=1$ hour. These findings motivated an experiment with a much longer time scale.

Figure 14. Transient evolution of the measured fluid thickness at nine points along the centreline from experiment 1. The pouring of gel upstream stopped at $\hat {t}=0$ , but it takes approximately 100 s for this to substantially alter the fluid thickness in the vicinity of the barrier. The time evolution follows power-law behaviour; note the log-log scale.

Figure 15. Steady-state flow curve for the hair gel and paint mixture. Data shown for up- and down-stepped tests, and the best-fit Herschel–Bulkley constitutive law (yellow line).

Figure 16. (a) The fluid thickness at $\hat {t}=17$ hours from experiment 3. (b) Theoretical prediction for the static shape with no fitting parameters.

5.2. Experiment 3

We attribute the excess material in experiments 1 and 2 to the slump not reaching its final stationary state within the time of the experiment. In experiment 3, the fluid was allowed to drain and approach the arrested state for a significantly longer time period of 17 hours.

The constitutive relation for the hair gel–paint mixture from this experiment was measured using a rheometer, which gave $\tau _Y=57.05$ Pa, consistency $k=27.57$ Pa $\,{\rm s}^{0.3}$ , and Herschel–Bulkley power-law exponent $n=0.3$ ; see figure 15. Apart from at low strain rates, there is good agreement between the up- and down-stepped tests. There was a thirty-second rest period between the up and down data, suggesting that thixotropic effects are negligible.

The measured height of the free surface after 17 hours is shown in figure 16, with transects shown in figure 17. This is compared with our analytical predictions from (3.10a ), (3.10b ) and (3.15). For this experiment, $b=1.81$ (see table 1), which is obtained directly from the Anton Paar MCR 702 rheometer; there are no fitting parameters here (unlike in experiments 1 and 2).

The comparison in figures 16 and 17 shows reasonable agreement, with the theory under-predicting the fluid thickness closer to the obstruction, and over-predicting the thickness further away from the obstruction by up to 11 %. The maximum height and excess volume of the fluid slump are shown in figure 8 by a black $\times$ , and are in excellent agreement with the lubrication model prediction. Again, the measured slump does not exhibit sharp cusps at the centreline, in contrast with the theoretical prediction. This is consistent with past experiments (see e.g. Balmforth et al. Reference Balmforth, Craster, Rust and Sassi2006), and could be due to the finite time of the experiments, since the kink is only realised asymptotically as $\hat {t} \to \infty$ , surface wrinkling obscuring this phenomenon, or perhaps surface tension effects localised to the seam at $\hat {y}=0$ .

Figure 17. The measured height of the free surface after 17 hours (solid blue line), and height predicted by the theory (dotted black line) for experiment 3 along (a) the centreline, and along the transects (b) $\hat {x}=-150$ mm and (c) $\hat {x}=-200$ mm. (d) The relative height of the yield surface after 17 hours, as predicted by (1.5), where white regions denote missing data due to the free surface being obscured by the fluid.

We briefly comment on the expected approach of the fluid to its stationary state. The time scale defined in (2.1) is approximately $10$ seconds for experiment 3. Thus 17 hours corresponds to a dimensionless time $t\approx 6\times 10^3$ ; compare to figure 4(d). The algebraic decay to the stationary state dictated by (2.5) and our numerical results predict that the fluid thickness after 17 hours is approximately 8 % thicker than the final stationary state, in agreement with the observed discrepancy. To improve on this, an even longer experiment (say over a week) would be required. However, this is impractical as other effects such as evaporation would begin to become significant over such time scales. Indeed, dam-break experiments using Carbopol have found that evaporation becomes significant after time scales of 10 hours (Cochard Reference Cochard2007; Dubash et al. Reference Dubash, Balmforth, Slim and Cochard2009).

We note that the lubrication parameter for experiment 3, being the ratio of the far-field height to the half-width of the obstruction, is $\epsilon =\hat {H}_\infty /L=0.34$ . This indicates significant departure from the lubrication approximation, with non-lubrication effects contributing errors of order $\epsilon ^2$ , being $\approx 12\,\%$ . We also note that the experimental behaviour for $\hat {x}\lt -300$ mm is attributed to edge effects from the pouring of the fluid onto the board.

Figure 17(d) shows the height of the yield surface, $\hat {Y}$ after 17 hours. This was obtained from the measured height of the free surface by smoothing $\hat {h}$ in the $\hat {x}$ direction using a smoothing spline, and the height of the yield surface, $\hat {Y}$ , is then inferred from (1.5). Figure 17(d) suggests that the height of the yield surface does not exceed 20 % of the free surface height in the bulk of the domain. However, near the upstream obstruction corner (at $\hat {y}=\hat {x}=-100$ mm), the inferred height of the yield surface is non-negligible, which indicates that the lubrication approximation may not be accurate in this region. In particular, the requirement that shear stresses far exceed deviatoric normal stresses ( $|\tau _{zx}|,|\tau _{zy}| \gg |\tau _{xx}|, |\tau _{zz}|$ ), underpinning the lubrication approximation, likely does not hold near the corners, where the stress balance requires a full three-dimensional treatment. A similar observation has been made for viscoplastic material accumulating near the wall of a moving squeegee (Ribinskas et al. Reference Ribinskas, Ball, Penney and Neufeld2024; Taylor-West & Hogg Reference Taylor-West and Hogg2024).

6. Conclusion

We have analysed the shape of the rigid viscoplastic slump supported by a barrier on an inclined plane, as well as the transient convergence to this shape for a Herschel–Bulkley fluid. An exact expression for the rigid shape was found, which depends on a single dimensionless parameter that measures the relative magnitude of the yield stress. The arrested slump shape is independent of the viscous rheology and the nature of the upstream source of fluid.

The reported theory for the stationary slump shape can be extended to obstructions of different cross-section, provided that the length scale of variation in the barrier shape is larger than the flow thickness so that lubrication theory applies. The case of a circle and a wedge are straightforward extensions, discussed in Appendix D.

Laboratory experiments of the slumping of hair gel were reported. It was found that the convergence to the final rigid state is very slow, in accordance with the expectations of the model. Upon leaving the slump for 17 hours, we found that its shape can be approximated by our analytical solution. Our results demonstrate that the late-time shape of a supported slump could be used to infer the yield stress of the material provided that the static state has been obtained. This could provide a tool for estimating the yield stress of various hazardous environmental flows by observing the remnant material left upstream of obstructions.

Figure 18. (a) Maximum thickness of the slump and (b) magnitude of the force on the half of the barrier in $y\geqslant 0$ for the lava flow example discussed in § 6, as a function of barrier angle to the downslope direction, $\alpha$ .

The hydrostatic force on the barrier is given by (4.7) for a square barrier, and by (D7) for a barrier aligned at angle $\alpha$ to the downslope direction (for the half of the barrier in $y\geqslant 0$ ). We apply these results to an example lava flow barrier, which is motivated by those used during recent eruptions in Iceland (Sigtryggsdottir et al. Reference Sigtryggsdóttir, Hrafnsdóttir, Steingrímsson and Gudmundsson2025). The example barrier has length scale $L=100$ m, the underlying slope is $2^\circ$ , and the lava has yield stress $\tau _Y=2000$ Pa and density $\rho =3000$ kg m $^{-3}$ . These values furnish $b=0.56$ and upstream thickness $1.95$ m. The maximum flow thickness at the upstream wall and the magnitude of the force, $|(\hat {F}_x,\hat {F}_y)|$ , on the half of the barrier in $y\geqslant 0$ are shown in figure 18 as a function of barrier angle $\alpha$ , where $\alpha =\unicode{x03C0} /2$ corresponds to the square barrier. Although a barrier of the same extent that is closer to perpendicular to the oncoming flow protects a larger area downstream, it experiences a much larger force owing to the greater accumulation of material.

Acknowledgements

The authors are grateful to D.R. Hewitt and J.J. Taylor-West for valuable discussions. We would also like to thank R. Cavalida for carrying out the rheometer tests and providing the associated data concerning the hair gel–paint mixture. Three anonymous referees made insightful suggestions, particularly concerning the interpretation of the experimental results. E.M.H. acknowledges support from the Australian Research Council through a Discovery Early Career Researcher Award (DECRA:DE240100755).

Declaration of interests

The authors report no conflict of interest.

Appendix A. Numerical method

Equation (2.3) is numerically integrated using the method of lines. The spatial domain is discretised using an evenly spaced mesh, and derivatives are approximated using a first-order finite-volume method. The height at $x=-6$ is fixed to $H=1$ .

The resulting set of coupled, nonlinear ordinary differential equations for the height of the free surface at each grid point is integrated in time using a forward-Euler scheme, until the height of the yield surface is less than $10^{-4}$ over the whole spatial domain, i.e. $\max(Y(x,y))\lt 10^{-4}$ . Typically, this required integration over the time domain $t\in [0,10^4]$ . The number of points in the spatial mesh is doubled until the final stationary height of the free surface changes by less than 1 % between successive meshes. Use of higher-order time integration schemes (we used MATLAB’s ode45 solver, implementing fourth- and fifth-order Runge–Kutta methods) gives a final height of the free surface within 0.1 % of that calculated with the forward-Euler scheme.

Appendix B. Integrating Charpit’s equations

Charpit’s equations (3.8a )–(3.8e ) are solved following Balmforth et al. (Reference Balmforth, Craster and Sassi2002). Combining (3.8c ) and (3.8d ), we obtain (Balmforth et al. Reference Balmforth, Craster and Sassi2002)

(B1) \begin{align} p=a q, \end{align}

where $a=a(s)$ is a constant along each characteristic. We eliminate $p$ and $q$ by substituting (B1) into (3.2), giving

(B2) \begin{align} p & = \frac {a^2 H\pm a \sqrt {1+a^2-H^2}}{(1+a^2)H b} , \quad q = \frac {a H \pm \sqrt {1+a^2-H^2}}{(1+a^2)H b} . \end{align}

Substituting into (3.8c ) and (3.8d ), and dividing each by (3.8e ) furnishes

(B3a) \begin{align} \frac {{\rm d}y}{{\rm d}H} &= \pm \frac {{\textit{bH}}}{\sqrt {1+a^2-H^2}} , \end{align}
(B3b) \begin{align} \frac {{\rm d}x}{{\rm d}H} &= \pm \left (\frac {bH}{aH\mp \sqrt {1+a^2-H^2}}\right ) \left (\frac {H}{\sqrt {1+a^2-H^2}}\pm a\right )\!. \end{align}

Integrating (B3a ) and (B3b ) gives

(B4) \begin{align} y-y_0 & = \mp b\big[\sqrt {1+a^2-H^2}-\sqrt {1+a^2-H_0^2}\big], \end{align}
(B5) \begin{align} x-x_0 & = b \left [H-H_0(s)\right ] -b \log [\varDelta (H)/\varDelta (H_0)], \end{align}

with

(B6) \begin{align} \varDelta (H) & = \frac {a\pm \sqrt {1+a^2-H^2}}{H-1}, \end{align}

and $\varDelta (1)=1$ . This analysis applies to any obstruction shape; see also Appendix D. Substituting the initial data (subscript ‘0’) for each set of characteristics gives the reported solutions: (3.10a ) and (3.10b ) for the characteristics emanating from the obstruction face, and (3.15) for those emanating from the obstruction upstream corners.

We note that for characteristics emanating from the face of a square barrier, at $\tau =0$ , we obtain

(B7) \begin{equation} p_0=\frac {1}{b}, \quad q_0 = \pm \frac {1}{{\textit{bH}}_0}, \quad a= \frac {p_0}{q_0} = \pm H_0. \end{equation}

Appendix C. Small and large $\textbf{b}$ regimes

The slump shape and its volume can be derived analytically in the regimes $b\ll 1$ and $b \gg 1$ .

C.1. Small $b$ regime

For $b\ll 1$ , the rescaled free-surface height on the obstruction, (3.6), is

(C1) \begin{align} H_0(s) = \sqrt {\frac {2}{b}}\left (1-|s|\right ) + O(\sqrt {b}) . \end{align}

Combining this with the no-flux boundary condition (3.3) gives the scale for the extent of slump as $x+1\sim \sqrt {b}$ . This motivates defining the rescaled $x$ coordinate

(C2) \begin{align} X=\frac {x+1}{\sqrt {b}}, \end{align}

and expanding the slump height in powers of $\sqrt {b}$ :

(C3) \begin{align} H= \frac {1}{\sqrt {b}}\left (\bar {H}_0+\sqrt {b} \bar {H}_1 + O(b)\right )\! . \end{align}

Substituting (C2) and (C3) into (3.2), and expanding for small $b$ , we obtain the governing equations for $\bar {H}_0$ and $\bar {H}_1$ :

(C4) \begin{align} \left (1-\frac {\partial\kern-1pt\bar{H}_0}{\partial X}\right )^2=0 , \quad \frac {\partial\kern-1pt\bar{H}_1}{\partial X}=\pm \sqrt {\frac {1}{\bar {H}_0^2}-\left (\frac {\partial\kern-1pt\bar{H}_0}{\partial y}\right )^2}, \end{align}

subject to the boundary conditions

(C5) \begin{align} \bar {H}_0(-1,y) = \sqrt {\frac {2}{b}}\left (1-|y|\right ), \quad \bar {H}_1(-1,y) = 0 . \end{align}

Equation (C4) is solved via integration, giving

(C6a) \begin{align} \bar {H}_0&= \max \left(0,X+\sqrt {2(1-|y|)}\right)\! , \end{align}
(C6b) \begin{align} \bar {H}_1 &= \log \left (\frac {X+\sqrt {2(1-|y|)}}{\sqrt {-X\left (2+\dfrac {X}{\sqrt {2(1-|y|)}}\right )}-\sqrt {2(1-|y|)}}\right )\notag\\&\qquad\quad -\sqrt {-\frac {X}{\sqrt {2(1-|y|)}}\left (2+\frac {X}{\sqrt {2(1-|y|)}}\right )}. \end{align}

The asymptotic expansion breaks down in the vicinity of $X=-\sqrt {2(1-|y|)}$ , i.e. further from the barrier, due to the leading-order term vanishing. An internal boundary layer whose width is $O(b)$ arises at $X=-\sqrt {2(1-|y|)}$ , which describes the slump height in this region. Nonetheless, substituting (C6a ) and (C6b ) into (4.3) yields the excess volume for $b\ll 1$ , reported in (4.4a ). Analysis of the internal boundary layer is required to obtain the volume at $O(b)$ and at higher orders, but is not needed for the results presented here.

C.2. Large $b$ regime

For $b\gg 1$ and $x+1=O(1)$ , the free surface height is expanded in powers of $b^{-1}$ as

(C7) \begin{align} H = 1+b^{-1} \tilde {H}_1+b^{-2} \tilde {H}_2 + O(b^{-3}) . \end{align}

Substituting (C7) into (3.2), and expanding for large $b$ , we obtain the governing equation for $\tilde {H}_1$ ,

(C8) \begin{align} \left (1-\frac {\partial \tilde {H}_1}{\partial x}\right )^2 + \left (\frac {\partial \tilde {H}_1}{\partial y}\right )^2 = 1 , \end{align}

with the boundary condition

(C9) \begin{align} \tilde {H}_1(-1,y) = 1-|y| . \end{align}

This system is equivalent to the granular problem in Tregaskis et al. (Reference Tregaskis, Johnson, Cui and Gray2022). It is solved to obtain

(C10) \begin{align} \tilde {H}_1 = x+1+\sqrt {(x+1)^2+(1-|y|)^2}. \end{align}

However, (C10) does not describe the height of the free surface when $x+1=O(b)$ ; a boundary layer structure pertains far from the barrier. The outer region, where $x+1=O(b)$ , is analysed by defining the scaled coordinate

(C11) \begin{align} X = \frac {x+1}{b}. \end{align}

For large $b$ , the exact solution far from the obstruction is described by the corner characteristics (3.15). Substituting (C11) into (3.15) and expanding for large $b$ , we obtain the free-surface height solution

(C12) \begin{align} H=1+\frac {1}{b^2} \left [\frac {(1-|y|)^2}{2 \left (\exp (-X)-1\right )}\right ]\!. \end{align}

Combining (C10) and (C12), we construct a uniformly valid composite asymptotic solution for the free-surface height in the regime $b\gg 1$ :

(C13) \begin{align} H_{\textit{uni}} = 1 &+ \frac {1}{b}\left [x+1 +\sqrt {(x+1)^2+(1-|y|)^2}\right ] + \frac {1}{b^2} \left [\frac {(1-|y|)^2}{2\left(\exp \left(-\dfrac{(x+1)}{b} \right)-1\right)}\right ]\notag\\&+ \frac {1}{b}\left [\frac {1}{2}\frac {(1-|y|)^2}{x+1}\right ]\! . \end{align}

Substituting (C13) into the excess volume equation (4.3) yields the asymptotic result reported in (4.4b ).

Appendix D. Viscoplastic slump supported by general obstruction shapes

The method used to obtain the final rigid slump supported by a square obstruction (§ 3) can be generalised to an arbitrary obstruction shape. The shape of the upstream boundary of the obstruction is parametrised by $(x_0(s),y_0(s))$ , where $s$ gives an arc length of the boundary. Combining the governing equation and no-flux boundary condition, (3.2) and (3.3), the height of the free surface on the upstream wall of the obstruction face satisfies

(D1) \begin{align} H_0(s)\left \lvert b \frac {{\rm d} H_0}{{\rm d} s} - \frac {{\rm d}x_0}{{\rm d}s} \right \rvert = 1, \end{align}

with appropriate boundary conditions at the obstruction edges (typically $H_0=1$ ). For a cylindrical obstruction (see figure 19), the height of the free surface at the upstream face satisfies

(D2) \begin{align} H_0(\theta )\left \lvert b \frac {{\rm d} H_0}{{\rm d}\theta } + \sin \theta \right \rvert = 1 , \end{align}

where $\theta$ is the angle made with the positive $x$ axis. To obtain $H_0(\theta )$ requires numerical integration. The initial data can be written in terms of $H_0(\theta )$ as follows:

(D3a) \begin{align} x_0(\theta ) &= \cos \theta, \end{align}
(D3b) \begin{align} y_0(\theta ) &= \sin \theta, \end{align}
(D3c) \begin{align} p_0(\theta ) &= \frac {1}{b} \mp \frac { \tan \theta }{b\,H_0(\theta )\, \sqrt {1+\tan ^2 \theta }}, \end{align}
(D3d) \begin{align} q_0(\theta ) &= \frac {\pm 1}{b\,H_0(\theta )\, \sqrt {1+\tan ^2 \theta }}, \end{align}

for $\unicode{x03C0} /2 \lt \theta \lt 3 \unicode{x03C0} /2$ . Substituting this initial information into (B4) and (B5) then furnishes the stationary slump shape throughout the domain.

For a wedge cross-section, with the $y\gt 0$ part of the boundary having length unity and angle $\alpha$ to the $x$ axis, the shape on the upstream face is given implicitly by

(D4) \begin{equation} \left (\frac {s-1}{b} \right ) \cos ^2 \alpha = (H_0 -1) \cos \alpha + \log \left ( \frac {H_0 \cos \alpha -1}{\cos \alpha -1} \right )\!, \end{equation}

where $s=1$ is the end of the boundary where $H_0=1$ . The initial data can be written as follows:

(D5a) \begin{align} x_0(s) &= s \cos \alpha , \end{align}
(D5b) \begin{align} y_0(s) &= s \sin \alpha , \end{align}
(D5c) \begin{align} p_0(s) &= \frac {1}{b}-\frac {\cos \alpha }{b\,H_0(s)}, \end{align}
(D5d) \begin{align} q_0(s) &= -\frac {\sin \alpha }{b\, H_0(s)} . \end{align}

We can again use (B4) and (B5) with this initial information to obtain the stationary slump shape throughout the domain.

The dimensional force $(\hat {F}_x,\hat {F}_y)$ on the half of the barrier in $y\geqslant 0$ is given by

(D6) \begin{equation} (\hat {F}_x,\hat {F}_y)=(\sin \alpha , -\cos \alpha )\int _0^1 \frac {H_0(s)^2}{2} \, \mathrm{d} s\, \hat {H}_\infty ^2 L \rho g \cos \beta , \end{equation}

which furnishes the magnitude of the force,

(D7) \begin{equation} |\hat {F}|= \left [\frac {6}{b}-2(H_0(0)^3-1) \cos \alpha -3(H_0(0)^2-1) \right ]\frac {b}{12 \cos ^2 \alpha }\hat {H}_\infty ^2 L \rho g \cos \beta , \end{equation}

where

(D8) \begin{equation} H_0(0)=\frac {W\left [(\cos \alpha -1)\exp (\cos \alpha -1-b^{-1}\cos ^2 \alpha ) \right ] +1}{\cos \alpha }, \end{equation}

and $W$ is the Lambert W function.

Examples of the stationary slump solutions supported by obstructions with circular and wedge cross-sections are shown in figure 19. The excess volume $V(b)$ for a cylindrical obstruction is non-monontonic with $b$ (figure 20), so it is informative to analyse the regimes of small and large $b$ for a cylindrical obstruction.

Figure 19. Stationary slumps supported by (a) circular and (b) wedge-shaped obstructions for $b=1$ .

Figure 20. Dimensionless excess volume of static material supported by an obstruction with circular cross section. The $b\ll 1$ (D16) and $b\gg 1$ (D25) regimes are shown with dashed yellow and dotted red lines, respectively.

D.1. Cylindrical obstruction, small $b$ regime

For $b \ll 1$ , the solution on the boundary is given by (see (D2))

(D9) \begin{align} H_0(\theta ) = \frac {1}{|\sin \theta |} , \end{align}

for $\unicode{x03C0} /2 \lt \theta \lt 3 \unicode{x03C0} /2$ . The no-flux boundary condition then furnishes the scaling $r \sim b$ for the radial extent of the supported rigid material. This region therefore contributes ${O}(b)$ to the excess volume.

However, near the centreline where $\theta =\unicode{x03C0}$ , (D9) becomes singular, and for a balance in (D2), we scale

(D10) \begin{equation} \theta = \unicode{x03C0} - b^{1/3} \phi , \quad r=1+b^{2/3} \xi , \quad H=b^{-1/3} \tilde {H}, \end{equation}

so (D2) becomes

(D11) \begin{equation} \tilde {H}_0 \left (\frac {\mathrm{d} \tilde {H}_0}{\mathrm{d} \phi }- \phi \right ) = -1, \end{equation}

which has solution (cf. Hinton & Hogg Reference Hinton and Hogg2022; Hinton et al. Reference Hinton, Hewitt and Hogg2023)

(D12) \begin{equation} \phi = -2^{2/3} \frac {\text{Ai}'\big[(\phi ^2/2 -\tilde {H}_0)/2^{1/3}\big]}{\text{Ai}\big[(\phi ^2/2 -\tilde {H}_0)/2^{1/3}\big]}, \end{equation}

where $\text{Ai}$ is the Airy function, and we have used the condition that $\tilde {H}_0 \sim 1/\phi$ as $\phi \to \infty$ for matching with (D9). The maximum thickness occurs at $\phi =0$ and is given by $H=1.2836 b^{-1/3}$ . The leading-order solution in the region $\unicode{x03C0} -\theta =\mathcal{O}(b^{1/3})$ is obtained from (3.2) as

(D13) \begin{equation} \tilde {H} = \tilde {H}_0(\phi ) - \xi . \end{equation}

This region provides the dominant contribution to the excess volume:

(D14) \begin{align} V=\iint (H-1) \, \mathrm{d} x\,\mathrm{d} y. \end{align}

Substituting for $H$ , we find that to leading order,

(D15) \begin{align} V&= 2b^{2/3} \int _0^{\infty } \int _0^{\tilde {H}_0} (\tilde {H}_0- \xi ) \, \mathrm{d} \xi\, \mathrm{d} \phi \end{align}
(D16) \begin{align} &= 2b^{2/3} \int _0^{\infty } \frac {\tilde {H}_0^2}{2}\, \mathrm{d} \phi \nonumber \\&= 1.78 b^{2/3}. \end{align}

This result is plotted as a yellow dashed line in figure 20, where it is compared to the volume obtained from integrating Charpit’s equations (blue line). Although figure 20 indicates that the dimensionless excess volume is a non-monotonic function of $b$ , the dimensional excess volume is still an increasing function of the yield stress (as was the case for a square obstruction). The scaling for the excess volume is different to that for a square obstruction owing to the curvature at the upstream wall; similar variation between cross-sections was observed for free-surface viscoplastic flows past obstructions (see Hinton et al. Reference Hinton, Hewitt and Hogg2023).

D.2. Cylindrical obstruction, large $b$ regime

We restrict the analysis to $\unicode{x03C0} /2 \lt \theta \lt \unicode{x03C0}$ , which is sufficient due to the symmetry of the problem. For $b \gg 1$ , we write

(D17) \begin{align} H= 1 +b^{-1} \tilde {H}_1. \end{align}

On the upstream face, $r=1$ , (D2) becomes

(D18) \begin{align} \frac {\mathrm{d} \tilde {H}_1}{\mathrm{d} \theta } + \sin \theta = 1, \end{align}

with solution

(D19) \begin{align} \tilde {H}_1(r=1) = \cos \theta + \theta - \frac {\unicode{x03C0} }{2}. \end{align}

In $r\geqslant 1$ , $\tilde {H}_1$ satisfies the governing equation (3.2), which becomes

(D20) \begin{align} \left (\cos \theta - \frac {\partial \tilde {H}_1}{\partial r} \right )^2 + \left (\sin \theta + \frac {1}{r} \frac {\partial \tilde {H}_1}{\partial \theta } \right )^2=1, \end{align}

with solution

(D21) \begin{align} \tilde {H}_1 = r \cos \theta + \theta - \frac {\unicode{x03C0} }{2}+\sqrt {r^2 -1}-\tan ^{-1}\left ( \sqrt {r^2 -1} \right )\!. \end{align}

As for the square obstruction, a different asymptotic expansion is required distant from the obstruction. For this scenario, the expansion is required when $r=\mathcal{O}(b)$ and $\unicode{x03C0} - \theta = \mathcal{O}(1/b)$ , and then $H=1+\mathcal{O}(b^{-2})$ . On substituting and integrating, we find

(D22) \begin{align} H=1+b^{-2} \frac {(1+R \phi )^2}{2({\rm e}^R-1)}, \end{align}

with $R=r/b$ and $\theta = \unicode{x03C0} + \phi /b$ .

Combining the expansions in the two regions, we construct a uniformly valid composite asymptotic solution for the height of the free surface:

(D23) \begin{align} H=1&{}+b^{-1}\left [r \cos \theta + \theta - \frac {\unicode{x03C0} }{2}+\sqrt {r^2 -1}-\tan ^{-1}\left ( \sqrt {r^2 -1} \right )\right ] \end{align}
(D24) \begin{align} &{}+b^{-2} \frac {[1+r(\theta -\unicode{x03C0} )]^2}{2({\rm e}^{r/b}-1)} -b^{-1} \frac {[1+r(\theta -\unicode{x03C0} )]^2}{2r}. \end{align}

The excess volume is calculated as

(D25) \begin{align} V=2 \int _{\unicode{x03C0} /2}^{\unicode{x03C0} } \int _{1}^{1/\sin \theta } (H-1) r \, \mathrm{d} r\, \mathrm{d} \theta = \frac {0.126+\log b}{3b}. \end{align}

This result is plotted as a red dotted line in figure 20, where it is compared to the volume obtained from integrating Charpit’s equations.

References

Ancey, C. & Cochard, S. 2009 The dam-break problem for Herschel–Bulkley viscoplastic fluids down steep flumes. J. Non-Newtonian Fluid Mech. 158 (1–3), 1835.10.1016/j.jnnfm.2008.08.008CrossRefGoogle Scholar
Balmforth, N.J. & Craster, R.V. 1999 A consistent thin-layer theory for Bingham plastics. J. Non-Newtonian Fluid Mech. 84 (1), 6581.10.1016/S0377-0257(98)00133-5CrossRefGoogle Scholar
Balmforth, N.J., Craster, R.V., Perona, P., Rust, A.C. & Sassi, R. 2007 Viscoplastic dam breaks and the Bostwick consistometer. J. Non-Newtonian Fluid Mech. 142 (1–3), 6378.10.1016/j.jnnfm.2006.06.005CrossRefGoogle Scholar
Balmforth, N.J., Craster, R.V., Rust, A.C. & Sassi, R. 2006 Viscoplastic flow over an inclined surface. J. Non-Newtonian Fluid Mech. 139 (1), 103127.10.1016/j.jnnfm.2006.07.010CrossRefGoogle Scholar
Balmforth, N.J., Craster, R.V. & Sassi, R. 2002 Shallow viscoplastic flow on an inclined plane. J. Fluid Mech. 470, 129.10.1017/S0022112002001660CrossRefGoogle Scholar
Barberi, F., Brondi, F., Carapezza, M.L., Cavarra, L. & Murgia, C. 2003 Earthen barriers to control lava flows in the 2001 eruption of Mt. Etna. J. Volcanol. Geotherm. Res. 123 (1–2), 231243.10.1016/S0377-0273(03)00038-6CrossRefGoogle Scholar
Bjarnason, E. & Bubola, E. 2024 Iceland faces new chapter of seismic activity as lava menaces town. The New York Times (Digital Edition), January 15.Google Scholar
Chevrel, M.O., Harris, A., Ajas, A., Biren, J., Gurioli, L. & Calabrò, L. 2019 Investigating physical and thermal interactions between lava and trees: the case of Kılauea’s July 1974 flow. Bull. Volcanol. 81, 119.10.1007/s00445-018-1263-8CrossRefGoogle Scholar
Cochard, S. 2007 Measurements of time-dependent free-surface viscoplastic flows down steep slopes. PhD thesis, EPFL, Lausanne.Google Scholar
Cui, X. & Gray, J.M.N.T. 2013 Gravity-driven granular free-surface flow around a circular cylinder. J. Fluid Mech. 720, 314337.10.1017/jfm.2013.42CrossRefGoogle Scholar
Dietterich, H.R., Cashman, K.V., Rust, A.C. & Lev, E. 2015 Diverting lava flows in the lab. Nat. Geosci. 8 (7), 494496.10.1038/ngeo2470CrossRefGoogle Scholar
Dubash, N., Balmforth, N.J., Slim, A.C. & Cochard, S. 2009 What is the final shape of a viscoplastic slump? J. Non-Newtonian Fluid Mech. 158 (1–3), 91100.10.1016/j.jnnfm.2008.08.004CrossRefGoogle Scholar
Griffiths, R.W. 2000 The dynamics of lava flows. Annu. Rev. Fluid Mech. 32 (1), 477518.10.1146/annurev.fluid.32.1.477CrossRefGoogle Scholar
Hewitt, D.R. & Balmforth, N.J. 2013 Thixotropic gravity currents. J. Fluid Mech. 727, 5682.10.1017/jfm.2013.235CrossRefGoogle Scholar
Hinton, E.M. 2022 Inferring rheology from free-surface observations. J. Fluid Mech. 937, R4.10.1017/jfm.2022.157CrossRefGoogle Scholar
Hinton, E.M., Hewitt, D.R. & Hogg, A.J. 2023 Obstructed free-surface viscoplastic flow on an inclined plane. J. Fluid Mech. 964, A35.10.1017/jfm.2023.389CrossRefGoogle Scholar
Hinton, E.M. & Hogg, A.J. 2022 Flow of a yield-stress fluid past a topographical feature. J. Non-Newtonian Fluid Mech. 299, 104696.10.1016/j.jnnfm.2021.104696CrossRefGoogle Scholar
Hinton, E.M., Hogg, A.J. & Huppert, H.E. 2020 Viscous free-surface flows past cylinders. Phys. Rev. Fluids 5, 084101.10.1103/PhysRevFluids.5.084101CrossRefGoogle Scholar
Hogg, A.J. & Matson, G.P. 2009 Slumps of viscoplastic fluids on slopes. J. Non-Newtonian Fluid Mech. 158 (1), 101112.10.1016/j.jnnfm.2008.07.003CrossRefGoogle Scholar
Ionescu, I.R. 2013 Viscoplastic shallow flow equations with topography. J. Non-Newtonian Fluid Mech. 193, 116128.10.1016/j.jnnfm.2012.09.009CrossRefGoogle Scholar
van der Kolk, J., Tieman, D. & Jalaal, M. 2023 Viscoplastic lines: printing a single filament of yield stress material on a surface. J. Fluid Mech. 958, A34.10.1017/jfm.2023.118CrossRefGoogle Scholar
Matson, G.P. & Hogg, A.J. 2007 Two-dimensional dam break flows of Herschel–Bulkley fluids: the approach to the arrested state. J. Non-Newtonian Fluid Mech. 142 (1–3), 7994.10.1016/j.jnnfm.2006.05.003CrossRefGoogle Scholar
Moore, H.J. 1982 A geologic evaluation of proposed lava diversion barriers for the NOAA Mauna Loa Observatory. Tech. Rep. US Geological Survey, Mauna Loa volcano, Hawaii.Google Scholar
Ockendon, J.R., Howison, S., Lacey, A. & Movchan, A. 2003 Applied Partial Differential Equations. Oxford University Press.10.1093/oso/9780198527701.001.0001CrossRefGoogle Scholar
Pedersen, G.B.M., Pfeffer, M.A., Barsotti, S., Tarquini, S., de’Michieli Vitturi, M., Óladóttir, B.A. & Prastarson, R.H. 2023 Lava flow hazard modeling during the 2021 Fagradalsfjall eruption, Iceland: applications of MrLavaLoba. Nat. Hazards Earth Syst. Sci. 23 (9), 31473168.10.5194/nhess-23-3147-2023CrossRefGoogle Scholar
Pihler-Puzović, D., Juel, A., Peng, G.G., Lister, J.R. & Heil, M. 2015 Displacement flows under elastic membranes. Part 1. Experiments and direct numerical simulations. J. Fluid Mech. 784, 487511.10.1017/jfm.2015.590CrossRefGoogle Scholar
Ribinskas, E., Ball, T.V., Penney, C.E. & Neufeld, J.A. 2024 Scraping of a viscoplastic fluid to model mountain building. J. Fluid Mech. 998, A58.10.1017/jfm.2024.908CrossRefGoogle Scholar
Roussel, N. & Coussot, P. 2005Fifty-cent rheometer’ for yield stress measurements: from slump to spreading flow. J. Rheol. 49 (3), 705718.10.1122/1.1879041CrossRefGoogle Scholar
Sigtryggsdóttir, F.G., Hrafnsdóttir, H., Steingrímsson, J.H. & Gudmundsson, A. 2025 Experience in diverting and containing lava flow by barriers constructed from in situ material during the 2021 Geldingardalir volcanic eruption. Bull. Volcanol. 87, 28.10.1007/s00445-025-01806-3CrossRefGoogle Scholar
Taylor-West, J.J. & Hogg, A.J. 2024 Scraping of a thin layer of viscoplastic fluid. Phys. Rev. Fluids 9 (5), 053301.10.1103/PhysRevFluids.9.053301CrossRefGoogle Scholar
Thompson, R.L. & Soares, E.J. 2016 Viscoplastic dimensionless numbers. J. Non-Newtonian Fluid Mech. 238, 5764.10.1016/j.jnnfm.2016.05.001CrossRefGoogle Scholar
Tokpavi, D.L., Jay, P., Magnin, A. & Jossic, L. 2009 Experimental study of the very slow flow of a yield stress fluid around a circular cylinder. J. Non-Newtonian Fluid Mech. 164 (1–3), 3544.10.1016/j.jnnfm.2009.08.002CrossRefGoogle Scholar
Tregaskis, C., Johnson, C.G., Cui, X. & Gray, J.M.N.T. 2022 Subcritical and supercritical granular flow around an obstacle on a rough inclined plane. J. Fluid Mech. 933, A25.10.1017/jfm.2021.1074CrossRefGoogle Scholar
Figure 0

Figure 1. Schematics of viscoplastic flow around an obstruction. (a) Steady flow. (b) Shape of the stationary slump following termination of the upstream source. (c) Cross-section of the steady flow along the centreline. (d) Photo of laboratory experiment 1; see § 5 (obstruction has extent 200 mm).

Figure 1

Figure 2. Steady flowing solutions for $b=1$ and: (a) Bingham fluid, $H_\infty =2.53$ (corresponding to $q_\infty =7.1$); (b) Bingham fluid, $H_\infty =1.68$ ($q_\infty =1$); (c) Herschel–Bulkley fluid, $n=0.5$, $H_\infty =1.30$ ($q_\infty =0.03$); see (1.3). When the flux source is switched off, all of these converge to the same final stationary slump, shown in (d). Note the different scale bars used across the panels.

Figure 2

Figure 3. Late-time slumping of a Bingham fluid. (a,b) The height of the free surface (dashed black lines) at the upstream wall of the obstruction ($x=-1$) at times $t\in \{\textrm{e}^1,\textrm{e}^2,\textrm{e}^3,\textrm{e}^4,\textrm{e}^5,\textrm{e}^6\}$ after the upstream flux (at $x=-6$) is terminated. (c,d) Corresponding results along the centreline ($y=0$). In (a) and (c), parameter values are $b=1,\ q_\infty =1$, whilst for (b) and (d), $b=1,\ q_{\infty }=1/2$. The final shape is shown with a blue line (from § 3).

Figure 3

Figure 4. Late-time decay of (a,b,c) a Bingham fluid to the final state, and (d,e,f) a Herschel–Bulkley fluid with $n=1/4$. (a,d) Transient evolution of the relative flow thickness at four points along the centreline ($y=0$); these show the late-time decay rate $\sim 1/t^n$. (b,e) The relative difference between the fluid thickness at $t=20$ and the final rigid thickness ($t \to \infty$). (c, f) The relative height of the yield surface at $t=20$. Parameter values are $b=1$ and $H_\infty =1.68$.

Figure 4

Figure 5. Streamlines of the flow (white lines) and magnitude of the flux (colour bar) for a Bingham fluid at $t=20$. Parameter values are $b=1$ and $H_\infty =1.68$.

Figure 5

Figure 6. (a) Characteristics for the thickness of the rigid shape supported by the obstruction for $b=1/2$. (b) Characteristic projections in the $(x,y)$ plane corresponding to (a). (c) Characteristic projections for $b=1$. (d) Characteristic projections for $b=2$. The solid black lines represent characteristics emanating from the upstream wall, whilst the red dashed lines represent those emanating from the corners. See also figure 7.

Figure 6

Figure 7. The exact solution for the stationary slump ((3.10a), (3.10b) and (3.15)) for (a) $b=1/4$, (b) $b=1/2$, (c) $b=1$, (d) $b=2$.

Figure 7

Figure 8. (a) The dimensionless maximum slump height as function of $b$. (b) The dimensional maximum height as a function of the yield stress $\tau _Y$ for $\rho g = 10^4$$\mathrm{kg\,m^{-2}\,s^{-2}}$, $L=0.1$$\mathrm{m}$ and $\beta =5 ^{\circ }$. (c) Dimensionless excess volume of the slump ((4.3), solid blue line), compared with the small and large $b$ asymptotic solutions (see (4.4a)–(4.4b)). (d) Dimensional volume as a function of $\tau _Y$ for the above-mentioned parameter values. The scaled maximum height and volume obtained in experiments 1 and 3 (see § 5) are shown in (a) and (c), where the error bars are smaller than the extent of the markers.

Figure 8

Figure 9. The streamlines of the numerically calculated late-time flow (solid blue lines) compared with the characteristic projections from Charpit’s method (dashed black lines), for $b=1$. The obstruction is shown in grey.

Figure 9

Figure 10. (a,c) Flowing and (b,d) stationary solutions for Bingham fluid ($n=1$) with $b=1/4,\ q_\infty =1$, for (a,b) off-centre vent source, (c,d) infinite line source.

Figure 10

Figure 11. Fluid thickness at $x=1$ (i.e. the downstream face of the obstruction), at various times after the source has been switched off: (a) infinite line source, (b) vent source, as described in figure 10. The dashed black lines indicate the edges of the obstruction’s downstream face.

Figure 11

Figure 12. (a) The fluid thickness at $\hat {t}=60$ minutes from experiment 1 (square obstruction located in $(\hat {x},\hat {y}) \in [-100 \text{ mm},100\text{ mm}]^2$), which is compared with the lubrication prediction for the arrested state in (b). (c,d) The fluid thickness for experiment 2 (circular obstruction) and the arrested state prediction, respectively. Parameter values are given in table 1. White in (c) denotes a region with no data due to the circle obstructing the camera’s view.

Figure 12

Figure 13. Fluid thickness along the centreline measured in (a) experiment 1 and (b) experiment 2 after 60 minutes, and the final arrested state predicted by (3.10a)–(3.15).

Figure 13

Table 1. Parameters for the three experiments conducted. For experiments 1 and 2, the upstream height $\hat {H}_\infty$ was measured at $\hat {x}=-370$ mm, with $b$ calculated using (1.6), and the yield stress using (1.2). For experiment 3, the yield stress was measured using a rheometer, and $b$ and $\hat {H}_\infty$ were calculated using (1.6) and (1.2), respectively.

Figure 14

Figure 14. Transient evolution of the measured fluid thickness at nine points along the centreline from experiment 1. The pouring of gel upstream stopped at $\hat {t}=0$, but it takes approximately 100 s for this to substantially alter the fluid thickness in the vicinity of the barrier. The time evolution follows power-law behaviour; note the log-log scale.

Figure 15

Figure 15. Steady-state flow curve for the hair gel and paint mixture. Data shown for up- and down-stepped tests, and the best-fit Herschel–Bulkley constitutive law (yellow line).

Figure 16

Figure 16. (a) The fluid thickness at $\hat {t}=17$ hours from experiment 3. (b) Theoretical prediction for the static shape with no fitting parameters.

Figure 17

Figure 17. The measured height of the free surface after 17 hours (solid blue line), and height predicted by the theory (dotted black line) for experiment 3 along (a) the centreline, and along the transects (b) $\hat {x}=-150$ mm and (c) $\hat {x}=-200$ mm. (d) The relative height of the yield surface after 17 hours, as predicted by (1.5), where white regions denote missing data due to the free surface being obscured by the fluid.

Figure 18

Figure 18. (a) Maximum thickness of the slump and (b) magnitude of the force on the half of the barrier in $y\geqslant 0$ for the lava flow example discussed in § 6, as a function of barrier angle to the downslope direction, $\alpha$.

Figure 19

Figure 19. Stationary slumps supported by (a) circular and (b) wedge-shaped obstructions for $b=1$.

Figure 20

Figure 20. Dimensionless excess volume of static material supported by an obstruction with circular cross section. The $b\ll 1$ (D16) and $b\gg 1$ (D25) regimes are shown with dashed yellow and dotted red lines, respectively.