Obstructed free-surface viscoplastic flow on an inclined plane

Abstract

The interaction of steady free-surface flows of viscoplastic material with a surface-piercing obstruction of square cross-section on an inclined plane is investigated theoretically. The flow thickness increases upstream of the obstruction and decreases in its lee. The flow depends on two dimensionless parameters: an aspect ratio that relates the flow thickness, the obstruction width and the plane inclination; and a Bingham number that quantifies the magnitude of the yield stress relative to the gravitationally induced stresses. Flows with a non-vanishing yield stress always form a static ‘dead’ zone in a neighbourhood of the upstream and downstream stagnation points. For relatively wide obstructions, a deep ‘ponded’ region develops upstream with a small dead zone, while the deflected flow reconnects over relatively long distances downstream. The depth of the upstream pond increases with both the dimensionless yield stress and width of the obstruction, while the unyielded dead zone varies primarily with the yield stress. Both are predicted asymptotically by balancing the volume flux of fluid into and out of the ponded region. When the obstruction is narrow, the perturbation to the depth of the oncoming flow is reduced. It exhibits fore–aft antisymmetry, while the dead zone is symmetric to leading order. Increasing the yield stress leads to larger dead zones that eventually encompass all of the upstream- and downstream-facing boundaries of the obstruction and fully divert the flow. Results for obstructions with circular and rhomboidal cross-sections are also presented and illustrate the effects of boundary shape on the properties of the steady flow.

1. Introduction

Free-surface flows of fluids that have a yield stress are common in the environment. Examples include lava flows, fluidised mine tailings, mud slides and avalanches (Coussot & Proust 1996; Fink & Griffiths 1998; Hu & Bürgmann 2020). Often these flows interact with topography and obstructions such as trees and constructed barriers (Tai et al. 2001; Barberi et al. 2003; Dietterich et al. 2015; Bernabeu, Saramito & Harris 2018; Chevrel et al. 2019). However, the all-important role of the yield stress that influences these interactions is not well-understood. In this paper, we analyse the steady interaction between a free-surface viscoplastic flow (with negligible inertia) and a cylindrical obstruction on an inclined plane (see figure 1a).

Figure 1. (a) Schematic of the steady free-surface flow. (b) Schematic of the flow along the centreline indicating the kinematically plugged region, which comes to encompass the entire layer near the stagnation points. (c) Steady flow thickness for $B=1$ and $\mathcal {L}=1$. (d) Corresponding height of the yield surface, $Y$.

There has been substantial research on Newtonian free-surface flows around obstructions, which provides a foundation for analysing the viscoplastic case. Investigating these flows is challenging owing to their nonlinear governing equations, even when a lubrication approximation is applied. This has motivated a multitude of numerical methods to simulate the motion including boundary integral formulations and finite element techniques (Baxter et al. 2009; Sellier et al. 2009). These works focused on very thin films, in which the effects of interfacial surface tension play a key role. Some of the observed features, such as a capillary ridge upstream of a step-down in the substrate, appear only when surface tension is dominant (Kalliadasis, Bielarz & Homsy 2000). However, other results, including the horseshoe-like flow deflection, are universal to the interaction of slow flows with obstructions, even when gravity dominates surface tension.

Steady viscous free-surface flows around cylinders of various cross-sections have recently been studied by Hinton, Hogg & Huppert (2020). They showed that the width of the cylinder relative to the thickness of the oncoming flow (and the slope inclination) was the key parameter determining the nature and magnitude of the interaction, and they constructed asymptotic solutions in both the regimes of very narrow and very wide obstructions. The latter can be encapsulated in a simple balance of volume fluxes that enables the prediction of the upstream flow thickness in the relatively deep flow just upstream of the obstruction. The key idea was that the flux in the normal direction to the obstruction vanishes throughout the upstream deep region resulting in a ‘pond’-like free-surface.

In this contribution, we focus on steady flows of viscoplastic materials, which are characterised by a ‘yield stress’. For deviatoric stresses below the yield stress, the material is rigid and when the yield stress is exceeded, the material yields and flows as a fluid albeit with an adjusted relationship between the stress and rate of strain tensor. Free-surface viscoplastic flows generally consist of two distinct regions. Near the substrate over which the material flows, there is typically a sufficiently large shear stress so that the material yields with a significant velocity gradient perpendicular to the substrate. In the upper portion of the layer the material is rigid since the stress vanishes at the free surface. The layer is partitioned into a yielded region below a rigid region with the boundary known as the ‘yield surface’. This simple description can require a subtle alteration to be compatible with lubrication theory. In cases where the thickness varies along the flow, the stress in the upper region just exceeds the yield stress with the vertical velocity gradient being negligible but non-zero (Balmforth & Craster 1999). In this case, the upper region behaves kinematically as a rigid (pseudo) plug but its slight yielding allows extensional variation in the flow. It is also worth noting that when the (pseudo) plug extends down to the substrate, the entire layer is truly rigid.

Previously, obstructed viscoplastic flows have mostly been studied in the context of simplified geometries. Examples include flow in the narrow gap between two plates with an obstacle present in the conduit, known as a Hele-Shaw cell (Hewitt et al. 2016; Daneshi et al. 2020), and two-dimensional flow past a cylinder (Mitsoulis 2004; Tokpavi, Magnin & Jay 2008). In the Hele-Shaw flow, the fluid in the far field is driven by a background pressure gradient so that it is yielded across part of the gap thickness. However, the material becomes rigid across the entire thickness in stagnant ‘dead’ zones upstream and downstream of the obstruction.

For the flow set-up considered in the present paper (figure 1a), the far-field flow must be yielded near the inclined plane whilst fully rigid, ‘dead’ zones may appear both immediately upstream and downstream of the obstruction. The steady interaction of such flows with an isolated topographic mound was studied by Hinton & Hogg (2022). They showed that the yield stress of the fluid reduces the diversion of flux around the mound and also leads to a thicker accumulation of fluid upstream of the mound. In this study we investigate the flow around a surface-piercing cylindrical obstruction, which is substantially different to flow over a mound with smoothly varying topography. However, the downstream reconnection of the flow in the lee of tall and wide mounds and obstructions is similar.

There are some analogies with obstructed free-surface granular flows despite the importance of inertia in that context. Examples include the long dead zone downstream of wide obstructions and that in the stagnant dead zone upstream, the material is on the verge of flowing (Cui & Gray 2013; Tregaskis et al. 2022).

The paper is structured as follows. The model is described in § 2 with the flow governed by two dimensionless groups. Numerical results for the flow thickness and the dead zones are reported in § 3, with full details of the numerical implementation given in Appendix A. In § 4, we analyse flow around relatively wide obstructions. A flux balance argument is deployed to quantify the flow thickness and this identifies the role played by the yield stress. Some of the peripheral asymptotic details are given in Appendix B. The case of relatively narrow obstructions is then investigated in § 5, where we show that there is an equivalence with flow in a Hele-Shaw cell (Hewitt et al. 2016), a connection that allows the shape of the dead zone to be approximated. Concluding remarks are given in § 6. Throughout we focus on obstructions with square cross-section but give the key differences for a rhombus-shaped cross-section in Appendix C and circular cross-section in Appendix D.

2. Theoretical model

We analyse the steady free-surface flow of a Bingham fluid on an inclined plane at angle $\beta$ to the horizontal (figure 1a). We consider a Bingham fluid because it has the simplest constitutive law that exhibits a yield stress enabling us to identify the influence of this phenomenon. In Bingham's model, the fluid is rigid for stresses below the yield stress, $\hat {\tau }_Y$ and the fluid yields with (constant) shear viscosity $\hat {\mu }$ when the yield stress is exceeded (Bingham 1916; Balmforth & Craster 1999). The fluid has constant density $\hat {\rho }$. The coordinate axes are oriented as follows: the $\hat {x}$ axis is directed downslope, the $\hat {y}$ axis cross-slope and the $\hat {z}$ axis is perpendicular to the plane (with $\hat {z}=0$ on the plane surface). The steady flow thickness is denoted by $\hat {h}(\hat {x},\hat {y})$. The flow is obstructed by a rigid obstruction with square cross-section lying in $(\hat {x},\hat {y}) \in [-\hat {L},\hat {L}]\times [-\hat {L},\hat {L}]$, which is assumed to be sufficiently tall that the flow never surmounts the obstruction.

The flow is supplied by a line source, with infinite extent in the $\hat {y}$ direction, which is located far upstream of the obstruction and delivers a sustained flux $\hat {Q}_{\infty }$ per unit length. We assume that the flow is relatively thin, inertia-less and that surface tension is negligible. A lubrication approximation is applied and so no-slip is not imposed on the obstruction boundary; the effects of friction on the obstruction walls is confined to a small unimportant boundary layer that we neglect (see Balsa 1998). The boundary conditions are therefore no flux on the obstruction boundaries and that the flow thickness returns to its unperturbed constant value far from the obstruction, which we denote $\hat {H}_{\infty }$.

We non-dimensionalise in plane lengths with the obstruction half-width, $\hat {L}$, and flow thicknesses with the thickness of Newtonian flow with viscosity $\hat {\mu }$ (Nusselt 1916)

(2.1)\begin{equation} {\hat{H}_N = \left(\frac{3 \hat{\mu} \hat{Q}_{\infty}}{\hat{\rho} \hat{g} \sin \beta} \right)^{1/3}}. \end{equation}

This choice of thickness scale ensures that the dimensionless upstream flux is unity for any value of the yield stress. It is equivalent to scaling the flux with $\hat {Q}_{\infty }$. The dimensionless variables are

(2.2a,b)\begin{equation} {(x,y)=(\hat{x},\hat{y})/\hat{L}}, \quad (z,h)=(\hat{z},\hat{h})/\hat{H}_N. \end{equation}

The governing equation for the flow is derived by balancing in-plane hydrostatic pressure gradients with the divergence of shear stress. The stress vanishes at the free surface and so the material is always unyielded to leading order in a neighbourhood of the free surface (Balmforth & Craster 1999). If the stress at the substrate ($z=0$) exceeds the yield stress then the layer is partitioned into a yielded region below a pseudo-plugged region; otherwise the entire layer is rigid and stationary (see figure 1b). Then using the Bingham constitutive law, the velocity field is obtained. Integrating over the flow thickness enables the dimensionless volume flux parallel to the plane, $\boldsymbol {q}$, to be calculated in terms of the flow thickness and its gradient (see Balmforth, Craster & Sassi 2002)

(2.3)\begin{equation} \boldsymbol{q} = \frac{1}{2}Y^2(3h-Y)\left(1- \mathcal{L}^{-1} \frac{\partial h}{\partial x}, -\mathcal{L}^{-1} \frac{\partial h}{\partial y}\right), \end{equation}

where the height of the yield surface is

(2.4)\begin{equation} Y(x,y) = \max \left(0,h - \frac{B}{\sqrt{\left(1-\mathcal{L}^{-1} \dfrac{\partial h}{\partial x}\right)^2 +\left(\mathcal{L}^{-1} \dfrac{\partial h}{\partial y}\right)^2}}\right), \end{equation}

and the flow is governed by two dimensionless parameters

(2.5a,b)\begin{equation} {B = \frac{\hat{\tau}_Y}{\hat{\rho} \hat{g} \hat{H}_N \sin \beta}, \quad \mathcal{L}=\frac{\hat{L} \tan \beta}{\hat{H}_N}.} \end{equation}

The former is the Bingham number, which is the ratio of the yield stress to a characteristic gravitationally induced stress, whilst the latter is an aspect ratio, which is associated with the relative width of the obstruction. The terms in the large brackets in (2.3) can be written as $(1,0)-\mathcal {L}^{-1}\boldsymbol {\nabla } h$, which reveals the effect of hydrostatic pressure gradients associated with the inclined plane and the free-surface shape, respectively. The latter acts to smooth free-surface gradients. The flow is yielded in $0< z< Y$ and pseudo-plugged in $Y< z< h$. In regions where $Y=0$, the entire thickness of the layer is truly plugged (e.g. figure 1b). The steady flow satisfies

(2.6)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q} =0, \end{equation}

with $\boldsymbol {q} \boldsymbol {\cdot }\boldsymbol {n} = 0$ on the obstruction boundaries.

Under this formulation, the dimensionless flux per unit width is unity far upstream; $\boldsymbol {q} \to (1,0)$ as $x \to -\infty$. The dimensionless flow thickness satisfies $h\to h_{\infty }=\hat {H}_{\infty }/\hat {H}_N$ far from the obstruction with $h_{\infty }$ given by the solution to the following (far-field) flux balance

(2.7)\begin{equation} 1= (h_{\infty} - B)^2(h_{\infty} +B/2). \end{equation}

The relationship between $h_{\infty }$ and $B$ is shown in figure 2 as a black line. For all $B$, we have $h_{\infty }>\max (1,B)$, and for small $B$,

(2.8)\begin{equation} h_{\infty} = 1+ B/2 + \ldots, \end{equation}

whilst for large $B$,

(2.9)\begin{equation} h_{\infty} = B + \left(\frac{2}{3B} \right)^{1/2} + \ldots. \end{equation}

These predictions are included in figure 2 as red and blue dashed lines.

Figure 2. The dimensionless far-field flow thickness $h_\infty$ as a function of the Bingham number $B$ (black line). The predictions with $B \ll 1$ (red dashed line) and $B \gg 1$ (blue dashed line) are also plotted (from (2.8) and (2.9)).

3. Numerical results and general observations

The system consisting of (2.3), (2.4) and (2.6), with boundary conditions of no-flux on the obstruction boundaries and $h\to h_{\infty }$ as $x^2+y^2 \to \infty$ is integrated numerically to find the steady flow thickness $h(x,y)$. Full details of the numerical method are given in Appendix A. The flow thickness calculated from the numerical method for $B=1$ and $\mathcal {L}=1$ is shown in figure 1(c). The domain of the numerical computation extends far beyond the region plotted in the figures. The steady flow deepens upstream of the obstruction and thins downstream of the obstruction. Our numerical results confirmed that the streamwise extent of the obstruction has a negligible influence on the upstream and downstream flow. We focus on square obstructions noting that almost identical results apply to rectangular cross-sections with the same width.

Throughout this paper, we analyse how the flow direction and thickness is influenced by the obstruction, with particular focus on the dependencies on the relative obstruction width (via $\mathcal {L}$) and the relative magnitude of the yield stress (via $B$). We find that the flow behaviour is delineated into four asymptotic regimes, shown in parameter space in figure 3. These regimes are partitioned into the case that the obstruction is wide relative to the flow thickness ($\mathcal {L} \gg h_{\infty } \approx \max (1,B)$) and the case that the obstruction is relatively narrow, on the right-hand and left-hand side of parameter space, respectively. The ‘wide’ regime is analysed first in § 4 and the ‘narrow’ regime in § 5.

Figure 3. The four regimes of flow around a square obstruction in ($\mathcal {L}$, $B$)-parameter space. The colour maps represents the steady flow thickness, $h(x,y)$. The wide regimes are shaded grey.

However, before launching into the analysis for each of these regimes, we comment on two general features of the flow; ‘dead’ zones at the stagnation points (§ 3.1), and the potential for non-unique solutions (§ 3.2).

3.1. Dead zones at the stagnation points

The height of the yield surface in figure 1(d) illustrates that ‘dead’ zones, in which the material is entirely rigid and static, develop upstream and downstream of the obstruction. Although the yield stress is regularised somewhat in our numerical scheme (see Appendix A), the log scale enables the distinction of the dead zone to be clearly seen as black regions in figure 1(d) (here we have chosen the dead zone to correspond to $Y \lesssim 2 \times 10^{-4}$ though this arbitrary threshold has little influence on the interpretation of the flow structure).

For $B>0$, there is always a dead zone in a neighbourhood of the two stagnation points at $(\pm 1,0)$; see, e.g., figure 4. This is easily argued by contradiction. Suppose that the material is yielded at $z=0$ at $(\pm 1,0)$ with $Y>0$. Then, by symmetry, the transverse flux $\boldsymbol {q}\boldsymbol {\cdot } \boldsymbol {e}_y$ vanishes, and by the boundary condition on the wall, $\boldsymbol {q}\boldsymbol {\cdot } \boldsymbol {e}_x$ vanishes (where $\boldsymbol {e}_x$ and $\boldsymbol {e}_y$ are unit vectors in the $x$ and $y$ directions, respectively). With $Y>0$, (2.3) furnishes $\partial h/\partial y=0$ and $\partial h/\partial x= \mathcal {L}$. The yield stress term in (2.4) is then singular implying that $Y=0$, a contradiction. The material is always rigid at these locations because the flow stagnates there and there is no hydrostatic pressure gradient to induce any shear stress in the layer.

Figure 4. (a) Steady flow thickness, $h(x,y)$, and (b) yield surface height, $Y(x,y)$, for $B=1$ and $\mathcal {L}=4$ (calculated numerically); (c) $h(x,y)$ and (d) $Y(x,y)$ for $B=1$ and $\mathcal {L}=16$.

The size of the two dead zones increases with the dimensionless yield stress $B$. The upstream dead zone becomes smaller for wider obstructions (greater $\mathcal {L}$) whilst the downstream dead zone becomes larger (e.g. figure 4). This is because for greater $\mathcal {L}$, the magnitude of the ‘smoothing’ term, $-\mathcal {L}^{-1} \boldsymbol {\nabla } h$ in (2.3) is diminished. Hence, the cross-slope gradient of the free surface at the upstream boundary is larger in order to divert the flow (see § 4), and the flow reconnects downstream over a greater distance (see Hinton & Hogg 2022).

A simple corollary is that there can be a ‘cusp’ at the upstream stagnation point, i.e. $\partial h/ \partial y$ is discontinuous. This is a common feature of plugged viscoplastic material (Balmforth et al. 2002). The cusp is borne out by our asymptotic analysis for relatively wide obstructions. It is slightly smoothed in the numerical results owing to the regularisation of the yield stress.

3.2. Non-uniqueness of the steady flow

Throughout this paper, we primarily focus on the upstream flow, noting that for $\mathcal {L} \gg 1$ there is a long dead zone downstream of the obstruction and the flow will reconnect over a length scale proportional to $\mathcal {L}$ (Hinton & Hogg 2022). Details of the flow structure at the edge of the downstream dead zone and the possibility of a non-unique steady state were investigated in Hinton & Hogg (2022) for flow around large mounds corresponding to $\mathcal {L} \gg 1$. The steady state downstream will depend on the initial condition prior to the supply of constant flux from upstream. For example, the plane may have been entirely dry with $h=0$ everywhere in which case the downstream dead zone mostly corresponds to a dry zone with $h=0$ (provided that $\mathcal {L} \gg 1$). If instead, the plane was coated with a uniform, rigid and stationary layer of fluid with thickness $h=B$ prior to the line source initiation, then the downstream dead zone would mostly be a region in which $h=B$. In the numerical results presented in this paper, we have assumed that this is the case; see, for example, figure 4. It should also be noted that any $0\leqslant h \leqslant B$ is allowed in the downstream dead zone provided that the layer is entirely rigid. Importantly, the downstream behaviour has negligible effect on the upstream thickness.

4. Flow past a relatively wide obstruction ($\mathcal {L} \gg h_{\infty }$)

For a relatively wide obstruction ($\mathcal {L} \gg h_{\infty } \approx \max (1,B)$), there is significant flow deepening at the upstream boundary; see figures 3 and 4. In addition, there is a long dead zone downstream. This regime is consistent with lubrication theory provided that

(4.1)\begin{equation} \mathcal{L} \gg \max(h_\infty, h_\infty \tan \beta). \end{equation}

Upstream of a wide obstruction, the fluid forms a deep pond. The extent of the pond in the $x$ direction is much smaller than the width of the obstruction; e.g. figure 4(c). Hence, $\partial /\partial x \gg \partial /\partial y$ and since there is no flux into the obstruction, we have $\boldsymbol {q} \boldsymbol {\cdot } \boldsymbol {e}_x \approx 0$ over the extent of the pond, where $\boldsymbol {e}_x$ is the unit vector in the $x$ direction. This corresponds to a free surface that is horizontal in the streamwise direction (i.e. perpendicular to the direction of gravity), as has been observed previously for Newtonian flow past obstructions (Hinton et al. 2020). These observations motivate writing the flow thickness in the deep region as

(4.2)\begin{equation} h=(1+x)\mathcal{L} + G(y), \end{equation}

where $G(y) \gg 1$ is the flow thickness on the upstream boundary, $x=-1$, which is to be determined by equating the flux arriving into the deep pond from upstream with the transverse flux along the pond. In the pond, the yield surface is located at

(4.3)\begin{equation} Y=\max\left(0, (1+x) \mathcal{L}+G(y) - \frac{B}{\mathcal{L}^{-1} |G'(y)|}\right). \end{equation}

The volume flux arriving into the pond from upstream in $[0,y]$ is $y$ since the dimensionless flux far upstream is unity per unit length. The transverse flux at $y$ is calculated from (2.3) over the region in which $Y >0$. The steady flux balance takes the form (cf. Lister & Hinton 2022)

(4.4)\begin{equation} y=\int_{pond} \boldsymbol{q} \boldsymbol{\cdot} \boldsymbol{e}_y \, \mathrm{d}\kern0.7pt x= \int_{Y>0} - \textstyle{\frac{1}{2}} \mathcal{L}^{-1} Y^2 (3h-Y) G'(y)\, \mathrm{d}\kern0.7pt x. \end{equation}

We recast this as an integral over $Y$ using (4.2) and (4.3)

(4.5)\begin{equation} y=- \mathcal{L}^{-2} G'(y) \int_0^{Y_0} Y^2\left(Y+ \frac{3B}{2 \mathcal{L}^{-1} |G'(y)|}\right) \mathrm{d} Y. \end{equation}

Note that $Y$ is a linearly increasing function of $x$, and

(4.6)\begin{equation} Y_0=Y(x=-1)=\max\left(0,G(y) - \frac{B}{\mathcal{L}^{-1}|G'(y)|}\right) \end{equation}

is its maximum value at the wall $x=-1$. Carrying out the integration furnishes

(4.7)\begin{equation} y = \frac{-\mathcal{L}^{-2} G'(y)}{4} \left( G(y) - \frac{B}{\mathcal{L}^{-1} |G'(y)|}\right)^3 \left( G(y) + \frac{B}{\mathcal{L}^{-1} |G'(y)|}\right). \end{equation}

The left-hand side represents the flux from upstream and the right-hand side represents the transverse flux. Equation (4.7) is a nonlinear first-order ordinary differential equation (ODE) with boundary condition $G(1)=0$ corresponding to the flow thickness returning to lower order at the sides of the obstruction (e.g. figure 4c). The dominant balance in this equation differs depending on the relative magnitude of $\mathcal {L}$ and $B$.

For a Newtonian fluid ($B=0$), (4.7) gives the pond depth as $G \sim \mathcal {L}^{2/5}$ (Hinton et al. 2020). This suggests that the relative importance of the yield stress terms in (4.7) is given by the quantity $B \mathcal {L}^{1/5}$. These observations motivate defining

(4.8a,b)\begin{equation} \lambda = B \mathcal{L}^{1/5} \quad \text{and} \quad G(y) = \mathcal{L}^{2/5} \bar{G}(y). \end{equation}

Equation (4.7) becomes

(4.9)\begin{equation} y = \frac{-\bar{G}'}{4} \left( \bar{G} - \frac{\lambda}{ |\bar{G}'|}\right)^3 \left( \bar{G}+ \frac{\lambda}{|\bar{G}'|}\right). \end{equation}

The rescaled flow thickness at the wall, $\bar {G}(y)$, depends only on the single parameter, $\lambda$, which represents the relative importance of the yield stress. Equation (4.9) is integrated numerically using the boundary condition $\bar {G}(1-\delta )=(20 \delta )^{1/5}$, where $\delta = 10^{-8}$. At each step in the negative $y$ direction, $\bar {G}'$ is obtained by applying Newton's method to (4.9), with the restriction that the terms in brackets are positive.

The solution for $\bar {G}(y)=h(x=-1,y)/\mathcal {L}^{2/5}$ is shown as black lines in figure 5(a) for $\lambda =0.1$, $1$ and $10$. The magenta stars in figure 5(a) represent the thickness at the wall ($x=-1$) from the full numerical integration of the governing partial differential equation, with $\mathcal {L}=20$ and $B$ chosen to give the required value of $\lambda$. The red dashed and blue dot-dashed lines are the small and large yield stress asymptotic results obtained in §§ 4.1 and 4.2. Figure 5(b) shows the rescaled maximum thickness (i.e. $h(-1,0)$) as a function of $\lambda$ with the lines calculated as in figure 5(a).

Figure 5. (a) Rescaled flow thickness at the wall $x=-1$ as a function of cross-slope position. (b) Rescaled maximum flow thickness as a function of $\lambda$. The black lines are obtained by integrating (4.7). The results from the full numerical method with $\mathcal {L}=20$ are shown as magenta stars. The red dashed and blue dot-dashed lines are the small and large yield stress asymptotic predictions, respectively (see §§ 4.1 and 4.2).

The flow thickness increases for larger $\lambda$. A steeper free surface is required to divert the upstream flux at higher yield stresses because a greater proportion of the layer is rigid. The relative importance of the plugged material in the overall flow depends primarily on $B$ but also weakly on $\mathcal {L}$; wider square obstructions exacerbate the plug size because the term in the stress associated with the transverse free-surface gradient, $-\mathcal {L}^{-1} \partial h/\partial y$ is smaller. Note that this does not occur for obstructions with walls that are angled to the downslope direction because the free-surface gradient is not the dominant contribution in the yield stress term; see Appendix C. We analyse the low and large yield stress regimes in §§ 4.1 and 4.2, respectively.

4.1. Relatively low yield stress ($\lambda \ll 1$)

In terms of the two original dimensionless groups, $B$ and $\mathcal {L}$, the present $\lambda \ll 1$ regime is (see figure 3)

(4.10a,b)\begin{equation} {B \ll \mathcal{L}^{-1/5}, \quad \mathcal{L} \gg 1.} \end{equation}

This regime corresponds to quasi-Newtonian flow; see figure 6. We seek a regular asymptotic expansion for the flow thickness at the wall of the form

(4.11)\begin{equation} \bar{G}(y) = \bar{G}_0(y) + \lambda \bar{G}_1(y) + \dots. \end{equation}

The terms are obtained from the ODE (4.9). The first term is given by the Newtonian solution (corresponding to $B=0$, see Hinton et al. 2020)

(4.12)\begin{equation} \bar{G}_0(y) = 10^{1/5} (1-y^2)^{1/5}. \end{equation}

The next term accounts for the yield stress (for details, see Appendix B.1)

(4.13)\begin{equation} \bar{G}_1(y)=\frac{0.9549-p(y)}{(1-y^2)^{4/5}}, \quad \text{where} \ p(y) = \frac{10^{4/5}}{5} \int_0^{|y|} (1-s^2)^{3/5} \, \mathrm{d} s. \end{equation}

The asymptotic prediction for $\lambda \ll 1$ is compared with the numerical results in figures 5(a) and 6(c), where good agreement is observed. The height of the yield surface in the deep region (4.3) is given by

(4.14)\begin{equation} Y=\max\left(0, \mathcal{L} (x+1) +\mathcal{L}^{2/5} \bar{G}_0(y) +B \mathcal{L}^{3/5}\left(\bar{G}_1(y)- \frac{1}{ |\bar{G}_0'(y)|}\right) + \dots \right), \end{equation}

which is compared with the numerical result along $y=0.6$ in figure 6(d). In this regime, the flow thickness, $h$, and the height of the yield surface, $Y$, are of the same order of magnitude in the deep ponded region, and the lateral extent of the deep region is similar to the lateral extent of the flowing region (where $Y>0$); see figure 6. Equation (4.14) can be used to obtain the boundary of the region in which $Y=0$, up to and including second-order terms,

(4.15)\begin{equation} x_B(y) =-1+ \min \left(0, \mathcal{L}^{-3/5} \bar{G}_0(y) + B \mathcal{L}^{-2/5} \left(\bar{G}_1(y)- \frac{1}{ |\bar{G}_0'(y)|}\right) \right), \end{equation}

which is shown as a blue dashed line in figure 6(b). This boundary intersects the obstruction wall near the centreline where $\bar {G}_0'(y)$ becomes small, and a different asymptotic expansion is needed in $|y| \sim \lambda$. This inner region arises because the yield stress plays a more significant role near the centreline owing to the relatively smaller gradients of flow thickness there. Indeed, as argued above, there is always a dead zone at the stagnation point. Details of this region and the matching with the outer expansion are given in Appendix B.1. Importantly, the flow thickness at $y=0$ is unchanged at first and second order from the outer expansion (4.11). Hence the maximum flow thickness, which occurs at $x=-1$, $y=0$, is given by

(4.16)\begin{equation} h_{max} = 10^{1/5} \mathcal{L}^{2/5} + 0.9549 B \mathcal{L}^{3/5}+ \dots, \end{equation}

which is shown as a red dashed line in figure 5(b). The present analysis is based upon $\lambda \ll 1$ and $h \gg h_{\infty }$ so that the expansion (4.11) remains asymptotic. In the next subsection, we analyse the other regime, of large $\lambda$.

Figure 6. Results for a wide obstruction and relatively low yield stress; $\mathcal {L}=40$ and $B=0.1$. (a) Steady flow thickness and (b) yield surface from the full numerical method. The boundary of the flowing region predicted by (4.15) is shown as blue dashed lines. (c) Flow thickness along $y=0$ (blue) and $y=0.6$ (black). The red dashed lines show the asymptotic prediction (4.2) with $\bar {G}(y)$ given by (4.11). (d) Corresponding height of the yield surface; red dashed line is (4.14).

Figure 7. Results for a wide obstruction and relatively high yield stress; $\mathcal {L}=40$ and $B=2$. (a) Steady flow thickness and (b) yield surface from the full numerical method. (c) Flow thickness along $y=0$ and $y=0.6$. The red dashed lines show the asymptotic prediction (4.2) with $\bar {G}(y)$ given by (4.18). (d) Flow thickness at the wall comparing the large yield stress prediction (red dashed line) with the result from the ODE (4.7) (blue dot-dashed line) and the full numerical result (black line).

4.2. Relatively large yield stress ($\lambda \gg 1$)

In terms of the two original dimensionless groups, $B$ and $\mathcal {L}$, the present $\lambda \gg 1$ regime is (see figure 3)

(4.17)\begin{equation} {\mathcal{L}^{-1/5}\ll B \ll \mathcal{L}.} \end{equation}

In this regime, the role of the yield stress is non-negligible and appears in the leading order of the expansion for $\bar {G}$. The fluid is plugged to leading order in the deep region upstream of the obstruction; see figure 7 and note the qualitatively different contour structure to the low yield stress case in figure 6(a). Hence, the first term in brackets in (4.7) becomes small. This leading-order behaviour provides no transverse flux. This motivates the following expansion

(4.18)\begin{equation} {\bar{G} = \lambda^{1/2}\left( \tilde{G}_0(y) + \lambda^{-5/6} \tilde{G}_1(y) + \dots \right),} \end{equation}

where the second term provides the transverse flux to balance the incoming flow from upstream. The leading-order plugged shape is then calculated from (4.9)

(4.19)\begin{equation} {\tilde{G}_0(y) = \left[2(1-|y|) \right]^{1/2}.} \end{equation}

The next term, which accounts for the flux balance (and hence the left-hand side of (4.7)), is (for details, see Appendix B.2)

(4.20)\begin{equation} {\tilde{G}_1(y) = \frac{1.0600 - r(y)}{(1-|y|)^{1/2}}} \quad \text{where} \ r(y) = \int_0^{|y|} \frac{s^{1/3}}{2^{2/3} (1-s)^{1/2}} \, \mathrm{d} s. \end{equation}

This prediction is shown as red dashed lines in figures 7(c) and 7(d) with $B=2$ and $\mathcal {L}=40$; the numerical results are shown as continuous lines and the blue dot-dashed line is from integrating (4.7). The height of the yield surface in the deep region is found to be (from (4.3))

(4.21)\begin{equation} Y=\mathcal{L}(1+ x) +\mathcal{L}^{1/3} B^{-1/3} 2^{1/3} |y|^{1/3} + \dots. \end{equation}

As expected, the height of the yield surface ($\sim \mathcal {L}^{1/3} B^{-1/3}$) is asymptotically smaller than the height of the free surface ($\sim B^{1/2} \mathcal {L}^{1/2}$) in this regime; see the colourscales in figures 7(a) and 7(b). The extent of the deep region is $x \sim \mathcal {L}^{-1/2} B^{1/2}$, owing to the quasi-plugged term, $\tilde {G}_0$. The extent of the region over which there is appreciable transverse flux is asymptotically smaller (in contrast to the case of relatively small yield stress) and is given by

(4.22)\begin{equation} x_B(y) =-1 -\mathcal{L}^{-2/3} B^{-1/3} 2^{1/3} |y|^{1/3}. \end{equation}

The maximum flow thickness is

(4.23)\begin{equation} h_{max} = 2^{1/2} B^{1/2} \mathcal{L}^{1/2} + 1.0600 B^{-1/3} \mathcal{L}^{1/3}+ \dots, \end{equation}

which is shown as a blue dot-dashed line in figure 5(b). The second term in (4.23) is smaller than the first provided that $\lambda \gg 1$, as expected.

For large fixed $\mathcal {L}$, the present analysis will, in fact, break down for very large $B$. This is because the far-field flow thickness is $h_{\infty } \approx B$ for large $B$, and this exceeds the pond flow thickness, ($\sim B^{1/2} \mathcal {L}^{1/2}$) for $B \gg \mathcal {L}$, which invalidates the assumption of a deep pond. Hence the expansion (4.18) is valid when (4.17) applies. The case $B \gg \mathcal {L} \gg 1$ instead corresponds to a ‘deep’ oncoming flow relative to the obstruction width, which is analysed in § 5 (see also the parameter space in figure 3). For example, in figure 7(d), there is excellent agreement as $B=2$, $\mathcal {L}=40$ and $\lambda = 4.183$, whereas the agreement in figure 5(a) for $\lambda =10$ is not so good, especially near the edges of the upstream wall, because $\mathcal {L}=20$ and $B=5.493$ with a far-field flow thickness of $h_{\infty }=5.834$.

5. Flow past a relatively narrow obstruction ($\mathcal {L} \ll h_{\infty }$)

In this section, we analyse the regime in which

(5.1)\begin{equation} {\mathcal{L} \ll h_{\infty} \approx \max(1,B)} \end{equation}

corresponding to flow past a relatively narrow square obstruction. This regime is consistent with lubrication theory provided that

(5.2)\begin{equation} h_\infty \tan \beta \ll \mathcal{L} \ll h_\infty. \end{equation}

This requires a shallow inclination angle; $\beta \ll 1$. Numerical results for the flow thickness and yield surface height are shown in figure 8 for $\mathcal {L}=0.1$ and three values of $B$. The perturbation to the flow thickness is approximately antisymmetric about $x=0$ whilst the yield surface and dead zones are symmetric about this line. For a low yield stress, the dead zone is localised to the stagnation points. At larger yield stresses, the dead zone encompasses the entire upstream and downstream boundaries and becomes progressively longer up and down slope. The perturbation to the flow thickness reduces at larger yield stresses. In what follows, we analyse the flow in the narrow obstruction regime and draw analogy with flows of viscoplastic fluids in a Hele-Shaw cell in order to understand these observations (see § 1 and Hewitt et al. 2016).

Figure 8. (a) Steady flow thickness and (b) the associated height of the yield surface with the analytic dead zone prediction (5.9), shown as blue lines, for flow around a square obstruction with $\mathcal {L}=0.1$ and $B=0.5$. (c) Flow thickness and (d) yield surface height for $\mathcal {L}=0.1$ and $B=2$. (e) Flow thickness and (f) yield surface height for $\mathcal {L}=0.1$ and $B=8$ (a larger domain is shown to capture the extent of the dead zones).

5.1. Analogy with flow in a Hele-Shaw cell ($\mathcal {L} \ll h_{\infty }$)

To dissect the flow physics in the narrow obstruction regime, we seek an asymptotic expansion of the form

(5.3)\begin{equation} h=h_{\infty} + \mathcal{L} \left[1+x+p(x,y) \right] + O\left( \mathcal{L}^2 \right), \end{equation}

where $p(x,y)$ is associated with the hydrostatic pressure that drives the flow. To leading order in $\mathcal {L}$, the flux (2.3) is

(5.4)\begin{equation} \boldsymbol{q} =- \frac{1}{2} Y^2 \left(3h_{\infty}-Y \right) \boldsymbol{\nabla} p, \quad \text{where} \ Y=\max\left(0, h_{\infty}-\frac{B}{|\boldsymbol{\nabla} p|} \right). \end{equation}

This form of the flux is identical to that for viscoplastic flow in a Hele-Shaw cell of fixed gap width, $h_{\infty }$, with $p(x,y)$ representing the pressure that drives the flow with constant flux (which is supplied far upstream) (Hewitt et al. 2016). Hence, to leading order, the streamlines and the pressure gradient in free-surface flow around a relatively narrow obstruction are given by those in a Hele-Shaw cell with a blockage of the same cross-section. The similarity arises because in this narrow obstruction regime, the thickness of the free-surface flow is very minorly perturbed and the flux is approximately independent of this thickness. The leading-order flux depends only on the gradients of the free surface (see (5.4)) and so the system behaves in an analogous fashion to flow in a Hele-Shaw cell of fixed gap width. This equivalence was previously observed for viscous Newtonian flows (Hinton et al. 2020).

These observations conveniently allow us to exploit previous research on obstructed viscoplastic flow in a Hele-Shaw cell (Hewitt et al. 2016; Daneshi et al. 2020). For example, the symmetry of the two dead zones (figure 8) is rationalised by the reversibility of steady flow in a Hele-Shaw cell. There is a slight difference in that here the far-field flow thickness (corresponding to the gap width) varies with $B$ but this merely leads to some extra factors of $h_\infty$ in the results (this alters some of the scalings, particularly for large $B$).

To analyse the flow, we use a hodograph transformation. The important results are quoted in the following, with full details given by Neuber (1961), Entov (1970) and Hewitt et al. (2016). We transform $(x,y)$ space to the hodograph plane of $(Q, \theta )$, where $Q=|\boldsymbol {q}|$ is the magnitude of the flux and $\theta$ is the angle that the streamlines make with the downslope direction. In this way using the results from Hewitt et al. (2016) and correcting a sign error, we find

(5.5)\begin{equation} \mathrm{d}\kern0.7pt x + \mathrm{i}\, \mathrm{d}y=-{\rm e}^{\mathrm{i} \theta} \left(\frac{\mathrm{d} p}{S} + \mathrm{i} S^{-2}\frac{\partial S}{\partial Q} \frac{\partial p}{\partial \theta} \,\mathrm{d} Q -\mathrm{i} \frac{Q}{S} \frac{\partial p}{\partial Q}\, \mathrm{d} \theta \right) \quad \text{where} \ S=|\boldsymbol{\nabla} p|. \end{equation}

The magnitude of (5.4)a relates $Q$ and $S$ via

(5.6)\begin{equation} Q=\left(h_{\infty}-\frac{B}{S} \right)^2 \left(h_{\infty} + \frac{B}{2 S} \right) S. \end{equation}

Local mass conservation, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {q}=0$, is recast as a linear elliptic partial differential equation for $p$,

(5.7)\begin{equation} \frac{S^2}{Q({\partial S}/{\partial Q})}\frac{\partial }{\partial Q} \left( \frac{Q^2}{S} \frac{\partial p}{\partial Q} \right) + \frac{\partial^2 p}{\partial \theta^2} =0, \end{equation}

where $S=S(Q)$ (see (5.6)). Some separable solutions to (5.7) exist, which enables analytical progress in § 5.2. The flow behaviour is qualitatively different in the regimes of a relatively low and relatively high yield stress, which are analysed in §§ 5.2–5.3, respectively.

5.2. Relatively low yield stress ($B \ll 1, \mathcal {L} \ll 1$)

For $0< B \ll 1$, $h_\infty \approx 1$ and so (5.1) becomes $\mathcal {L} \ll 1$; the present regime occupies the bottom-left of parameter space in figure 3. In this regime, viscoplastic stagnation point flow arises; see figure 9(a). A key feature is the appearance of a dead zone in the vicinity of the stagnation point; see figures 8(a) and 8(b).

Figure 9. Direction and magnitude of the flux $\boldsymbol {q}$ obtained from the numerical method with $\mathcal {L}=0.1$ and (a) $B=0.5$ and (b) $B=8$.

The shape of the dead zone can be determined analytically by noting that $Q=0$ on the boundary of the dead zone and there is stagnation point flow further away (with larger values of $Q$) with

(5.8)\begin{equation} p \sim (1+x)^2 -y^2 = \frac{Q^2}{4} \cos 2 \theta. \end{equation}

Note that this stagnation flow solution for $p$ is defined up to an arbitrary multiplicative and additive constant (and, hence, the same applies to (5.10)). Equation (5.8) provides the boundary condition for the partial differential equation in the hodograph plane (5.7). A separable solution is sought for (5.7) and the behaviour for $Q \ll 1$ is obtained by shooting numerically. The leading-order shape of the dead zone boundary and the pressure at the boundary are then parameterised by (see Hewitt et al. 2016)

(5.9)$$\begin{gather} (x+1,y) = 0.242 B h_{\infty}^2 (-\cos^3 \theta, \pm \sin^3 \theta), \end{gather}$$

(5.10)$$\begin{gather}p \sim 0.181 B^2 h_{\infty} \cos 2 \theta, \end{gather}$$

where $\theta \in [0, {\rm \pi}/2]$. The dead zone boundary (5.9) is compared with the numerical result for $B=0.5$ in figure 8(b). It extends a distance $0.242 B h_{\infty }^2$ away from the origin along both the $x$ and $y$ axes. An identical dead zone occurs at the downstream boundary.

The dead-zone shape (5.9) is valid provided that there is stagnation point flow at the upstream wall. This requires that the dead zone occupies only a small fraction of the wall. That is, $B h_{\infty }^2 \ll 1$ or, equivalently, $B \ll 1$ since $h_{\infty } \sim 1$ for small $B$. For larger values of $B$, the dead zone encompasses the corners of the square obstruction invalidating the assumption of stagnation point flow; see figures 8 and 9. Different analysis is needed in this case; see § 5.3.

For $0< B\ll 1$, the presence of the dead zone only alters the flow thickness at order $B^2$ (see (5.10)). Away from the dead zone, there is an outer flow that is quasi-Newtonian for which the yield stress influences the flow thickness at order $B$. In this outer flow, we seek the following expansion for the hydrostatic pressure

(5.11)\begin{equation} p(x,y) = p_N(x,y) + B p_1(x,y) + O(B^2), \end{equation}

where $p_N(x,y)$ is Newtonian solution. Using (5.3) and (5.4), we find that the pressure terms satisfy

(5.12a,b)\begin{equation} \nabla^2 p_N =0, \quad \nabla^2 p_1 = \frac{3}{2} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\boldsymbol{\nabla} p_N}{|\boldsymbol{\nabla} p_N|}\right), \end{equation}

with boundary conditions,

(5.13)$$\begin{gather} \frac{\partial p_N}{\partial n} =0, \quad \frac{\partial p_1}{\partial n} =0, \quad \text{for}\ (x,y) \in \partial \varOmega, \end{gather}$$

(5.14)$$\begin{gather}p_N \to -x-1, \quad p_1 \to 0, \quad \text{as}\ x^2+y^2 \to \infty, \end{gather}$$

where $\partial \varOmega$ is the boundary of the square and $\partial /\partial n$ denotes the normal derivative. The linear system for $p_N(x,y)$ and $p_1(x,y)$ is integrated numerically using the same finite-difference method as described in Appendix A. The solutions are shown in figures 10(a) and 10(b).

Figure 10. Steady flow around a relatively narrow obstruction with $B \ll 1$. (a) Leading-order pressure, $p_N(x,y)$ for Newtonian flow (see (5.12a,b)). (b) The pressure perturbation $p_1(x,y)$ associated with $B \ll 1$. Flow thickness for $\mathcal {L}=0.05$ (c) along the centreline, $y=0$, upstream of the obstruction, and (d) at the wall, $x=-1$ showing numerical results (black lines) and prediction (5.11) (dot-dashed lines).

The prediction for the flow thickness from (5.11) is compared with the full numerical results along the centreline ($y=0$) upstream of the obstruction, and at the wall ($x=-1$) in figures 10(c) and 10(d), respectively. The form of $p_1(x,y)$ predicted by the full numerical result for $\mathcal {L}=0.025$ and $B=0.5$ is shown in figure 11(a), whilst the difference between the numerical and asymptotic thickness is shown in figure 11(b); the error is at most $3.5 \times 10^{-3}$.

Figure 11. Comparison between the full numerical results and the asymptotic expansion given by (5.3) and (5.11) with $\mathcal {L}=0.025$ and $B=0.5$. (a) Pressure perturbation $p_1(x,y)$ inferred from the full numerical results (cf. figure 10b). (b) Difference in predicted flow thickness; $\mathrm {error}=h_{numerical} - h_{asymptotic}$.

As the stagnation point, $(-1,0)$, is approached the functions take the following limits

(5.15a,b)\begin{equation} p_N \to 1.353, \quad p_1 \to -0.679. \end{equation}

Although both $p_N(x,y)$ and $p_1(x,y)$ are finite everywhere, the present asymptotic analysis breaks down near the stagnation point because $|\boldsymbol {\nabla } p_N|$ vanishes there. The rigid region encompasses the entire thickness of the layer (see (5.4)) leading to the dead zone analysed previously. However, the pressure (and, hence, the flow thickness) is accurately approximated by $p=p_N(x,y) + B p_1(x,y)$ everywhere because the effect of the dead zone only appears at $O(B^2)$ in the pressure.

The maximum flow thickness is

(5.16)\begin{equation} h_{max}=h_{\infty}(B) + \mathcal{L}[1.353 - 0.679 B] + O(\mathcal{L}^2,\mathcal{L}B^2). \end{equation}

This is compared favourably to the numerical results in figure 12. The perturbation to the flow thickness reduces with increasing yield stress (this also occurs for the large $B$ regime; see § 5.3 and figure 13a). This contrasts with relatively wide obstructions for which the flow thickness perturbation increased with greater yield stress. For wide obstructions, the flow is primarily diverted around the obstruction by gradients in the free surface inducing transverse flux and for a fixed flux, greater gradients are needed at higher yield stresses.

Figure 12. The maximum flow thickness, $h_{max}$, as a function of $\mathcal {L}$ for $B=0$, $B=0.25$ and $B=0.5$. The lines show the small-$B$ prediction (5.16), and the plus signs are from the full numerical method; see Appendix A.

Figure 13. Maximum perturbation to the flow thickness, $h_{max}-h_{\infty }$ for $\mathcal {L}=0.1$. Analytic scaling (black line, $h_{max}-h_{\infty } \sim B^{-3/4}$ (5.25)) is compared with numerical results (blue crosses). (b) Flow thickness along the centreline for $\mathcal {L}=0.1$ and $B=8$. Numerical result (black line; see also figure 8e) is compared with the analytic prediction for the gradient (red dashed line; (5.24)).

The decrease in the flow thickness perturbation with $B$ for narrow obstructions can be interpreted through the relationship (5.6) between the volume flux, $Q$ and the pressure gradient $S=|\boldsymbol {\nabla } p|$. For $B \ll 1$, (5.6) becomes

(5.17)\begin{equation} Q = S + \tfrac{3}{2} B (S-1) + O(B^2). \end{equation}

Hence, the volume flux (relative to the far-field volume flux) of a Bingham fluid varies more rapidly with the pressure gradient, $S$, than a Newtonian fluid (this is because the pressure gradient alters the proportion of the layer that is yielded). Smaller variations in pressure are required to divert the volume flux around the obstruction for a Bingham fluid and so there is a smaller variation in the flow thickness.

5.3. Relatively high yield stress ($B \gg 1$, $\mathcal {L} \ll B$)

For large yield stresses, $h_\infty \approx B$ and (5.1) becomes $\mathcal {L} \ll B$; this case is depicted in the top-left of parameter space in figure 3. In this regime, the dead zone encompasses the entire upstream and downstream boundaries of the obstruction (see figure 8f). The streamwise extent of the dead zone increases with yield stress. Although it is not possible to use the stagnation point flow to obtain an analytic solution for the shape of the dead zone (see figure 9b), we can find scalings for the extent of the dead zone and the flow thickness. Fluid from upstream in $-1< y<1$ must be diverted around the obstruction via boundary layers at the edge of the dead zone in which $Q \sim 1$ (Hewitt et al. 2016). Since $h_{\infty } \sim B \gg 1$, (5.6) may be written as

(5.18)\begin{equation} S(Q) = \frac{B}{h_{\infty}} + \frac{1}{h_{\infty}^2} \left( \frac{2 B Q}{3} \right)^{1/2} + \dots. \end{equation}

Equation (5.7) then becomes, to leading order,

(5.19)\begin{equation} 6^{1/2} B^{1/2} h_{\infty} Q^{-1/2} \frac{\partial }{\partial Q} \left(Q^2 \frac{\partial p}{\partial Q} \right) + \frac{\partial^2 p}{\partial \theta^2} =0. \end{equation}

Since the flux from upstream is diverted by the flow in a neighbourhood of the dead zone, we have $Q=O(1)$ and balancing the two terms in (5.19) furnishes a scaling for the angle of the dead zone boundary (cf. Hewitt et al. 2016),

(5.20)\begin{equation} \theta \sim B^{-1/4} h_{\infty}^{-1/2}. \end{equation}

Since the edge of the dead zone meets the obstruction at $y=1$ and $\theta \sim y/x$, the dead zone extent is given by

(5.21)\begin{equation} x \sim B^{1/4} h_{\infty}^{1/2} \sim B^{3/4}. \end{equation}

These scalings are different to those in a Hele-Shaw cell (Hewitt et al. 2016) because the effective gap width, $h_{\infty }$ varies with $B$. Within the dead zone, $Y=0$ and using both (2.4) and (5.21) furnishes the scalings $x=B^{3/4} \chi$ and $h=h_{\infty } + \mathcal {L} B^{-3/4} h_1$. Upon substituting these into $Y=0$, we find that

(5.22)\begin{equation} h_{\infty} - \frac{B}{\sqrt{1- 2 B^{-3/2} \dfrac{\partial h_1}{\partial \chi} + B^{-3/2} \left(\dfrac{\partial h_1}{\partial y}\right)^2}} + O\left(\mathcal{L} B^{-3/4}, B^{-2}\right)=0. \end{equation}

Recalling that $h_{\infty }=B+[2/(3B)]^{1/2}+\dots,$ the thickness in the dead zone satisfies

(5.23)\begin{equation} \frac{\partial h_1}{\partial \chi} - \frac{1}{2} \left(\frac{\partial h_1}{\partial y} \right)^2 = \left(\frac{2}{3} \right)^{1/2}, \end{equation}

to leading order. With analytical boundary data at the edge of the dead zone, this equation could be solved via Charpit's method. Unfortunately, such data are not available (see Hewitt et al. 2016). However, we can make some progress by noting that along the centreline $\partial h_1/\partial y=0$ (there is no cusp to leading order in the $\mathcal {L} \ll 1$ regime) and, hence,

(5.24)\begin{equation} \frac{\partial h}{ \partial x}=B^{-3/2} \mathcal{L} \sqrt{2/3}, \end{equation}

which shows good agreement with the full numerical result (see figure 13b). The perturbation to the maximum flow thickness has the scaling

(5.25)\begin{equation} h_{max}-h_{\infty} \sim \mathcal{L} B^{-3/4}, \end{equation}

which is compared favourably to the numerical results in figure 13(a).

6. Conclusion

We have analysed steady free-surface viscoplastic flow around an obstruction of square cross-section on an inclined plane. The flow is characterised by two dimensionless parameters that measure the relative width of the obstruction, $\mathcal {L}$, and the relative magnitude of the yield stress, $B$. We have focused on the flow upstream of the obstruction and the dependence of the flow thickness upon the dimensionless parameters. Example parameter values for lava flow past a tree and a constructed barrier for two possible values of the yield stress are listed in table 1. The four possibilities span the four asymptotic regimes obtained; see figure 3 (although the tree, in particular, might be better modelled as a circular obstruction, as outlined in Appendix D).

Table 1. Dimensionless groups corresponding to lava flow past a tree or a constructed barrier for two values of the yield stress, $\hat {\tau }_Y$. We use $\hat {\rho } =2700\ \mathrm {kg}\ \mathrm {m}^{-3}$, $\hat {H}_N=6$ m and $\beta =15^\circ$. Parameter values from Fink & Griffiths (1998), Barberi et al. (2003) and Chevrel et al. (2019).

If the obstruction is relatively wide, then the fluid forms a deep pond upstream, which provides sufficient transverse flux to divert the oncoming flow. This deep region may be mostly plugged or mostly yielded depending on the magnitude of the yield stress relative to the viscous stresses and the obstruction width relative to the flow thickness and the inclination angle. These dependencies are altered for obstructions with different cross-section (see Appendix C and D). Free-surface flow around a relatively narrow obstruction is controlled by the same physics as pressure-driven flow in a Hele-Shaw cell, and so the viscoplastic free-surface flow is analogous to that in a Hele-Shaw cell with the same obstruction cross-section. The streamlines and pressure perturbation are identical in the two settings, to leading order, with fore–aft symmetry in the stagnant dead zones. Moreover, this equivalence holds for any obstruction cross-section and any generalised Newtonian fluid.

Our general results could be extended to more complicated rheologies such as Herschel–Bulkley fluids, noting that the qualitative flow structure and dependencies will be predominantly controlled by the yield stress and, hence, will be similar to the present analysis. Finally, it would be interesting to study the environmentally important case of a time-varying source of fluid and in particular the evolution when the source flux vanishes. A key question in this scenario is what stagnant fluid is left behind by the unsteady flow down the plane.