3.1. General observations

Figure 3(a) shows snapshots of a typical viscoplastic fluid injection flow experiment. Upon injection, the heavier viscoplastic fluid forms a jet with considerable inertial effects. The heavy viscoplastic fluid then slumps beneath the light ambient Newtonian fluid, initially maintaining a near-constant velocity. Beyond a certain point, the interfacial front velocity decreases, causing the viscoplastic fluid to accumulate near the inlet as the injection continues. Meanwhile, the interfacial front motion continuously slows down and eventually stops, with the interface becoming stationary. This corresponds to a complete interfacial arrest. Under continuous injection, once the interfacial front is arrested, the injected viscoplastic fluid redirects towards the outlet. Each experiment spans $O(10^2)$ seconds to confirm flow cessation and stabilisation.

Figure 3(b) presents the spatiotemporal representation of the normalised depth-averaged concentration field for the viscoplastic fluid in the experiment corresponding to figure 3(a). The colourmap illustrates the mean concentration, $\overline {\overline {c}}(z,t)$ , averaged over the pipe diameter, as a function of the spatial position, $z$ , and time, $t$ . Here,

(3.1) \begin{align} \overline {c}\left ( {y,z,t} \right ) &= \int _{ -1/2}^{ +1/2} {c\left ( {x,y,z,t} \right )} \,\mathrm{d}x, \end{align}

(3.2) \begin{align} \overline {\overline {c}}\left ( {z,t} \right ) &= \int _0^1 {\overline c\left ( {y,z,t} \right )}\, \mathrm{d}y, \end{align}

where the integration limits $-1/2$ and $+1/2$ correspond to the pipe’s half-width. The function $c(x, y, z, t)$ symbolises the true concentration field, while ${\overline {c}}(y, z, t)$ refers to the concentration captured in side-view experimental images. In general, the spatiotemporal map in figure 3(b) provides insights into the time-dependent flow dynamics, interfacial behaviours and front dynamics. For example, the inverse of the slope of the boundary between the coloured and white regions (identified by the black dashed lines) represents the dimensionless interfacial front velocity, $V_f$ . As time progresses and the interface develops, two interfacial flow phases emerge: at short times, inertia from the imposed flow dominates, causing the interfacial front to advance at a velocity nearly equal to the mean imposed flow velocity ( $V_f \approx 1$ in dimensionless form), defining the inertial-dominated phase. Later, the flow inertia diminishes, and the balance between viscoplasticity (encompassing effective viscosity and yield stress) and longitudinal buoyancy governs the forward motion of the interface, marking the viscoplastic-dominated phase. The solid black line in figure 3(b) and its inset marks the initial inertial-dominated phase ( $z \lesssim 2$ , $t \lesssim 2$ ), which transitions to a viscoplastic-dominated phase as inertial stresses diminish. In the latter phase, buoyancy drives the interface, but the front velocity slows down while the interface shape changes, indicating effective viscous stresses balancing buoyant stresses and causing the deceleration, i.e. a signature of viscous exchange flows (Séon et al. Reference Séon, Znaien, Salin, Hulin, Hinch, E.J. and & Perrin, B2007; Taghavi et al. Reference Taghavi, Alba, Séon, Wielage Burchard, Martinez and Frigaard2012b ). As indicated by the black arrow, the depth-averaged concentration behind the front increases as $V_f$ drops below the injection velocity, resulting in viscoplastic fluid accumulation. The vertical dashed lines at longer times exhibit the interfacial front arrest and stationary condition ( $V_f = 0$ ), suggesting a complete overcoming of all stresses by the yield stress. The unaffected region shown in white is indicative of a pure light fluid beyond the arrested interfacial front. The inset in figure 3(b) provides a detailed view of the inertial-dominated phase and the transition to the viscoplastic-dominated phase, highlighting the critical transition length ( $\mathcal{L}$ ) and time ( $\mathcal{T}$ ).

Figure 4(a) compares the time-dependent interfacial front position from the experiment (lines) with CFD simulation results (symbols) for three different cases. The CFD results exhibit reasonable agreement with the experimental results, validating the observed experimental phases, including the inertial- and viscoplastic-dominated phases. During the initial phase ( $z \lesssim 2$ ), i.e. the inertial-dominated phase, both CFD simulations and experiments show similar front velocities ( $V_f \approx 1$ ), as the flow is dominated by the imposed flow inertia. This observation is consistent with previous studies on buoyant miscible injections of Newtonian (Akbari & Taghavi Reference Akbari and Taghavi2021) and viscoplastic fluids (Faramarzi et al. Reference Faramarzi, Akbari and Taghavi2024; Akbari & Taghavi Reference Akbari and Taghavi2022a ), across various flow regimes that can generally occur in such injection flows, including mixing, separation, and stable and unstable slumping (Faramarzi et al. Reference Faramarzi, Akbari and Taghavi2024). Figure 4(b) provides simulation snapshots of the concentration profile, interface and front development, which align well with experimental data, albeit with some discrepancy in interface height, likely due to uncertainties in experimental conditions. While the experiment provides side-view concentration data, the CFD simulations reveal cross-sectional variations, as shown in figure 4(c), with contours of the heavy viscoplastic fluid. These results show a horizontal interface where the dense viscoplastic fluid slumps stably beneath the lighter fluid and progressively forms an unyielded zone as stress falls below the yield threshold. Note that although we focus on stable flows with negligible secondary flows, at ranges of the Péclet number considered ( $Pe = 10^5$ – $10^7$ ), thin mixing layers naturally develop at the interface, as observed experimentally (figure 3 a) and consistently captured in CFD simulations (figure 4 b).

Before proceeding, note that, from a CFD perspective, the analysis of concentration and velocity profiles (omitted for brevity) reveals that, despite the no-slip condition enforcing zero velocity at the wall, at longer times, the concentration near the ‘triple contact point’ approaches that of the pure viscoplastic fluid. The near-wall concentration becomes non-zero immediately after the viscoplastic front arrival, indicating no persistent Newtonian fluid computational ‘cell’ along the wall. This behaviour results from: (i) density-driven instability, where the initially unstable configuration – heavier viscoplastic fluid overlying the lighter Newtonian phase near the wall – rapidly dominates; and (ii) fluid miscibility, where diffusion and interfacial mixing, enhanced by the mesh and cell size near the wall, prevent the persistence of an unmixed Newtonian layer. Consequently, no permanent Newtonian pocket remains at the lower wall at long times.

Figure 4. CFD and experimental results. (a) Time-dependent interfacial front position from CFD simulations (symbols) versus experiments (lines) in the pipe’s centre plane ( $z$ – $y$ plane at $x=0$ ): ( , ) for $Re=1.3$ , $M=0.001$ , $B=0.74$ , $\beta =88^\circ$ , $Fr=0.55$ ; ( , ) for $Re=1$ , $M=0.0009$ , $B=0.75$ , $\beta =82^\circ$ , $Fr=0.47$ ; ( , ) for $Re=1.8$ , $M=0.001$ , $B=0.56$ , $\beta =82^\circ$ , $Fr=0.73$ . (b) CFD snapshots of concentration profiles, at the pipe centre plane (i.e. $z$ – $y$ plane at $x=0$ ), for $Re = 1.8$ , $M = 0.001$ , $B = 0.56$ , $\beta = 82^\circ$ and $Fr = 0.73$ . (c) Cross-sectional concentration contours in the $x$ – $y$ plane at $z = 6$ at different times, marked by arrows in panel (b). Colourbar shows heavy viscoplastic fluid concentration.

Figure 5. Experimental results. (a) Variation of interfacial front distance from the gate valve over time in log-log coordinates. Symbols represent experiments with: $\circ$ ( $Re = 3.4$ , $M = 0.0016$ , $B = 0.67$ , $\beta = 82^\circ$ , $Fr = 0.91$ ), $\square$ ( $Re = 0.24$ , $M = 0.0005$ , $B = 0.83$ , $\beta = 75^\circ$ , $Fr = 0.22$ ), $\lozenge$ ( $Re = 0.1$ , $M = 0.0003$ , $B = 0.75$ , $\beta = 82^\circ$ , $Fr = 0.14$ ) and $\triangle$ ( $Re = 1.8$ , $M = 0.001$ , $B = 0.56$ , $\beta = 82^\circ$ , $Fr = 0.73$ ). Red lines fit early-time data with slope ${\sim} 1$ in the inertial-dominated phase and brown lines fit the initial stage of the viscoplastic-dominated phase. (b) Dimensional interfacial front velocity in the inertial-dominant phase ( $\hat {V}_{fi}$ ) versus $\hat {V}_0$ ; dashed lines indicate 1 : 1 relationship. (c) Critical length, $\mathcal{L}$ (transition from an inertial-dominated phase to a viscoplastic-dominated phase) versus $Fr^2$ . (d) Power-law trend length $\mathcal{L}^\dagger$ versus $Re_b$ , with the dashed line marking a linear fit. (e) Sharp decay length ( $\mathcal{L}^{\dagger \dagger }$ ) versus ${\chi \cos \beta }/{B}$ , with the dashed line marking a linear fit. Marker sizes indicate $\cos \beta$ in panel (b–e); colours in panel (b) show $\hat {\mu }_H$ magnitude, while in the rest, they represent $Re$ . Symbols indicate $\Delta \hat {ho }$ : circles ( $1 \ \text{kg}\,{\textrm m^{-3}}$ ); squares ( $4 \ \text{kg}\,{\textrm m^{-3}}$ ); triangles ( $5 \ \text{kg}\,{\textrm m^{-3}}$ ); diamonds ( $10 \ \text{kg}\,{\textrm m^{-3}}$ ); pentagrams ( $20 \ \text{kg}\,{\textrm m^{-3}}$ ).

3.2. Flow, interface and front dynamics

This section analyses the flow, interface and front dynamics, highlighting distinct flow phases and regimes during the viscoplastic fluid injection. Figure 5(a) presents the temporal evolution of the interfacial front position, illustrating the inertial-dominated, viscous–buoyant and yield-stress-dominated phases. The latter two regimes can be collectively described as the viscoplastic-dominated phase, during which the effects of viscous and yield stresses increasingly control the front dynamics. In the initial inertial-dominated phase, the front advances at a constant velocity that matches the injection speed, resulting in a linear trend with a slope of approximately one in dimensionless coordinates. This behaviour is consistently observed across both CFD simulations and experiments, where the early-time data collapse in the dimensionless $t$ – $z$ space (figure 4 a). The flow then transitions into a viscous–buoyant regime, where the front velocity decays following a power-law trend with exponents ranging from 0.76 to 0.9. This regime reflects a balance between longitudinal buoyant stresses and viscous resistance and appears largely independent of the yield stress. The power-law behaviour manifests as a straight-line segment on a log–log plot, as shown by the brown lines in figure 5(a). Finally, the front enters a yield-stress-dominated phase, where the velocity decays more rapidly, deviating from the earlier power-law behaviour. This regime ultimately leads to complete interfacial arrest. The increasing effective viscous stress and diminishing longitudinal buoyant stress (due to interface flattening) result in progressive slowing until the yield stress fully dominates the flow.

Figure 5(b) presents the constant interfacial front velocity (dimensional) in the inertial-dominated phase, $\hat {V}_{fi}$ , plotted against the mean imposed flow, $\hat {V}_0$ , confirming that $\hat {V}_{fi}$ closely matches $\hat {V}_0$ . Minor discrepancies arise from measurement errors in $\hat {V}_0$ and $\hat {V}_{fi}$ . Our results identify a critical development length from the injection origin ( $\mathcal{L}$ ), i.e. the effective jet length, marking the transition from an inertial-dominated phase to a viscoplastic-dominated phase. Our findings suggest that $\mathcal{L}$ may be governed by a dimensionless group balancing the imposed inertial stress with the resisting buoyant stress, expressed as

(3.3) \begin{equation} \hat {\overline {ho }} \hat {V}_0^2 \sim \Delta \hat {ho } \hat {g} \hat {D} \quad ightarrow \quad \frac {\hat {\overline {ho }} \hat {V}_0^2}{\Delta \hat {ho } \hat {g} \hat {D}} \propto Fr^2, \end{equation}

where $\hat {\overline {ho }}={ ( {{{\hat ho }_H} + {{\hat ho }_L}} )}/{2}$ . Equation (3.3) shows consistency with the results presented in figure 5(c), demonstrating that $\mathcal{L}$ increases with increasing imposed inertial stress (i.e. increasing $Re$ , as indicated by the colourbar) or decreasing buoyant stress. A correlation can accordingly be suggested for $\mathcal{L}$ as

(3.4) \begin{equation} \mathcal{L}\approx 1.26 Fr^{0.14}, \end{equation}

which, considering $V_f\approx 1$ in the inertial-dominated phase, results in a critical transition time $\mathcal{T}\approx \mathcal{L}$ . Figure 5(c) also demonstrates that the critical length exhibits minimal dependence on $\beta$ . Furthermore, for all inclination angles, cases with larger $Fr$ correspond to greater $\mathcal{L}$ .

Our results suggest that the power-law trend length, $\mathcal{L}^{\dagger }$ , observed after the inertial phase (brown lines in figure 5 a), is governed by a balance between characteristic longitudinal buoyant stresses and the corresponding viscous stresses:

(3.5) \begin{equation} \Delta \hat ho \hat g\hat D\cos \beta \sim \frac {{{{\hat \mu }_H}{{\hat V}_i}}}{{\hat D}} \equiv \frac {{{{\hat \mu }_H}{{\left ( {\frac {{\Delta \hat ho }}{{\hat {\overline ho } }}\hat g\hat D} \right )}^{1/2}}}}{{\hat D}}\quad ightarrow \quad R{e_b} = \frac {{{{\hat {\overline ho } }^{1/2}}{{ ( {\Delta \hat {ho } } )}^{1/2}}{{\hat g}^{1/2}}{{\hat D}^{3/2}}\cos \beta }}{{{{\hat \mu }_H}}}, \end{equation}

where ${\hat V}_i = ( ({\Delta \hat {ho }}/{\hat {\overline ho }}) \hat {g} \hat {D} )^{1/2}$ is the characteristic buoyancy-induced inertial velocity and $Re_b$ is the corresponding Reynolds number. This balance between longitudinal buoyancy forces and viscous resistance governs the gradual deceleration of the front in this intermediate regime. Figure 5(d) shows that $\mathcal{L}^\dagger$ increases nearly linearly with $Re_b$ , highlighting the balance between buoyancy and the corresponding viscous stresses in this region. Therefore, $\mathcal{L}^\dagger$ characterises the extent of the viscous–buoyant regime prior to the onset of yield-stress control.

Figure 6. Experimental results. (a) Time-dependent variation in interface height ( $h$ ) versus $z$ . (b–e) Effect of similarity parameters on interface height profiles: (b) $z/t$ ; (c) $(z-t)/t$ ; (d) $z{-}t$ ; (e) $z/\sqrt {t}$ and (f) $z-V_ft$ (i.e. $z-z_f(t)$ , where $z_f$ is the front position). Across all panels, the experimental parameters correspond to figure 3. The colourbars represent $t$ , with darker colours indicating longer times. Dash-dotted and solid line profiles represent the inertial- and viscoplastic-dominated phases, respectively.

Moreover, our results suggest that the sharp decay length, $\mathcal{L}^{\dagger \dagger }$ , may be estimated, albeit very crudely, by balancing the characteristic longitudinal buoyant stress and yield stress:

(3.6) \begin{equation} \Delta \hat {ho } \hat {g} \hat {D} \cos {\beta } \sim \hat {\tau }_{y} \quad ightarrow \quad \frac {\Delta \hat {ho } \hat {g} \hat {D} \cos {\beta }}{\hat {\tau }_{y}} \propto \frac {{Re}\cos (\beta )}{{Fr}^2 B} \equiv \frac {\chi \cos (\beta )}{B}. \end{equation}

Figure 5(e) shows that $\mathcal{L}^{\dagger \dagger }$ increases with axial buoyancy ( $\chi \cos \beta$ ) and decreases with yield stress ( $B$ ), indicating that ${\chi \cos (\beta )}/{B}$ is a relevant governing dimensionless group, although deviations from linearity suggest the proposed form is not fully adequate.

Analysing interface height profiles ( $h$ ) provides insights into interfacial front dynamics. Figure 6(a) illustrates the evolution of $h$ throughout a typical experiment, spanning the inertial-dominated phase to the viscoplastic-dominated phase. During the inertial-dominated phase ( $z \lesssim 2$ ), the first three profiles (identified by dash-dotted lines) maintain similar shapes with consistent spacing at predefined time intervals ( $\Delta t$ ). Upon entering the viscoplastic-dominated phase (identified by solid lines), the interface height increases rapidly before slowing, as reflected by the decreasing profile spacing. This stage corresponds to reduced front velocity as inertial effects wane and effective viscous and buoyant stresses govern the flow. Longitudinal buoyancy stresses, influenced by density differences, geometric inclination and interface slope, drive the flow, while increasing effective viscosity impedes motion, progressively reducing the shear rate. The interfacial motion experiences significant deceleration as the interface nears arrest. The interface evolution finally halts when shear stresses within the viscoplastic layer fall below the yield stress. As discussed in figures 3 and 5, this deceleration and arrest result from decreasing shear rates, increasing effective viscosity, and, as will be seen, expanding unyielded regions within the viscoplastic fluid layer. These factors collectively slow the interface until it halts, as visually seen in figures 5(a) and 6(a).

The interface height profiles in figure 6(a) also reveal outward and upward spreading of the viscoplastic fluid as time progresses. A small inertial bump appears at the interfacial front from the onset, similar to those reported in previous studies, e.g. Taghavi et al. (Reference Taghavi, Alba, Séon, Wielage Burchard, Martinez and Frigaard2012b ), which examined buoyant miscible displacement flows with Newtonian fluids. This inertial bump is a flow feature driven by the balance of buoyancy and local inertia, likely unrelated to the viscoplastic nature of the flow.

We now analyse whether the interface height profiles in both the inertial- and viscoplastic-dominated phases exhibit self-similarity, indicating the dominant physical flow mechanism and simplifying the interfacial flow into a universal model. Several moving reference frames were evaluated: ${z}/{t}$ (figure 6 b); $(z-t)/t$ (figure 6 c); $z {-} t$ (figure 6 d); ${z}/{\sqrt {t}}$ (figure 6 e); and $z - V_f t$ (figure 6 f). The dimensionless frame ${z}/{t}$ suggests linear propagation, with the early-time interface heights (inertial-dominated phase) ‘almost’ collapsing onto one another, especially at the frontal region, but no self-similarity is observed in the longer-time interface heights (viscoplastic-dominated phase). Similarly, $(z-t)/t$ shows a collapse in the inertial-dominated phase for $(z-t)/t \in [-1, 0]$ , with unsteady travelling waves emerging at longer times. The $z {-} t$ frame confirms uniform interfacial front advancement in the inertial-dominated phase, but no self-similarity in the viscoplastic-dominated phase. The ${z}/{\sqrt {t}}$ frame reveals no self-similarity, indicating that a macroscopically diffusive interface does not govern the flow. Only $z - V_f t$ achieves convergence to self-similar curves in both phases, each with distinct master curves. In the inertial-dominated phase, self-similarity is observed roughly in ${z}/{t}$ and $(z-t)/t$ , and more completely in $z {-} t$ and $z - V_f t$ , which coincide due to $V_f \approx 1$ in this phase. In the viscoplastic-dominated phase, only $z - V_f t$ maintains self-similarity. Since $V_f$ varies with time, this indicates an unsteady travelling wave at the interfacial front.

When a heavier fluid is placed above a lighter one in a closed-end pipe, purely Newtonian viscous fluids and viscoplastic fluids exhibit distinct behaviours. For near-horizontal inclinations with Newtonian fluids, density differences drive the heavier fluid downward, displacing the lighter one. This causes slow slumping, eventually leading to a stratified configuration with the heavier fluid below the lighter one (Frigaard & Ngwa Reference Frigaard and Ngwa2004; Taghavi et al. Reference Taghavi, Seon, Martinez and Frigaard2009), as purely viscous stresses cannot halt motion. Even highly viscous Newtonian fluids continue moving, albeit slowly (Etrati & Frigaard Reference Etrati and Frigaard2018). However, with viscoplastic fluids, two key differences may arise. First, due to the yield stress, there are scenarios in which the heavy viscoplastic fluid can remain atop the lighter fluid with a stable interface. Second, even if the viscoplastic fluid slumps beneath the lighter fluid, it may stop at a certain point as the interface elongates, e.g. as observed in our experiments. Thus, as we place the heavy viscoplastic fluid atop the light Newtonian fluid via injection, two critical questions arise: (1) Under what conditions does the interface between the heavy viscoplastic and light Newtonian fluids remain stable or become unstable? (2) What is the maximum slumping extent/length during the interface motion before it is arrested?

Regarding the first question, several foundational experimental and theoretical studies (Frigaard Reference Frigaard1998; Frigaard & Scherzer Reference Frigaard and Scherzer1998, Reference Frigaard and Scherzer2000; Crawshaw & Frigaard Reference Crawshaw and Frigaard1999; Frigaard & Crawshaw Reference Frigaard and Crawshaw1999) have examined the conditions under which viscoplastic fluids initiate slow, viscous slump flows following interface failure, i.e. when the interface becomes unstable. These studies demonstrated that the yield stress can inhibit fluid movement, preventing further motion following an initial disturbance. Through dimensional analysis, they concluded that any criterion defining when the flow becomes stationary must depend solely on the inclination angle, $\beta$ , and the yield number, which is the ratio of the yield stress to the buoyancy stress:

(3.7) \begin{equation} Y = \frac {\hat {\tau }_{y}}{\Delta \hat {ho } \hat {g} \hat {D}} \equiv \frac {B}{\chi } ,\end{equation}

where $Y$ represents the yield number. Both experimental and analytical studies have characterised the yield number ranges governing the initiation and sustained flow of viscoplastic fluids in quiescent states across various flow configurations. Pertinent to this study’s flow configuration, Frigaard & Crawshaw (Reference Frigaard and Crawshaw1999) and Crawshaw & Frigaard (Reference Crawshaw and Frigaard1999) studied interface stability in viscous exchange flows between miscible Bingham fluids with differing densities and rheological properties under density-unstable conditions in inclined, closed-end pipes. Using motion equations, and applying no-slip boundary condition and velocity and stress continuity at the interface, they developed a lubrication-type two-layer stratified flow model to establish a conservative boundary between flow and no-flow conditions:

(3.8) \begin{equation} {Y_m} = {\left ( {{{ [ {{\tau _{y,e}}(\phi )\cos \left ( \beta \right )} ]}^2} + {{ [ {{\tau _{y,n}}(\phi )\sin ( \beta )} ]}^2}} \right )^{\frac {1}{2}}} ,\end{equation}

where $Y_m$ is the modified theoretical yield number, $\tau _{y,e}$ is the tangential component of deviatoric stress (off-diagonal components), $\tau _{y,n}$ is the normal component of deviatoric stress (diagonal components), $\phi$ is defined in radians and varies between $-\pi /4$ and $\pi /4$ , and $\beta$ is the deviation from the vertical. A flow configuration where only one fluid has a yield stress corresponds to $\phi = \pm \pi / 4$ , $\tau _{y,e} = 0.3043$ and $\tau _{y,n} = 0.25$ . Applying (3.8) to this scenario produces the solid line in figure 7(a), which delineates flow and no-flow conditions. Projecting our experimental data onto this $Y$ – $\beta$ plane shows that, at various inclination angles, our results fall below the theoretical yield number separating flow and no-flow regions. Thus, our study remains below the theoretical condition, allowing the injected viscoplastic fluid to flow naturally within the lighter fluid, initiating buoyancy-driven interfacial slumping.

Figure 7. Experimental and theoretical results. (a) Experimental results on the $Y$ – $\beta$ plane. The dividing line marks the theoretical yield number for angles (0 $^\circ$ –90 $^\circ$ ) in a density-unstable configuration based on (3.8). (b) Experimental arrest length ( $L_s^{{Exp}} \times {Re}_L$ ) versus model-determined arrest length ( $L_s^{{Model}} \times {Re}_L$ ), with marker size and colour representing ${Re}_L$ and ${\chi \cos (\beta )}/{B}$ , respectively. Solid line is the 1:1 line. (c) Model-based classification (3.9) of the arrested and unhalted regimes on the $Y \times L$ and $\alpha \equiv {\chi \cos (\beta )}/{B}$ plane. In panel (c), pink and cyan areas represent the arrested and unhalted interfacial flow regimes, with red circles and blue squares indicating experimental data for each regime. Experiments with ${\chi \cos (\beta )}/{B} = 0$ are shown on the left of panel (c).

Regarding the above mentioned second question, Frigaard & Ngwa (Reference Frigaard and Ngwa2004) and Malekmohammadi et al. (Reference Malekmohammadi, Naccache, Frigaard and Martinez2010) studied exchange flows of two viscoplastic fluids in near-horizontal pipes with closed ends, using scaling from thin-film lubrication theory. They identified the minimal surface inclination needed to initiate flow, through applying no-slip conditions at walls, stress and velocity continuity across the interface, and negligible surface tension, assuming miscibility on a short time scale. Their results showed that buoyancy-driven slumping continues until a maximum static slump length is reached, where the shear stresses from density differences, interface slope and inclination balance the yield stresses. In our flow configuration, we adapt the analysis of Frigaard & Ngwa (Reference Frigaard and Ngwa2004), considering the simplification that only one fluid is viscoplastic (see Appendix A). The maximum static or upper-bound slump length for our density-unstable conditions, as per Frigaard & Ngwa (Reference Frigaard and Ngwa2004), can be transformed and expressed as

(3.9) \begin{equation} L_{{s}}= \frac {1}{Y} \int _{0}^{1} \frac {\beta _0(h)}{\beta _c(h) - \alpha \beta _0(h)} \, \mathrm{d}h, \end{equation}

where $\alpha$ is determined by writing a balance between longitudinal buoyant stress and the yield stress of the heavy fluid, i.e. (3.6):

(3.10) \begin{equation} \alpha \equiv \frac {{Re}\cos (\beta )}{{Fr}^2 B}\equiv \frac {\chi \cos (\beta )}{B}, \end{equation}

and the geometrical functions $\beta _0(h)$ and $\beta _c(h)$ , which depend solely on $h$ , are given in Appendix A.

Figure 7(b) compares the theoretical maximum static length – i.e. the axial length over which the interfacial front evolves before it stops (when the interface is arrested) – with corresponding experimental data. The data are adjusted by subtracting the inertial-dominated length ( $\mathcal{L}$ ) to match the inertialess theoretical model. Both axes are scaled by ${Re}_L = {\hat {\overline ho }_L \hat {V}_0 \hat {D}}/{\hat {\mu }_L}$ to enhance the visualisation of data distribution. As shown, our experimental static lengths follow the model, but generally exceed the model’s predictions.

Despite the model’s simplifying assumptions and the challenges associated with viscoplastic fluid experiments, its predictions show reasonable agreement with experimental results, although some discrepancies are observed, particularly where experimental values exceed model predictions. One likely cause is the yield stress value used in the analysis, derived from Herschel–Bulkley model fitting, which differs slightly from values obtained by other methods (see figure 2 and table 2, and compare figures 2 a and 2 b). This reflects the nuanced yielding behaviour of viscoplastic fluids and suggests that steady-state shear rheometry may not always provide the most appropriate yield stress for different flow conditions. This discrepancy influences the value of $L_s^{{Model}}$ in figure 7(b). Other experimental parameters, such as density, inclination and length, may also have minor effects on $L_s^{{Model}}$ . Another source of discrepancy arises from fundamental differences between the simplified model and the richer dynamics of the experimental flow, where several effects are omitted. For example, the model neglects the potential weak influence of injection velocity on the measured $L_s^{{Exp}}$ , which may contribute to deviations observed in figure 7(b) at higher ${Re}_L$ values, where the gap between experimental data and model predictions exceeds that at lower ${Re}_L$ , as seen in the log-log plot.

Note that, while our analysis builds on the framework of Frigaard & Ngwa (Reference Frigaard and Ngwa2004), our configuration involves continuous injection, where the finite volume emerges dynamically rather than being prescribed. This distinction is crucial: the final volume and front position arise from the interplay of inertia, buoyancy and viscoplasticity. Thus, the finite-volume model is used here as an idealised upper bound for viscoplastic arrest, not to fully capture the evolving dynamics. In fact, these differences can help explain the discrepancies between experimental interfacial arrest lengths and theoretical predictions (as shown in figure 7 b).

Following standard methods (Saltelli, Tarantola & Campolongo Reference Saltelli, Tarantola and Campolongo2000; Sobol Reference Sobol2001; Ebtehaj et al. Reference Ebtehaj, Bonakdari, Zaji, Azimi and Sharifi2015; Reed et al. Reference Reed2022), we conducted a sensitivity analysis revealing that model–experiment deviations increase significantly with rising density contrast, decreasing yield stress and steeper inclinations, indicating reduced accuracy when the Boussinesq approximation or yield-dominated assumptions are violated. Exceeding the small-density-difference regime or using very low yield stress leads to a well over 10-fold increase in error, while higher yield stress improves accuracy by a factor of ${\sim}3$ . Shifting inclination from near-horizontal to moderately steep causes a ${\sim}5$ -fold increase in error, consistent with the growing influence of transverse pressure gradients and secondary flows. Injection velocity, although not explicitly included in the model, also correlates with deviations, suggesting residual inertial effects in the viscoplastic-dominated regime. These trends demonstrate that model accuracy is governed primarily by density contrast and yield stress, with injection velocity acting as a secondary but important factor. The model remains quantitatively reliable only within a narrow regime where inertia is negligible, density differences are small and flow remains stably stratified. This aligns with our previous findings (Faramarzi et al. Reference Faramarzi, Akbari and Taghavi2024) that the stable slumping regime is bounded by ${\chi \cos \beta }/{B} \lesssim 40$ and ${Re}/{M} \lt 2800$ . Furthermore, as shown in figure 7(a), very high yield stress may induce a no-flow condition at small ${\chi \cos \beta }/{B}$ . Thus, the model is valid only in flowable regimes, characterised by stable stratification and the absence of secondary flow, with strong agreement for $0.1 \lesssim {\chi \cos \beta }/{B} \lesssim 3$ and ${Re}/{M} \lesssim 1000$ .

Let us now classify the flow regimes based on whether the interface motion becomes arrested within the finite length of our pipe. Considering that both (3.9) and our experimental results indicate that varying conditions produce different static lengths, a large number of experiments were conducted. Experiments where the interfacial motion arrest occurs before the front reaches the pipe end are categorised as the ‘arrested’ regime, while those where arrest does not occur before the front reaches the pipe end are classified as the ‘unhalted’ regime. Using (3.9), the theoretical classification of such regimes is achieved by plotting the $Y \times L_{s}^{{Model}}$ curve for various $\alpha \equiv {\chi \cos (\beta )}/{B}$ values, as shown in figure 7(c), delineating arrested (above the curve) and unhalted (below the curve) regime regions. Our experimental results show that reducing the flow domain length can shift an experiment from the arrested regime to the unhalted regime. Similarly, changes in density difference, yield stress or inclination angle (e.g. an increased inclination shifting a datapoint to the right) can reclassify an experiment from the arrested to the unhalted regime. In general, our experimental results align well with this classification. However, discrepancies arise where some unhalted regime experiments fall within the predicted arrested regime region. This occurs because, as illustrated in figure 7(b), the experimental static lengths are generally greater than those predicted by the model ( $L_s^{{Exp}} \gt L_s^{{Model}}$ ). Depending on the pipe length, this difference may cause the experimentally observed static length to exceed the pipe length, resulting in the classification of a datapoint as belonging to the unhalted regime, even though the interface in that experiment would be arrested in a slightly longer pipe.

Two limiting cases offer further insight into the flow behaviour. First, for a horizontal pipe ( $\beta = \pi /2$ ) with finite injection velocity ( $\hat {V}_0$ ), experiments can exhibit both ‘unhalted’ and ‘arrested’ regimes. In the model, this condition yields $\alpha \approx 0$ , reducing the governing equation to $L_s = {1}/{Y} \int _0^1 {\beta _0(h)}/{\beta _c(h)}\, \mathrm{d}h$ , indicating that, under flowable conditions, natural slumping occurs due to interface height (Frigaard & Ngwa Reference Frigaard and Ngwa2004). Consequently, the model predicts both flow regimes depending on pipe length and other parameters, consistent with experimental observations. Second, in the absence of imposed flow ( $\hat {V}_0 = 0$ ) but at high inclinations, distinct flow configurations emerge, including drainage flows (Akbari, Frigaard & Taghavi Reference Akbari, Frigaard and Taghavi2024) and buoyancy-driven exchange flows (Charabin & Frigaard Reference Charabin and Frigaard2024). In drainage flows, the viscoplastic fluid in the inner pipe drains under gravity into the surrounding Newtonian fluid, with flow onset controlled by density difference, rheology and inclination angle. In exchange flows, vertical pipes with a low-viscosity Newtonian fluid below and a viscoplastic fluid above exhibit flow cessation for large $Y$ , while lower $Y$ values lead to central upward motion and descending viscoplastic displacement, forming regimes such as helical finger, disconnected finger and slug flow. In contrast to these, our study focuses on flowable regimes with non-zero injection velocity, where the viscoplastic fluid penetrates the ambient phase via buoyant and inertial forces.

Figure 8. Experimental results (dimensional) at $\hat {x}=0$ and $\hat {z}= 0.06$ (m). (a) Axial velocity profiles versus $y$ , determined by UDV for the experiment shown in figure 3. Colour intensity represents time [5.2, 8.6, 11.4] (s), with darker curves indicating longer times. (b) Axial velocity profile versus $\hat {y}$ from additional UDV measurements (UDV probe at the lower pipe wall), taken at $\hat {x} = 0$ , $\hat {z}=0.12$ m and $\hat {t}=19$ s for $Re = 2.2$ , $M = 0.003$ , $B = 0.75$ , $\beta = 82^\circ$ and $Fr = 0.68$ . (c) Approximated axial shear stress ( $\hat {\tau }$ ) within the light, Newtonian layer. (d) Shear stress at the interface ( $\hat {\tau }_i$ ), as a function of time (pink curve with circular markers corresponding to the left axis) and $\hat {h}$ (cyan curve with square markers corresponding to the right axis). In panels (a) and (b), colour-matched solid lines on the curves indicate the interface position between the two fluids ( $\hat {h}$ ).

3.3. Arrest mechanism

To further elucidate the interfacial arrest mechanism, this section examines velocity, shear rate and shear stress fields as well as yield/unyielded zones, using a combination of experimental and CFD results.

Figure 8(a) presents the dimensional axial velocity profiles from UDV measurements, corresponding to the experiment shown in figure 3, within $0 \leqslant \hat {y} \leqslant 0.038$ m at $\hat {x} = 0$ and $\hat {z} = 0.06$ m. The profiles darken over time. Positive velocities denote the injection direction, while negative values indicate counterflow, as the light in-place Newtonian fluid moves upward when the heavy viscoplastic fluid slumps beneath it. Near the upper wall, velocities remain near zero due to the no-slip condition, although UDV measurement reflections at the lower wall hinder exact zero-velocity measurements (Alba et al. Reference Alba, Taghavi and Frigaard2013b ). However, a no-slip condition for the Carbopol solution, consistent with prior studies (Alba et al. Reference Alba, Taghavi, de Bruyn and Frigaard2013a ; Alba & Frigaard Reference Alba and Frigaard2016; Freydier et al. Reference Freydier, Chambon and Naaim2017; Liu et al. Reference Liu, Balmforth and Hormozi2018, Reference Liu, Balmforth and Hormozi2019; Ball & Balmforth Reference Ball and Balmforth2021; Valette et al. Reference Valette, Pereira, Riber, Sardo, Larcher and Hachem2021; Isukwem et al. Reference Isukwem, Godefroid, Monteux, Bouttes, Castellani, Hachem, Valette and Pereira2024), was validated through additional UDV measurements near the lower pipe wall, where a zero-velocity region was observed (figure 8 b), confirming negligible wall slip. As time progresses, figure 8(a) shows that an unyielded zone forms (seen as uniform velocity profiles) in the viscoplastic fluid layer due to its yield stress, consistent with previous findings (Kazemi et al. Reference Kazemi, Akbari, Vidal and Taghavi2024). Both injection and counterflow velocities decrease over time, with the pseudo-unyielded region expanding in the $y$ -direction, propagating towards the lower pipe wall, progressively enveloping the pipe’s entire cross-section and reducing axial velocity to approach zero.

Using experiments, shear stresses can be roughly estimated from the shear rate above the interface within the Newtonian layer (but not within the viscoplastic layer) due to the constant viscosity of the Newtonian layer, as seen in figure 8(c). This also allows a rough estimate of the shear stresses at the interface. Figure 8(d) shows that these interfacial shear stresses progressively decrease over time as the interface height grows, until $\hat {h} \approx 0.03$ m (i.e. $h \approx 3/4$ in dimensionless form), indicating heightened susceptibility to unyielded zone formation in this region. In addition, the interfacial shear stresses remain small overall and typically below the viscoplastic fluid’s yield stress (considering the shear stress continuity at the interface), suggesting that the unyielded zones begin to develop in the proximity of the interface, most likely expanding from the neighbourhood of $h \approx 3/4$ . This hypothesis will be further explored by our CFD simulation results.

Figure 9. CFD results. (a) Velocity vectors in the $z$ -direction and (b) cross-sectional $z$ -direction velocity profiles at various times and locations for $Re = 1.8$ , $M = 0.001$ , $B = 0.56$ , $\beta = 82^\circ$ and $Fr = 0.73$ . The colourbar shows velocity magnitude, with positive values for forward flow and negative for counterflow.

The CFD simulations provide detailed insights into the velocity field across the entire flow domain, capturing the advancing interfacial front, regions beyond the front, and both inertial- and viscoplastic-dominated phases, along with insights into the interfacial arrest phenomena. Figure 9(a) shows axial velocity vectors in the pipe centre plane (i.e. the y–z plane at $x = 0$ ). The velocity vectors primarily align in the $z$ -direction, indicating stable flow, especially far from the injection point and the interfacial front. Near the inlet, velocities are much higher, consistent with the injection velocity, and decrease as the front advances. This suggests that unyielded zones would be less pronounced in the inertial-dominated phase than in the viscoplastic-dominated phase. Near the front, counterflow interactions may induce instabilities, suggesting a small yielded region, which diminishes further back, as shown by velocimetry studies (Andreini, Epely-Chauvin & Ancey Reference Andreini, Epely-Chauvin and Ancey2012; Freydier et al. Reference Freydier, Chambon and Naaim2017; Akbari & Taghavi Reference Akbari and Taghavi2022a ,Reference Akbari and Taghavi b ; Faramarzi et al. Reference Faramarzi, Akbari and Taghavi2024). No velocity vectors appear beyond the front, showing that the light, ambient Newtonian fluid remains stationary until it is displaced by the interface as the viscoplastic fluid moves, as also confirmed by UDV measurements (omitted for brevity).

Figure 9(a) shows velocity deviations near the inlet and near-zero velocities downstream at late times, indicating unyielded plug formation following interfacial arrest (as later confirmed in figure 10). Continued injection leads to flow redirection due to mass conservation, with reversal near $z = 0$ observed in both simulations and experiments (figures 3 and 4). This reversal marks a deviation from hydrostatic conditions, with vertical motion implying non-negligible vertical accelerations and non-hydrostatic dynamics during the transition to arrest. Spatially, this localised behaviour occurs in the region where the inertial-dominated regime was initially observed. Also, UDV measurements near the inlet, taken after complete cross-sectional filling with the viscoplastic fluid, confirmed late-time flow reversal (figure omitted for brevity). Similar behaviour was reported by Akbari & Taghavi (Reference Akbari and Taghavi2021) for buoyant miscible Newtonian fluids, where UDV profiles 5 cm above the gate valve showed negative axial velocities towards the outlet.

Figure 9(b) illustrates velocity profiles within the cross-section, across multiple time intervals and axial positions, with the axial position increasing from top to bottom and time progressing from left to right. The profiles display lower velocities near the pipe walls due to the no-slip condition and reveal a significant velocity reduction over time at a given location or at extended times for a specific axial position.

Figure 10. CFD results. (a) Snapshots of numerically derived yield surfaces, with the red line indicating the interface, in the y– $z$ plane at $x=0$ for $Re=1.8$ , $M=0.001$ , $B=0.56$ , $\beta =82^\circ$ and $Fr=0.73$ . (b) Cross-sectional yield surfaces in the $x$ – $y$ plane at various times and positions ( $z=3.4$ , $5.1$ , $6.8$ and $7.7$ ), marked in panel (a) by arrows (dark green, blue, green and orange, respectively).

Figure 10(a) illustrates the evolution of unyielded zones (black) on the pipe’s centre plane (y– $z$ plane at $x = 0$ ), initially forming near the interface (red) where shear stresses are lowest, as also confirmed by UDV measurements (figure 8). These zones expand towards the pipe centre and eventually towards the pipe’s lower wall over time. Near the inlet, where the inertial-dominated phase appears at early times, unyielded regions are minimal. However, moving away from the inlet and as the flow transitions to the viscoplastic-dominated phase, unyielded zones form and progressively grow, increasing the effective viscosity of the viscoplastic layer. This growth causes the interfacial front to decelerate and eventually halt as unyielded regions encompass the entire viscoplastic layer.

An intriguing aspect of the unyielded zone formation onset is observed at $t = 5.5$ , where unyielded zones initially develop in regions with higher interface heights; this is influenced by the height variations along the pipe length ( $z$ ). Over time, the unyielded zones expand to regions with lower interface heights, becoming fully established. Meanwhile, in regions with higher interface heights, the unyielded zones extend towards the pipe centre, but small yielded zones remain near the lower wall, where shear stresses are relatively large. These observations indicate that the interface height influences unyielded zone formation, with higher heights being initially more susceptible, but lower-height regions ultimately exhibiting a higher unyielded-to-yielded ratio.

Figure 10(b) illustrates the evolution of unyielded areas in cross-sectional views at various times and locations, showing that these zones initially form at the interface (at $x \approx 0$ ), where the shear stress is minimal and farthest from the pipe walls, where it is maximal. Over time, and with increasing distance from the injection point, the proportion of unyielded regions grows, reflecting the rising influence of solid-like behaviour and effective viscosity in the viscoplastic fluid. In both the y–z (figure 10 a) and $x$ – $y$ (figure 10 b) planes, the unyielded zones expand continuously without breakage, indicating stratified flow stability within the viscoplastic fluid layer.

Before concluding the manuscript, let us summarise the interpretation of our observations in this section. The arrest mechanism in our buoyant viscoplastic injections arises from the interplay between buoyancy-driven motion, rheological behaviour and geometrical constraints within the inclined pipe. As the flow transitions from an inertial-dominated phase to a viscoplastic-dominated phase, inertia becomes less relevant and buoyancy drives the interface forward. Simultaneously, the yield stress gains relative significance, increasing the viscoplastic fluid’s effective viscosity and resisting the longitudinal buoyancy component. This leads to interface deceleration as the yield stress surpasses the diminishing shear stresses within the viscoplastic layer. Unyielded zones first emerge near the interface, where shear stresses are minimal and progressively expand towards the pipe centre, walls, confining yielded regions to zones of maximal shear near the walls. This progressive reduction in active flow area eventually arrests the interface once shear stresses throughout the viscoplastic layer drop below the yield stress, yielding static equilibrium.