2. Theoretical model

Figure 1.

Schematic of the slumping of a viscoplastic material atop an oscillating plate. (a) Initial shape, which is then oscillated (in (b)), and slumps towards its arrested shape (c). The height of the deposit has been exaggerated for visibility.

We investigate the gravity-driven flow of a two-dimensional deposit of Bingham fluid of area $A$ lying atop a solid horizontal plate located at $z=0$ , with free-surface profile $z=h(x,t)$ (see figure 1). The fluid is assumed to have completely slumped under gravity to its static profile (1.1) before the oscillatory plate motion commences. For $t\geqslant 0$ the basal plate oscillates with velocity $U_0 \sin (\omega t)$ in the $x$ -direction, where $\omega$ is the oscillation frequency. The deposit then occupies the region between its two contact points $x_{\kern-1pt f}^- (t) \leqslant x \leqslant x_{\kern-1pt f}^+ (t)$ . The fluid is incompressible and subject to Cauchy’s momentum equation and mass conservation, respectively,

(2.1a) \begin{align} \rho \left ( \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} \right ) &= -\boldsymbol{\nabla }p + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau } - \rho g \hat {\boldsymbol{z}}, \\[-12pt] \nonumber \end{align}

(2.1b) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0 , \\[10pt] \nonumber \end{align}

where $\boldsymbol{u}=(u, w)$ is the fluid velocity in the $x$ and $z$ directions, $\boldsymbol{\tau }=\tau _{\textit{ij}}$ is the deviatoric stress tensor and $p$ is the pressure. The fluid obeys the Bingham constitutive law

(2.2a)

(2.2b)

where $\mu$ is the shear viscosity of the yielded material, $\dot {\gamma }_{\textit{ij}}$ is the strain-rate tensor, and $|{\tau }|=\sqrt {\tau _{\textit{ij}} \tau _{\textit{ij}}/2}$ and $|{\dot {\gamma }}|=\sqrt { \dot {\gamma }_{\textit{ij}} \dot {\gamma }_{\textit{ij}}/2}$ are the second invariants of their respective tensors.

The dynamical regime of interest is that of relatively thin deposits, so characteristic height and width scales are identified from the gravitationally slumped initial profile (1.1)–(1.2)

(2.3) \begin{align} H=A^{1/3} \left (\frac {\tau _Y}{\rho g}\right )^{1/3}, \qquad L=A^{2/3} \left (\frac {\tau _Y}{\rho g}\right )^{-1/3}. \end{align}

The aspect ratio $\varepsilon$ defines the ratio of these disparate length scales

(2.4) \begin{align} \varepsilon =\left (\frac {\tau _Y}{\rho g \sqrt {A}}\right )^{2/3}, \end{align}

and henceforth $\varepsilon \ll 1$ is assumed. This condition, specified by $\rho g \sqrt {A} \gg \tau _Y$ , corresponds to a regime where the hydrostatic pressure gradients have sufficiently overcome the yield stress resulting in a relatively shallow free-surface profile under the action of gravity.

The initial conditions are (Liu & Mei Reference Liu and Mei1989)

(2.5) \begin{align} \begin{split} u(x,z,0)=0, \quad h(x,0)= \sqrt {\frac {2 \tau _Y}{\rho g} (x_{\kern-1pt f}{(0)} - |{x}|)}, \\ \tau _{xz}(x,z,0)= \tau _Y \text{sgn}(x)\left [1-\frac {z}{\sqrt {\dfrac{2 \tau _Y}{\rho g} (x_{\kern-1pt f}{(0)} - |{x}|)}}\right ]\!. \end{split} \end{align}

The boundary conditions are

(2.6) \begin{align} \begin{split} u(x,0,t)=U_0 \sin (\omega t), \quad w(x,0,t)=0, \quad h(x_{\kern-1pt f}^{\pm } (t),t)=0, \\ w(x,h,t)= \frac {\partial h}{\partial t} + u(x,h,t) \frac {\partial h}{\partial x}, \qquad \tau _{xz}(x,h,t)=0. \end{split} \end{align}

It is assumed that the deposit is large enough to ignore capillary effects. This regime is specified by $\sigma \ll \rho \omega U_0 L^3/H$ , where $\sigma$ is the surface tension. Additionally, the height of the deposit is assumed to be much larger than the slip length, and so the no-slip condition is applied at the plate. We also assume that the elastic relaxation time of the material is much more rapid than the shear rate: the Weissenberg number $\textit{Wi}=U_0\mu /\textit{GH} \ll 1$ , where $G$ is the elastic shear modulus of the material.

A balance in the continuity (2.1b ) reveals that the horizontal and vertical velocity scales are $U_0$ and $ U_0 H/L$ , respectively. The system is non-dimensionalised by scaling the $z$ -coordinate and free-surface profile with $H$ , the $x$ -coordinate with $L$ and time with the scale for oscillation $1/\omega$ . The pressure and stresses are non-dimensionalised according to their inertia-driven scales in (2.1a ), which are $\rho L \omega U_0$ and $\rho H \omega U_0$ , respectively. Altogether,

(2.7) \begin{align} \begin{split} \hat {u}=\frac {u}{U_0}, \quad \hat {w}=\frac {L w}{U_0 H}, \quad \hat {t}=\omega t, \quad \hat {h}=\frac {h}{H}, \quad \hat {x}=\frac {x}{L}, \quad \hat {z}=\frac {z}{H}, \\ \hat {p}=\frac {p}{\rho L \omega U_0}, \quad {\left \{\hat {\tau }_{xx},\hat {\tau }_{xz},\hat {\tau }_{zz}\right \}=\frac {1}{\rho H \omega U_0} \left \{{\tau }_{xx},{\tau _{xz}},{\tau }_{zz}\right \}}, \end{split} \end{align}

where the hatted symbols signify the corresponding dimensionless variables.

Henceforth, the hats are omitted and all variables are taken to be non-dimensional. The associated dimensionless groups are

(2.8) \begin{align} \delta =\frac {U_0}{L\omega }, \qquad B=\frac {\tau _Y}{\rho H \omega U_0}, \qquad \alpha =\frac {\rho H^2 \omega }{\mu }. \end{align}

The dimensionless oscillation amplitude $\delta$ defines the ratio of the oscillation amplitude $U_0 /\omega$ of the plate relative to the length of the initial deposit. The Bingham number $B$ defines the ratio of the yield stress to the internal shear stress scale, and the oscillatory Reynolds number $\alpha$ defines the time scale for momentum diffusion relative to the period of the oscillating plate. The quantity $\delta \alpha$ is the lubrication Reynolds number.

To simplify the analysis, the governing equations are expressed in coordinates co-moving with the plate

(2.9) \begin{align} \eta = x-\delta \left (1 - \cos t\right ),\qquad \mathscr{u}=u-\sin t. \end{align}

The deposit occupies the region $\eta _{\kern-1pt f}^{-} (t)\leqslant \eta \leqslant \eta _{\kern-1pt f}^{+} (t)$ , where $\eta _{\kern-1pt f}^{\pm } (t)=x_{\kern-1pt f}^{\pm } (t)-\delta (1 - \cos t )$ denotes the contact points in the $\eta$ frame. The dimensionless versions of the governing (2.1) in the moving frame are given by

(2.10a) \begin{align} \frac {\partial \mathscr{u}}{\partial t} + \cos {t} + \delta \mathscr{u} \frac {\partial \mathscr{u}}{\partial \eta } + \delta w \frac {\partial \mathscr{u}}{\partial z}&=-\frac {\partial p}{\partial \eta } + \frac {\partial \tau _{\eta z}}{\partial z}+{\varepsilon } \frac {\partial \tau _{\eta \eta }}{\partial \eta }, \\[-12pt] \nonumber \end{align}

(2.10b) \begin{align} \varepsilon ^2 \frac {\partial w}{\partial t} + \delta \varepsilon ^2 \mathscr{u} \frac {\partial w}{\partial \eta } + \delta \varepsilon ^2 w \frac {\partial w}{\partial z}&=-\frac {\partial p}{\partial z} + \varepsilon ^2 \frac {\partial \tau _{\eta z}}{\partial \eta }+{\varepsilon } \frac {\partial \tau _{z z}}{\partial z} - B, \\[-12pt] \nonumber \end{align}

(2.10c) \begin{align} \frac {\partial \mathscr{u}}{\partial \eta } + \frac {\partial w}{\partial z} &= 0. \\[10pt] \nonumber \end{align}

The velocity $\mathscr{u}$ now vanishes at the plate, with the oscillatory motion manifesting as a forcing term, $\cos t$ , on the left-hand side of (2.10a ). In this paper, the model will be interrogated to leading order in $\varepsilon$ . In this limit, the vertical momentum (2.10b ) identifies that the pressure is hydrostatic. Additionally, the constitutive law (2.2) reveals that the normal stresses are $\mathcal{O}{(\varepsilon )}$ smaller than the shear stress in yielded regions, so the yield criterion reduces to $|{\tau }|=|{\tau _{\eta z}}| \gt \tau _Y$ . Finally, the only stress contribution in the leading-order shallow system is from the shear stress $\tau _{\eta z}$ , which will henceforth be written as $\tau$ . The horizontal momentum (2.10a ) now reduces to

(2.11) \begin{align} \frac {\partial \mathscr{u}}{\partial t} + \cos t +\delta \mathscr{u} \frac {\partial \mathscr{u}}{\partial \eta } + \delta w \frac {\partial \mathscr{u}}{\partial z}=-B \frac {\partial h}{\partial \eta } + \frac {\partial \tau }{\partial z}, \end{align}

with the constitutive law

(2.12a)

(2.12b)

Depth integrating the continuity (2.10c ) across the deposit and applying the kinematic condition in (2.6) at the free surface furnishes

(2.13) \begin{align} \frac {\partial h}{\partial t}+\delta \frac {\partial }{\partial \eta } \int _0^h \mathscr{u} \, \text{d}z =0. \end{align}

The dimensionless initial and boundary conditions are

(2.14a) \begin{align} \mathscr{u}(\eta ,z,0) & =0, \ \ h(\eta ,0)= \sqrt {2 \left ( \eta _{\kern-1pt f}{(0)} - |{\eta }|\right )}, \nonumber \\ \tau (\eta ,z,0) & = B \, \text{sgn}(\eta )\left [1-\frac {z}{\sqrt {2 \left ( \eta _{\kern-1pt f}{(0)} - |{\eta }|\right )}}\right ] \!, \\[-28pt] \nonumber \end{align}

(2.14b) \begin{align} \mathscr{u}(\eta ,0,t)= 0, \quad h(\eta _{\kern-1pt f}^{\pm } (t),t)=0, \quad \tau (\eta ,h,t)=0, \\[0pt] \nonumber \end{align}

where the scaled half-width of the initial free-surface profile is given by ${\eta _{\kern-1pt f}{(0)}}= (9/32 )^{1/3}$ . Since the shear stress vanishes at the free surface, there must exist a rigid layer adjacent to the free surface in which $|{\tau }| \leqslant B$ . An equation for the uppermost yield surface $Y_1$ which bounds this rigid region from below may be derived by integrating the rigid momentum (2.11) over the region $Y_1 (\eta ,t) \lt z \lt h(\eta ,t)$ , giving (Sekimoto Reference Sekimoto1993)

(2.15) \begin{align} \frac {\partial \mathscr{u}}{\partial t}+\delta \mathscr{u} \frac {\partial \mathscr{u}}{\partial \eta }=\frac {\partial \mathscr{u}_R}{\partial t} +\delta \mathscr{u}_R \frac {\partial \mathscr{u}_R}{\partial \eta }=-\cos {t} - B \frac {\partial h}{\partial \eta } - \frac {B \, \text{sgn}{\left (\tau \right )}}{h-Y_1} \quad \text{at} \quad z=Y_1 (\eta ,t), \end{align}

where the subscript $R$ henceforth refers to the rigid region; that is, $\mathscr{u}_R (\eta ,t)=\mathscr{u}(\eta ,Y_1(\eta ,t),t)$ is the velocity of the pseudo-plugged flow, or equivalently, the velocity at the free surface. Equation (2.15) is valid throughout the uppermost rigid layer. A single yield surface emerges from the plate when the basal stress breaches the yield stress, initially propagating upwards and then eventually downwards through the deposit, colliding with the plate when the material rigidifies there. Multiple yield surfaces – and accompanying rigid regions – may coexist if $\alpha$ is large, and $B$ is sufficiently small (Hinton et al. Reference Hinton, Collis and Sader2023a ). In this regime, momentum diffusion is sufficiently slow relative to the plate oscillation, so an additional yield surface may emerge at the plate before the yield surfaces already present in the deposit have collided with the plate. The yield surfaces are denoted $z=Y_i (\eta ,t)$ for $i=1,2,\ldots , N$ , with $Y_1 (\eta ,t) \gt Y_2(\eta ,t) \gt \ldots \gt Y_N(\eta ,t)$ for a fixed $\eta$ . From (2.12), it follows that the strain rate $\dot {\gamma }=\partial \mathscr{u} / \partial z$ , must be continuous at each yield surface. The velocity must also be continuous everywhere, and requiring equality of $\mathscr{u}$ on each yield surface provides a condition coupling motion in the rigid and yielded regions, hence closing the initial-boundary-value problem (2.11–2.14).

3. Numerical results

In this section, we describe the numerical method and validate it against some of our asymptotic results. Numerical results are presented for a region of the $(B, \alpha , \delta )$ parameter space in which $\delta$ is small.

3.1. Numerical method

Equations (2.11–2.13) are solved numerically with respect to the initial and boundary conditions (2.14) for the free-surface profile $h(x,t)$ , shear stress $\tau (x,z,t)$ and strain rate $\dot {\gamma }=\partial \mathscr{u} / \partial z$ . The strategy relies on evolving $\dot {\gamma }$ in time, and calculating $\tau$ using the inverse of the constitutive law (2.12)

(3.1) \begin{align} \dot {\gamma }=\alpha \, \text{sgn}{(\tau )} \left (|{\tau }|-B\right )\varTheta {(|{\tau }|-B)}, \end{align}

where $\varTheta (\boldsymbol{\cdot })$ is the Heaviside step function. The numerical method solves the stress on a grid that is fully two-dimensional. The momentum (2.11) is differentiated with respect to $z$ to obtain

(3.2) \begin{align} \frac {\partial \dot {\gamma }}{\partial t}+\delta \mathscr{u}\frac {\partial \dot {\gamma }}{\partial \eta } + \delta w \frac {\partial \dot {\gamma }}{\partial z}=\frac {\partial ^2 \tau }{\partial z^2}, \end{align}

where the velocities are written in terms of the strain rate via

(3.3) \begin{align} \mathscr{u}=\int _0^z \dot {\gamma }(\eta ,z',t) \, \text{d} z', \qquad w=-\int _0^z \frac {\partial }{\partial \eta } \mathscr{u}(\eta ,z',t) \, \text{d} z'. \end{align}

By integrating the flux term in the conservation of mass (2.13) by parts, we obtain an equation for the evolution of the free surface in terms of $\dot {\gamma }$ , rather than $\mathscr{u}$

(3.4) \begin{align} \frac {\partial h}{\partial t}=-\delta \frac {\partial }{\partial \eta } \int _0^h (h-z) \dot {\gamma } \, \text{d}z. \end{align}

This formulation avoids the need to explicitly determine the yield surfaces (see also Balmforth et al. Reference Balmforth, Forterre and Pouliquen2009). The basal velocity boundary condition (2.14b ) can be expressed in terms of the shear stress by substituting the condition into the momentum (2.11), so all of the relevant boundary conditions are

(3.5) \begin{align} \frac {\partial \tau }{\partial z} \bigg {\rvert }_{z=0}=\cos t + \frac {\partial h}{\partial \eta }, \qquad \tau \rvert _{z=h}=0 . \end{align}

with the initial conditions given in (2.14a ). Full details of the numerical implementation can be found in Appendix A.

3.2. Physical features

Figure 2.

The free-surface profile at $t=1, 10, 100, 1000, 10\,000$ for twelve pairs of $(B,\alpha )$ , and $\delta =0.05$ . The black lines represent the numerical solution for $h$ , and the red dotted lines are the numerical integration of the low- $\alpha$ asymptotic (4.5). The initial free-surface profile $h{(\eta ,0)}$ is the blue dashed line and the arrested state $h_{\infty }{(\eta )}$ (whose solution is derived in § 5) is the solid blue line. As $\alpha$ increases beyond $\mathcal{O}(1)$ , the low- $\alpha$ free-surface profile deviates from the numerical solution at early times. At later times, even for $\alpha \gt \mathcal{O}{(1)}$ , the free surface is well described by the low- $\alpha$ solution because the inertia is less significant.

Results for the evolution of the free-surface profile are presented in figure 2, captured at $t=1$ , $10$ , $100$ , $1000$ and $10\,000$ for $(\alpha ,B)\in \{0.05,0.5,5,20\}\times \{0.3,0.8,1.5\}$ . The deposit slumps away from its initial shape (the dashed blue line) towards a symmetric arrested state (the solid blue line). Larger $\alpha$ values result in more rapid slumping, because momentum diffusion is suppressed and $\mathscr{u}=\mathcal{O}{(1)}$ .

Initially, the deformation of the free surface is confined to the centre of the deposit, driven by the discontinuity of the free-surface gradient at $\eta =0$ . Disturbances propagate laterally towards the contact points, reaching them when $t \sim 1/\delta \alpha$ . In the early-time regime, where $t \leqslant \mathcal{O}{(1/\delta \alpha )}$ , the contact points are only weakly advanced by these disturbances while the free surface reshapes behind them, approximately presenting as a ‘waiting time solution’ (Gratton, Minotti & Mahajan Reference Gratton, Minotti and Mahajan1999). This phenomenon can be described more specifically by analysing the deposit at each front, where the free-surface gradients dominate (2.11) and (2.15) and the leading-order contact point evolution is described by

(3.6) \begin{align} \frac {\text{d} \eta _{\kern-1pt f}^{\pm }}{\text{d} t} \sim \lim _{\eta \rightarrow \eta _{f}^{\pm }} \frac {B \alpha \delta }{h} \frac {\partial h}{\partial \eta } \left (\frac {h Y_{1}^2}{2} - \frac {Y_{1}^3}{6}\right ), \qquad Y_1 \sim \text{max}\left \{0, \, h-\frac {1}{|{{\partial h}/{\partial \eta }}|} \right \}. \end{align}

Thus, the free surface about the contact points must take the form $h \propto (\eta ^{\pm }_f(t)-|{\eta }|)^{1/3}$ to propagate at the same asymptotic order as the free surface evolution, which is a common result for creeping gravity currents; see Perazzo & Gratton (Reference Perazzo and Gratton2003), Matson & Hogg (Reference Matson and Hogg2007) and Appendix C. Evidently the front propagation is sub-leading order at early times, until the disturbances reach the contact points and the fronts steepen accordingly.

Prominent asymmetry is visible in the bulk of the free-surface profile during the first few oscillations when $\alpha$ is large. In addition, once the early-time regime is breached, additional asymmetry in the contact points is evident, driven by the ordering of each oscillation direction. This asymmetry becomes more pronounced with larger $\delta \alpha$ , as evident from the propagation rate (3.6). As an aside, we observe that if $\delta \alpha \gg 1$ , the initial ‘waiting time’ period is negligible and the contact points may advance during the first oscillation. A host of interesting nonlinear behaviour can occur in this regime. In particular, if $B$ is sufficiently large, the negative $\eta$ contact point may advance beyond the predicted arrested position, and the deposit appears to slump towards a separate, asymmetric arrested state. Analysis of this regime is beyond the scope of this paper, which we restrict to the case where $\delta \alpha$ is small.

In this case, a late-time regime is eventually reached wherein the free surface is close to its symmetric arrested state. Inertial forcing is almost entirely opposed by the yield stress, and little variation in the free surface occurs during a period of oscillation. Despite any pronounced early-time asymmetries, the deposit is only weakly asymmetric in this regime. Because inertial effects are relatively small, the behaviour is described by the dynamics in the low $\alpha$ case.

Figure 3.

Numerical solutions of the free-surface profile at $t=1, 10, 100, 1000$ for $B=0.8$ , $\alpha =5$ and $\delta =0.05$ , beginning from two different initial shapes with area of unity. The red curves indicate the behaviour with initial condition given in (2.14). The solutions for some elliptical initial conditions are given by the black curves: in (a) the ellipse is chosen to have the same initial width as in (2.14), while in (b) the ellipse is chosen to have the same maximum height as in (2.14). The dotted lines are the respective initial profiles.

Figure 4.

The yield surfaces plotted beneath the free-surface profile when $t=5 \pi /4$ for sixteen pairs of $(B,\alpha )$ , and $\delta =0.05$ . The solid black lines are the free surfaces, and the black dot-dash lines are the yield surfaces, obtained from the numerical solution. The red dotted lines are low- $\alpha$ results for the free surfaces and yield surfaces, obtained from numerical integration of (4.3)–(4.5). The material immediately beneath the free surface is always rigid. The vertical grey dashed lines are at $\eta =\eta _{\kern-1pt f}{(0)}/2$ , where the spatial snapshot is taken, shown in figure 5.

Figure 5.

The yield surfaces plotted beneath the free-surface profile for $\pi \leqslant t \leqslant 3\pi$ at $\eta =\eta _{\kern-1pt f}{(0)}/2$ for sixteen pairs of $(B,\alpha )$ , and $\delta =0.05$ . The solid black lines are the free surfaces, and the black dot-dash lines are the yield surfaces, obtained from the numerical solution. The red dotted lines are low- $\alpha$ results for the free surfaces and yield surfaces, obtained from numerical integration of (4.3–4.5). The material immediately beneath the free surface is always rigid material, and the yield surfaces partition this region from the yielded material. The vertical grey dashed lines are at $t=5\pi /4$ , where the temporal snapshot is taken, shown in figure 4.

The sensitivity of the slumping dynamics to variations in the initial condition is explored in figure 3. We observe that, for initial profiles ‘above’ the arrested state, early-time discrepancies are eventually forgotten as the slumping trajectories relax to the same late-time shape. The heights of the yield surfaces and free surface at $t=5\pi /4$ are plotted in figure 4. Early-time contact point asymmetry is visible in the $\alpha =20$ row, because $\delta \alpha =1$ and so the early-time regime had subsided. The height of the yield surfaces and the free surface at $\eta = \eta _{\kern-1pt f}{(0)}/2$ during an early period of oscillation are presented in figure 5. The regimes observed in figure 5 are qualitatively similar to those described by Hinton et al. (Reference Hinton, Collis and Sader2023a ) for a horizontal free surface, although the addition of free-surface gradients introduces a spatial dependence that may suppress or support the growth of a yield surface in some region.

For large $B$ , $h$ remains close to $h(\eta ,0)$ and the yield surface is confined to a thin region above the plate. For $B \gg 1$ , the deposit only yields in $\eta \leqslant 0$ when $\cos {t}\geqslant 0$ , and only yields in $\eta \geqslant 0$ when $\cos {t} \leqslant 0$ , because the oscillatory forcing is no longer large enough to overcome the yield stress when their signs are opposing (see figure 4 for $\alpha =0.05, 0.5$ and $B=0.8, 1.5$ ).

When $\alpha \gg 1$ and $B \geqslant \mathcal{O}(1)$ , momentum diffusion is slow relative to the plate oscillation, so the shear stress decays rapidly away from the plate. In this regime, the flow within the deposit is envisioned as a large plugged bulk moving atop a thin lubricating layer at the base (see figure 5 for $\alpha =5, 20$ and $B=0.8, 1.5$ ). The large rigid layer permits instantaneous momentum transfer to the free surface, allowing the velocity of the plug to adopt the $\mathcal{O}(1)$ values of the basal velocity and thus evolve the free surface more rapidly than the other regimes.

Generally, smaller values of $B$ result in a larger proportion of the material becoming yielded at early times. When $\alpha \gg 1$ and $B \lt \mathcal{O}(1)$ the velocity dissipates rapidly within the yielded region adjacent to the plate, dampening momentum transfer to the free surface. The thin rigid regions emanating from the plate take longer to reach the free surface, so multiple yield surfaces may exist in the bulk of the deposit if $B$ is sufficiently small (for example, see figure 5 for $\alpha =5, B=0.05$ ; for $\alpha =20$ , multiple yield surfaces can be observed for $B=0.05$ and $B=0.3$ ). This only occurs where the hydrostatic pressure gradients are negligible, however; in a region about the contact points the hydrostatic pressure gradients are larger than the oscillatory forcing, so the inertia is relatively small and there may only be one yield surface there at any point in time. Thus, simultaneous existence of regions with multiple yield surfaces, and regions with one yield surface can be observed in the deposit (see figure 4 for $\alpha =5, B=0.05$ ). This phenomenon creates closed regions of yielded material that propagate upwards through the deposit, which shrink and disappear as the free surface is approached. Convective inertia can be significant in this regime.

4. Low-frequency regime

In this section, the low oscillatory Reynolds number regime, wherein $\alpha \ll 1$ , is analysed to elucidate characteristic features of both the internal fluid motion and evolution of the free surface.

In this regime, inertia is small and momentum quickly diffuses from the plate, and there is at most one yield surface for each $\eta$ coordinate at any time, the uppermost yield surface which we denote $Y$ . The velocity field may then be calculated by partitioning the flow domain and momentum (2.11) into its constituent yielded state below $Y$ and rigid state above $Y$

(4.1a)

(4.1b)

In the limiting case where $\alpha =0$ , ${\partial \mathscr{u}}/{\partial z}=0$ in both the yielded and rigid regions, enforcing that $ \mathscr{u}=0$ and $w=0$ everywhere, i.e. the entire deposit moves in concert with the plate. This motivates an asymptotic expansion of the form

(4.2) \begin{align} \mathscr{u} = \alpha \mathscr{u}_1 + \mathcal{O}(\alpha ^2), \qquad w = \alpha w_1 + \mathcal{O}(\alpha ^2), \qquad Y=Y_{0} + \mathcal{O}(\alpha ). \end{align}

Substituting $\mathscr{u}=0$ into the uppermost yield surface (2.15) reveals an expression for $Y_{0}$

(4.3) \begin{align} Y_{0}= \text{max}\left \{0, \, h-\frac {B}{\left|{B\frac {\partial h}{\partial \eta }+\cos t} \right|} \right \}. \end{align}

Similarly to the flow considered in Hinton et al. (Reference Hinton, Collis and Sader2023a ), momentum propagates instantaneously in the $\alpha =0$ limit, and the material is thus immediately yielded upon commencement of oscillation.

The $\mathcal{O}(\alpha )$ velocity field in the yielded region may be determined by integrating (4.1a ), noting that the inertial terms on the left-hand side are negligible for $\delta \alpha ^2 \ll 1$ . Evaluating $\mathscr{u}$ at $z=Y$ and enforcing its continuity there then gives the velocity of the plugged region above it. Altogether,

(4.4a)

(4.4b)

Equation (4.4) provides an extension to the analogous formula of Hinton et al. (Reference Hinton, Collis and Sader2023a ) in the presence of free-surface gradients. Equation (4.4) may be substituted into the mass conservation (2.13) to close the system and determine the evolution of the free surface

(4.5) \begin{align} \frac {\partial h}{\partial t} = \alpha \delta \frac {\partial }{\partial \eta } \left [\left (\frac {h Y_{0}^2}{2} - \frac {Y_{0}^3}{6}\right )\left (B\frac {\partial h}{\partial \eta }+ \cos t\right )\right ] + \mathcal{O}(\alpha ^2). \end{align}

This equation is integrated numerically using the method of lines and a fourth-order Runge–Kutta method. Results for the long-time evolution of the deposit for $t=1, 10, 100, 1000, 10\,000$ are presented in figure 2. The flux driving the motion of the free surface is at most $\mathcal{O}{(\alpha )}$ and this slumping occurs slowly.

For $B \ll 1$ , the deposit is thin and the hydrostatic pressure gradients are weak, so the oscillatory forcing dominates the flow and $\mathscr{u}_R=\mathcal{O}{(\alpha )}$ . The yield surface periodically approaches the free surface at sufficiently early times; far from the contact points $Y_{0}=h-\mathcal{O}{(B)}$ , except at times that satisfy $|{\cos {t}}|=\mathcal{O}{(B)}$ . This can be clearly observed in figure 5 when $B=0.05$ . For $B \gg 1$ , hydrostatic pressure gradients are large and suppress the inertial forcing, so $Y_0=\mathcal{O}{(1/B)}$ and $\mathscr{u}_R=\mathcal{O}{(\alpha /B)}$ .

5. The arrested state

Figure 2 indicates the existence of an arrested profile which the deposit slumps toward. In this section it is explicitly calculated.

5.1. Criterion for the arrested state

When the material is rigid everywhere, it moves as a plug with velocity $\sin {t}$ (i.e. it is stationary in the $\eta$ frame). Differentiating the momentum (2.11) with respect to $z$ reveals that the shear stress is linear in $z$ throughout an entirely rigid deposit

(5.1) \begin{align} \frac {\partial ^2 \tau _R}{\partial z^2}=0. \end{align}

Integrating this equation with respect to $z$ and applying the boundary conditions (2.14) then gives the expression for the shear stress within a rigid layer of thickness $h$

(5.2) \begin{align} \tau _R=\left (z-h\right ) \left (B \frac {\partial h}{\partial \eta } + \cos t\right )\!. \end{align}

The shear stress within a rigid deposit $\tau _R$ takes on its maximum magnitude at the plate $z=0$ . The magnitude of the basal rigid stress $|{\tau _R}|_{z=0}$ must remain smaller than the yield stress, which provides a criterion for an arrested region within the deposit

(5.3) \begin{align} B \geqslant h \, \left |{B \frac {\partial h}{\partial \eta } + \cos t}\right |. \end{align}

If this condition holds throughout an entire oscillation and across the entire deposit, then the deposit has reached its arrested state. A corollary of (5.3) is that the deposit initially yields where $\eta \lt \text{max}{\{0, \eta _{\kern-1pt f}{(0)} - 2B^2\}}$ .

5.2. Calculation of the arrested state

The magnitude of the basal stress (5.2) attains its maxima when $\cos t = \pm 1$ , i.e. at $t= n \pi$ . Enforcing that the maximum basal stress is never above the yield stress provides a criterion for a fully rigid profile

(5.4) \begin{align} B \geqslant h \left |{B \frac {\partial h}{\partial \eta } + \text{sgn}{\left (\frac {\partial h}{\partial \eta }\right )}}\right |. \end{align}

Any profile which satisfies (5.4) is a fully rigid state. When the deposit is not released with significant inertia, it approaches one unique final arrested state which is symmetric about $\eta =0$ , and is denoted $h_{\infty } {(\eta )}$ . This final state is defined by equality of (5.4)

(5.5) \begin{align} B = h_{\infty } \left |{B \frac {\partial h_{\infty }}{\partial \eta } + \text{sgn}{\left (\frac {\partial h_{\infty }}{\partial \eta }\right )}}\right |. \end{align}

This equation may be integrated to obtain an implicit expression for the arrested profile in terms of the contact points $\pm \eta _{f \infty }$ , where $h_{\infty } (\pm \eta _{f \infty })=0$

(5.6) \begin{align} {|{\eta }| - \eta _{f \infty }}=B^2 \log \left (1-\frac {h_{\infty }}{B}\right ) + B h_{\infty }. \end{align}

An explicit formula is available in terms of principal branch of the Lambert-W function $W_0 (\boldsymbol{\cdot })$ , namely, $h_{\infty }=B+ B W_0{ [- \exp (-1-\{\eta _{f \infty }-|{\eta }|\}/B^2) ]}$ . The final contact points are determined by enforcing mass conservation throughout the deposit

(5.7) \begin{align} \int _{0}^{h_{\infty }(0)} \eta \, d h_{\infty } = \frac {1}{2}, \end{align}

by which an implicit formula for these contact points may be derived

(5.8) \begin{align} \frac {\eta _{f \infty }}{B^2}+\frac {1}{B}\sqrt {2\eta _{f \infty } -\frac {1}{B}} +\log {\left (1-\frac {1}{B}\sqrt {2\eta _{f \infty } -\frac {1}{B}}\right )}=0. \end{align}

Figure 6.

The arrested profile $h_{\infty }(\eta )$ , the solution to (5.6), for $B=0.15, 0.3, 0.5, 0.8, 2$ in solid blue. The limiting cases of $B \rightarrow \infty$ and $B=0$ are presented in the dashed blue lines.

The arrested profile $h_{\infty } {(\eta )}$ is plotted for various values of $B$ in figure 6. Explicit asymptotic formulae may be derived for $\eta _{f \infty }$ in the limits of small and large $B$ ; expansion of (5.8) in these limits yields

(5.9) \begin{align} \eta _{f\infty } \sim \begin{cases} \frac {1}{2 B} & B \ll 1 ,\\[9pt] \left (\frac {9}{32}\right )^{1/3} + \frac {1}{8 B} & B \gg 1. \end{cases} \end{align}

The maximum height of the arrested state may be obtained by expansion of (5.6)

(5.10) \begin{align} h_{\infty } {(0)} \sim \begin{cases} B - B \exp \left (-1-\frac {1}{2 B ^3}\right ) & B \ll 1 ,\\[9pt] \left (\frac {3}{2}\right )^{1/3} - \left (\frac {3}{16}\right )^{2/3}\frac {1}{B} & B \gg 1. \end{cases} \end{align}

The arrested contact point $\eta _{f \infty }$ and maximum height $h_{\infty } {(0)}$ , along with their asymptotic expressions, are plotted as functions of $B$ in figure 7. For large values of $B$ , inertial forcing is suppressed by strong hydrostatic pressure gradients, resulting in only a small deviation from the initial profile. For small values of $B$ the arrested profile tends towards a flat, uniform profile, supported near the edges by a thin region where the free-surface gradients are large enough to block the deposit from spreading. The cusp point at the origin also becomes less pronounced as $B$ decreases, owing to the smaller discontinuity in the shear stress and thus in the free-surface gradients required to keep the stress at or below $B$ .

Additionally, no specific form of constitutive law in the yielded region is required for this derivation, so this arrested shape is expected for any general viscoplastic material.

Figure 7.

(a) Contact point $\eta _{f \infty }$ of the arrested slump and (b) the maximum height of the arrested slump, $h_{\infty } {(0)}$ , shown as functions of $B$ . The implicit solutions to (5.8) and (5.6) are plotted in (solid) blue, respectively. The red (dashed) lines are the asymptotic expressions for small and large $B$ , (5.9) and (5.10).

5.3. Equivalence to the gravity-driven arrested state on a slope

Figure 8.

The equivalence between (a) the arrested state due to oscillation $u=U_0 \sin {(\omega t)}$ and (b) the gravitationally slumped arrested state on a triangular bed angled $\theta =\arcsin {({U_0 \omega }/{g})}$ to the horizontal. The height of each deposit is exaggerated for clarity.

Liu & Mei (Reference Liu and Mei1989) showed that a shallow deposit of Bingham material on a plane inclined at angle $\theta$ to the horizontal, with $x$ denoting the streamwise coordinate along the plane, has a down-slope arrested state $h(x)$ given by

(5.11) \begin{align} x-x_{\infty }=\frac {\tau _Y}{\rho g \sin ^2 \theta } \log {\left (1-\frac {\rho g h \sin \theta }{\tau _Y} \right )}+\frac {h}{\sin \theta }, \end{align}

in dimensional terms. The down-slope contact point $x_{\infty }$ is calculated by enforcing mass conservation; for the purposes of this subsection, $\int _0^{x_{\infty }} h \, d x=A/2$ ; see figure 8.

Curiously, the down-slope gravity-driven arrested profile (5.11) is the same shape as the oscillation-forced arrested state (5.6) for $\eta \geqslant 0$ (in dimensional terms). By symmetry, this set-up may be likened to a deposit on a triangular bed as shown in figure 8, where the bed is angled appropriately to the horizontal. Indeed, for the specific slope angle

(5.12) \begin{align} \sin \theta = \frac {U_0 \omega }{g}, \end{align}

which depends only on the parameters of the basal plate oscillation and gravity, (5.11) and (5.6) are exactly matched. This equivalence occurs because the shear stress associated with the hydrostatic pressure gradients of the slump on an inclined plane

(5.13) \begin{align} \tau =\rho g \left (z-h\right ) \left ( \frac {\partial h}{\partial x} -\sin {\theta }\right )\!, \end{align}

takes the exact form of the maximum rigid stress of the oscillatory problem for $\eta \geqslant 0$ , which is given by (5.2) in dimensional variables

(5.14) \begin{align} \tau _R=\rho g \left (z-h_{\infty }\right ) \left ( \frac {\partial h_{\infty }}{\partial x} - \frac {U_0 \omega }{g}\right )\!. \end{align}

Since (5.11) is valid for $\sin {\theta } \ll 1$ , this equivalence requires that $U_0 \omega /g$ is small, or equivalently, that $B \gg \varepsilon$ , which is implicit within this paper.

It is worth noting that this equivalence is only of the arrested states between these two problems; the different physical processes guiding the respective approaches to the arrested state have fundamentally different decay properties, which will be explored in § 6.

6. Approach to the arrested state

In this section, the late-time dynamics of a deposit close to the symmetric arrested state are analysed. A scaling argument for the temporal decay of a deposit of Herschel–Bulkley fluid is presented in Appendix B, which is found to be $h-h_{\infty } \sim (\delta \alpha t)^{-2N /(N+2)}$ . Consequentially, the contact point error $\eta _{f \infty } \mp \eta _{\kern-1pt f}^{\pm } \sim (\delta \alpha t)^{-2N /(N+2)}$ to conserve mass. The late-time solution is explicitly derived in this section for the Bingham fluid ( $N=1$ ) case.

At late times inertia is negligible, so the low- $\alpha$ formulation of § 4 applies to this regime. In addition, the deposit is close to its arrested state, so the material only yields for short time periods around the times of maximum basal stress when $\cos t=\pm 1$ , and is rigid throughout the remaining oscillation period. A fast time variable $T$ , where $|{T}| \ll 1$ , is introduced to describe this period

(6.1) \begin{align} t= n \pi + T, \end{align}

and $n \gg 1$ . The contact points are expanded about their arrested states in a small perturbation of the ‘effective time’ parameter $\delta \alpha t$

(6.2) \begin{align} \eta _{\kern-1pt f} (t)=\eta _{f \infty } - (\delta \alpha t)^{-2/3} \tilde {\eta }_f. \end{align}

In doing so, it is assumed that any asymmetries in the positions of the left and right contact points occur at higher asymptotic orders. A strained spatial coordinate $\xi$ is introduced

(6.3) \begin{align} \xi =\frac {\eta }{\eta _{\kern-1pt f} (t)}=\frac {\eta }{\eta _{f \infty } - (\delta \alpha t)^{-2/3}\tilde {\eta }_f}, \end{align}

where $h=0$ at $\xi =\pm 1$ . For ease of calculation, the arrested state is taken to be $H_{\infty }{(\xi )} :=h_{\infty }{(\eta _{f \infty } \xi )}$ in $\xi$ coordinates, noting that $H_{\infty }$ is the true arrested state with $\tilde {\eta }_f=0$ . A perturbation is taken about this arrested state

(6.4) \begin{align} h (\xi , t)=H_{\infty }(\xi ) + (\delta \alpha t)^{-2/3} \tilde {h} (\xi ), \end{align}

which is valid sufficiently far from the contact points. A condition for the deposit to yield within these short time periods is obtained by substituting (6.1–6.4) into the criterion for yielding (5.3), expanding in each small parameter and retaining only the first terms in each expansion

(6.5) \begin{align} \begin{split} H_{\infty }^2 \left [\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } + (-1)^n\right ]^2 + 2 (\delta \alpha n\pi )^{-2/3} H_{\infty } \left [\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } + (-1)^n\right ] \times & \\\left \{ \tilde {h} \left [ \!\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } + (-1)^n \!\right ] + H_{\infty } \! \left [ \! \frac {B}{\eta _{f \infty }} \frac {\text{d} \tilde {h}}{\text{d}\xi } + \frac {B \tilde {\eta }_f }{\eta _{f \infty }^2} \frac {\text{d} H_{\infty }}{\text{d}\xi } - (-1)^n(\delta \alpha n\pi )^{2/3} \frac {T^2}{2}\right ] \right \}& \gt B^2.\end{split} \end{align}

This equation provides two key pieces of information. Firstly, by rewriting (5.5) as

(6.6) \begin{align} H_{\infty } \left [\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } - \text{sgn}(\xi )\right ]=-B \text{sgn}(\xi ), \end{align}

it can be seen that the deposit only yields when $ (-1 )^n=-\text{sgn}(\xi )$ . That is, when $n$ is even, the plate acceleration is positive, so yield surfaces only appear in the $\xi \lt 0$ plane, and likewise when $n$ is odd. This regime is therefore akin to the large $B$ limit where the free-surface profile is always close to its arrested state. Any yielding at the origin $\xi =0$ is due to higher-order asymmetry, and so the positive and negative $\xi$ planes are disconnected into independent problems for odd and even $n$ , respectively. Secondly, bounds upon $T$ may be calculated, which describe the window of time in which the fluid yields

(6.7) \begin{align} |{T}| \lt \tilde {T}_n := (\delta \alpha t)^{-1/3}\sqrt {2 B \left [\frac {\tilde {h}}{H_{\infty }^2}+\frac {(-1)^n}{\eta _{f\infty }} \left (\frac {\text{d} \tilde {h}}{\text{d}\xi } +\frac { \tilde {\eta }_f}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi }\right )\right ]} + o{\left (t^{-1/3}\right )}, \end{align}

where $t$ and $n \pi$ may be used interchangeably to leading order. So, temporal windows of yielding are appropriately diminishing in time; as $h$ tends towards its arrested state, $\tilde {T}_n \sim (\delta \alpha t)^{-1/3}$ . The windows are symmetric about $t=n \pi$ to leading order. Clearly, as the front is approached and free-surface gradients become large, $\tilde {T}_n= O(1)$ and the linearisation (6.1) is no longer valid. The free surface is analysed independently in this ‘boundary layer’ in Appendix C. The remaining analysis in this section takes place sufficiently far from the front.

Identically to the $\alpha \ll 1$ regime, only one yield surface may exist at late times, which we denote $Y$ . Substitution of (6.1–6.4) into the expression for the yield surface (4.3) gives

(6.8) \begin{align} Y=(\delta \alpha t)^{-2/3} \tilde {h} +\frac {H_{\infty }^2}{B}\left [ \frac {B(-1)^n (\delta \alpha t)^{-2/3}}{\eta _{f\infty }} \left (\frac {\text{d} \tilde {h}}{\text{d}\xi }+\frac {\tilde {\eta }_f}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi }\right )-\frac {T^2}{2}\right ], \end{align}

for $|{T}| \lt \tilde {T}_n$ , so velocity variations are confined to a thin region of order $z=\mathcal{O}{(t^{-2/3})}$ near the plate. Further substitution of (6.1–6.4) and (6.8) into the low- $\alpha$ velocity expression (4.4) gives the velocity field, partitioned by $Y$ , to be

(6.9a) \begin{align} \mathscr{u}&=\alpha \left (\frac {z^2}{2} - Y z\right ) \left (\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } + (-1)^n\right ) & \kern-50pt \text{for} \quad z\lt Y&, \\[-12pt] \nonumber \end{align}

(6.9b) \begin{align} \mathscr{u}_R&=-\frac {\alpha Y^2}{2}\left (\frac {B}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi }+(-1)^n\right ) & \kern-50pt \text{for} \quad z \geqslant Y&. \\[10pt] \nonumber \end{align}

Finally, substituting (6.1–6.4) and (6.8) into mass conservation (4.5) and once again utilising (6.6) then provides the evolution of the free-surface profile

(6.10) \begin{align} \frac {\partial }{\partial t}h(\eta (\xi ,t),t) &= \frac { B \delta \alpha (-1)^{n}}{2 \eta _{f \infty }}\frac {\partial }{\partial \xi } Y^2 \\[-12pt] \nonumber \end{align}

(6.11) \begin{align} \frac {\partial }{\partial t} \left [(\delta \alpha t)^{-2/3}\tilde {h} \right ] +\frac {\partial \xi }{\partial t} \ \frac {\mathrm{d} H_{\infty }}{\mathrm{d} \xi } &= \frac { B \delta \alpha (-1)^{n}}{2 \eta _{f \infty }}\frac {\partial }{\partial \xi } Y^2. \\[10pt] \nonumber \end{align}

Only the rigid component of the velocity contributes because the flux from the yielded region is a higher-order correction. Importantly, the free surface may only evolve within the small time periods of yielding $|{T}| \lt \tilde {T}_n$ , which elucidates a separation of time scales that may be exploited analytically by averaging over the fast $T$ time scale. To this end, the period average of (6.11) about $t=n \pi$ is obtained by integrating it from $(n-1) \pi$ to $(n+1) \pi$ . By symmetry, only one side of the deposit is considered; $\xi \lt 0$ (if $n$ is even) or $\xi \gt 0$ (if $n$ is odd) is fixed and integrated independently of the other side. On the left-hand side of (6.11), we obtain

(6.12) \begin{align} \frac {1}{2 \pi }\int _{(n-1) \pi }^{(n+1)\pi } \frac {\partial }{\partial t}h(\eta (\xi ,t),t) \, \text{d}t &=\frac {1}{2 \pi } \! \int _{(n-1) \pi }^{(n+1)\pi }\frac {\partial }{\partial t} \! \left [ \!(\delta \alpha t)^{-2/3}\tilde {h} +\frac {\xi \tilde {\eta }_f (\delta \alpha t)^{-2/3} }{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } \!\right ] \text{d}t \nonumber\\[3pt]&=-\frac {2}{3} (\delta \alpha )^{-2/3} t^{-5/3}\left (\tilde {h} +\frac {\xi \tilde {\eta }_f }{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi }\right )\nonumber\\[3pt]& \quad + o{\left [(\delta \alpha )^{-2/3} t^{-5/3}\right ]}. \end{align}

We have used $t$ and $n \pi$ interchangeably, and have retained only the large $n$ asymptotic. On the right-hand side, the integration bounds may be simplified to only include times of yielding where the flux is non-zero. Utilising (6.8), we can take the time average of the right-hand side of (6.11) to obtain

(6.13) \begin{align} \begin{split} \frac {1}{2 \pi } \int _{(n-1) \pi }^{(n+1)\pi } \frac { B \delta \alpha (-1)^{n}}{2 \eta _{f \infty }}\frac {\partial }{\partial \xi } Y^2 \, \text{d}t &= \frac {1}{2 \pi } \int _{n \pi - \tilde {T}_n}^{n \pi + \tilde {T}_n} \frac { B \delta \alpha (-1)^{n}}{2 \eta _{f \infty }}\frac {\partial }{\partial \xi } Y^2 \, \text{d}t \\&= \frac {B \delta \alpha (-1)^{n}}{4 \pi \eta _{f \infty }}\frac {\partial }{\partial \xi } \int _{- \tilde {T}_n}^{ \tilde {T}_n} Y^2\, dT\\= - \frac {4 \sqrt {2} B (\delta \alpha )^{-2/3} t^{-5/3} (-1)^{n+1}}{15 \pi \eta _{f \infty }}& \frac{\text{d}}{\text{d}\xi } \left \{ \frac {H_{\infty }^4}{B^2} \left [\! \frac {B \tilde {h}}{H_{\infty }^2} + \frac {B (-1)^n}{\eta _{f\infty }}\! \left ( \! \frac {\text{d} \tilde {h}}{\text{d}\xi }+\frac {\tilde {\eta }_f}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } \!\right ) \! \right ]^{5/2}\right \} \!. \end{split} \end{align}

Equating (6.12) and (6.13) reveals the spatial equation for $\tilde {h}$ and $\tilde {\eta }_f$

(6.14) \begin{align} \tilde {h} +\frac {\xi \tilde {\eta }_f }{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi }=\frac {2 \sqrt {2} B (-1)^{n+1}}{5 \pi \eta _{f \infty }} \frac{\text{d}}{\text{d}\xi } \! \left \{ \frac {H_{\infty }^4}{B^2} \left [\frac {B \tilde {h}}{H_{\infty }^2} + \frac {B (-1)^n}{\eta _{f\infty }} \left (\frac {\text{d} \tilde {h}}{\text{d}\xi }+\frac {\tilde {\eta }_f}{\eta _{f \infty }} \frac {\text{d} H_{\infty }}{\text{d}\xi } \! \right ) \! \right ]^{5/2}\! \right \}\!. \end{align}

The form of (6.14) can be simplified slightly by observing that every $\xi$ derivative may be written as a derivative with respect to $H_{\infty }$ , and under this transformation (6.14) is

(6.15) \begin{align} \begin{split} &\tilde {h}+\frac {\tilde {\eta }_f}{\eta _{f \infty }} \frac {H_{\infty } -B}{H_{\infty }}\left [B \log {\left (1-\frac {H_{\infty }}{B}\right ) + H_{\infty } + \frac {\eta _{f \infty }}{B}}\right ]\\ &= \, \frac {2 \sqrt {2}}{5 \pi } \frac {H_{\infty } -B}{H_{\infty }} \frac{\text{d}}{\text{d}H_{\infty }} \left \{ \frac {H_{\infty }^4}{B^2} \left [\frac {B \tilde {h}}{H_{\infty }^2} + \frac {B- H_{\infty }}{H_{\infty }} \left (\frac {\text{d} \tilde {h}}{\text{d}H_{\infty }}+\frac {\tilde {\eta }_f}{\eta _{f \infty }}\right )\right ]^{5/2}\right \}. \end{split} \end{align}

The first boundary condition is obtained by enforcing that $h$ , and thus $\tilde {h}$ , must vanish at $H_{\infty }=0$ . It is shown in Appendix C that the dynamics near the contact point does not affect the leading-order problem for $\tilde {h}$ , and that this choice of boundary condition is justified. The second condition is obtained by substituting (6.2) and (6.4) into the volume conservation condition $\int _{-\eta _{\kern-1pt f}(t)}^{\eta _{\kern-1pt f}(t)} h \, \text{d}\eta =1$ . Respectively, these are

(6.16) \begin{align} \tilde {h}{(0)}=0, \qquad \int _{0}^{H_{\infty }(0)} \frac {\tilde {h} H_{\infty }}{B-H_{\infty }}\, \text{d}H_{\infty }=\frac {\tilde {\eta }_f}{2 B \eta _{f \infty }}. \end{align}

Figure 9.

(a) The scaled shape function $\tilde {h}{(\xi )}/(\tilde {\eta }_f/\eta _{f \infty }^2)$ and (b) the scaled auxiliary shape function $\mathcal{I}(H_{\infty })/(\tilde {\eta }_f H_{\infty }(0)/B)$ (plotted on a normalised domain), for ten logarithmically spaced values of $B$ between $0.1$ and $10$ .

Figure 10.

The contact point constant $\tilde {\eta }_f$ (6.2) in red – obtained from numerical integration of (D2) with respect to the boundary conditions (D3) – as a function of $B$ for $B \in [0.05, 100]$ . The dot-dashed lines are numerical approximations of the asymptotic forms of $\tilde {\eta }_f$ fitted to the red curve, namely $\tilde {\eta }_f=0.237/B$ for large $B$ (in black), and $\tilde {\eta }_f=0.423/B^3$ for small $B$ (in gold).

In addition, $\tilde {\eta }$ is directly solved for by taking the limit of (6.15) as $H_{\infty } \rightarrow 0$ . The details of the numerical integration of (6.15–6.16) may be found in Appendix D. The resulting plots of $\tilde {h}$ and $\mathcal{I}$ , the auxiliary shape function, are shown in figure 9. As $B$ decreases towards zero, the maximum of the shape perturbation traverses towards the contact point. Simultaneously, the value of the perturbation at the origin $\xi =0$ tends towards zero. Thus, when $B \ll 1$ , we expect $h$ to have flattened near the origin at earlier times, while away from the centre, the deposit is still deforming as the arrested state is approached. Observing figure 2 for $B=0.3$ and $\alpha =0.5, 5$ confirms this expectation. The plot of the resulting contact point constant $\tilde {\eta }_f$ is exhibited in figure 10. It may be observed that $\tilde {\eta }_f = \mathcal{O}{(1/B)}$ when $B \gg 1$ , and $\tilde {\eta }_f = \mathcal{O}{(1/B^3)}$ when $B \ll 1$ . The asymptotic expressions for $\tilde {h}(0)$ and $\tilde {\eta }_f$ are plotted against the full numerical solutions in time in figure 11. This figure delineates the presence of the early ‘waiting time’ regime – $\tilde {\eta }_f$ is virtually frozen when $t \leqslant \mathcal{O}{(1/\delta \alpha )}$ – and the late-time $(\delta \alpha t)^{-2/3}$ scaling is clearly observed.

Figure 11.

The late-time perturbations for $t \in [10, 10\,000]$ of (a) the height at the origin $(\delta \alpha t)^{-2/3}\tilde {h}(0)=h(0,t)-H_{\infty }(0)$ , and (b) the contact points $ (\delta \alpha t)^{-2/3} \tilde {\eta }_f= \eta _{f \infty } - |{\eta _{\kern-1pt f}^{\pm } (t)}|$ . Plots are for parameter values $\alpha =0.5, B=0.8, \delta =0.05$ . The numerical solutions are presented in solid black lines. Since $\delta \alpha =0.025$ , the contact points are only weakly asymmetric and are virtually indistinguishable on the scale in (b). The inset of (b) depicts the asymmetry of the contact points at early times as they deviate from the starting position: the left contact point $\eta ^- (t)$ is in grey, and the right contact point $\eta ^+(t)$ is in black. The red dotted lines are the corresponding asymptotic expressions at late times (6.4) and (6.2), obtained from numerical integration of (D2) with respect to the boundary conditions (D3).

Figure 12.

The spatial component of the free-surface perturbation $\tilde {h}$ (6.4). The numerical solutions for $\tilde {h}(\xi )=(\delta \alpha t)^{2/3}[h(\xi ,t)-H_{\infty }(\xi )]$ for $\alpha =0.5, B=0.8, \delta =0.05$ and $t=10, 100, 1000, 10\,000$ are presented in solid black lines. A lower value of $\alpha =0.05$ is presented in the grey dashed line. These are compared with the red dashed line – the solution to (6.14) – obtained from numerical integration of (D2) with respect to the boundary conditions (D3). Note that the full numerical solution for $h$ mispredicts the free surface curvature in the vicinity of the contact points (see Appendix C), and so the perturbation $\tilde {h}$ is erroneously calculated to be negative there. Only positive values of $\tilde {h}$ are shown in this figure.

The shape function $\tilde {h}$ is plotted against the full numerical solution in space in figure 12. The $\alpha =0.5$ curves for $t=10, 100, 1000$ are almost indistinguishable from the $\alpha =0.05$ curves for $t=100, 1000, 10\,000$ , respectively. This is because the effective time $\delta \alpha t$ is the same for the respective curves. The asymptotic expression for the yield surface at the station $\xi =1/2$ is presented against the full numerical simulation in figure 13, and as expected, the yield surfaces only occur for odd multiples of $\pi$ .

Figure 13.

The yield surface at $\xi =1/2$ at late times, for $\alpha =0.5, B=0.8, \delta =0.05$ . The black dot-dash curve represents the numerical result, and the red dashed curve is the asymptotic prediction obtained through substitution of $\tilde {h}$ and $\tilde {\eta }_f$ into (6.8). The yield surfaces only appear for (progressively shorter) times about $n \pi$ when $n$ is odd for $\xi \gt 0$ , as expected.

Hogg & Matson (Reference Hogg and Matson2009) showed that, for slumping on an inclined plane, the free-surface profile of a Bingham fluid asymptotically approaches its arrested state as $t^{-1}$ . For oscillatory forcing, $h-h_{\infty } = \mathcal{O}(t^{-2/3})$ , which is a slower decay towards the same arrested state. This difference can be attributed to the periodically varying stress in the oscillatory problem; when the free-surface profiles are close to the final shape, both the inclined deposit and oscillating deposit have the same leading-order maximum shear stress profile. In the inclined problem, down-slope gravitational forces impart this maximum stress on the deposit invariably in time. However, the oscillatory deposit only attains this stress transiently, yielding in the vicinity of this maximum stress for progressively shorter times about $t=n \pi$ as the arrested state is approached.