3.1. Exchange flow and no-flow regimes

The stresses generated by buoyancy in the absence of initial motion are strongly dependent on interface configuration (Frigaard Reference Frigaard1998; Frigaard & Scherzer Reference Frigaard and Scherzer1998; Frigaard & Crawshaw Reference Frigaard and Crawshaw1999; Frigaard & Scherzer Reference Frigaard and Scherzer2000). Experimentally, we have not found a reliable way to implement such a (computationally ideal) initial condition. The stresses and initial motion induced by opening of the gate valve are not well quantified, but the opening is automated/controlled by constant air pressure and hence consistent between experiments. Due to this uncertainty, a pragmatic approach is needed to delineate no flow from unstable exchange flows.

In our experiments, the yield-stress fluid is attached to the outer wall, either wholly (vertical) or partially (inclined slumping). In this case, in order to mobilise the fluid at the walls the wall shear stress must exceed the yield stress. Similarly, flow should cease when $\hat \tau _y \gt \hat \tau _w$ . The source of the wall shear stress is buoyancy, which is partitioned between the two fluids. Since we are interested primarily in flows that propagate along the pipe, we consider only the axial component of buoyancy to be dominant. Thus, from a scaling perspective, we expect no flow to be found if

(3.1) \begin{equation} \hat \tau _Y \gt \frac {1}{2} ({\hat \rho _H - \hat \rho _L} )\hat g{\hat R}\cos \beta . \end{equation}

Note that if all of the buoyancy gradient is applied to a single fluid, the right-hand side of (3.1) is the wall shear stress. A similar hypothesis and scaling law has been also used recently for the drainage of a viscoplastic fluid from an inner pipe into an outer closed-end pipe filled by a Newtonian fluid (Akbari et al. Reference Akbari, Frigaard and Taghavi2024).

Considering (3.1) as an approximate threshold between flow and no-flow conditions, the experiments can be categorised into two distinct flow regimes: flow and no flow. For the no-flow regime, the system was observed for at least 3 minutes after opening the gate valve to determine whether any exchange flow would occur. Figure 5 presents this classification in the plane of yield stress of the viscoplastic fluid versus the longitudinal buoyancy stress. In this figure, the circle and star data points represent experiment and simulation data, respectively, with white and yellow symbols, respectively, indicating flow and no-flow cases. The solid line corresponds to $\hat \tau _y = 0.5\Delta \hat \rho \hat g{\hat R} \cos \beta$ . As a guide, the dotted line, $\hat {\tau }_y = \Delta \hat {\rho } \hat {g} \hat {R}\cos \beta$ , or $Y = \cos \beta$ , is also plotted. The results demonstrate that this analysis effectively separates observed flow regimes into two distinct regions, highlighted by different colours. Experimental data points corresponding to the no-flow regime reported in Charabin & Frigaard (Reference Charabin and Frigaard2024) are added using blue symbols.

Figure 5.

Classification of flow (white) and no-flow (yellow) regimes based on analysis made in (3.1). The circle and star symbols represent experimental and simulation data, respectively. The solid line divides the plane into flow (red) and no-flow (blue) regions at $\hat \tau _y = 0.5\Delta \hat \rho \hat g{\hat R} \cos \beta$ ( $Y = 0.5 \cos \beta$ ). The blue data points correspond to no-flow-regime experiments from Charabin & Frigaard (Reference Charabin and Frigaard2024). The dotted line marks $Y= \cos \beta$ , i.e. $\hat \tau _y = \Delta \hat \rho \hat g{\hat R}\cos \beta$ .

As discussed in Charabin & Frigaard (Reference Charabin and Frigaard2024) and shown again here, classifications such as (3.1) provide a sufficient condition for no flow; consequently, falling below this yield-stress threshold establishes a necessary, but not sufficient, condition for flow initiation. Flow ultimately depends on the interface configuration, which generates the local stress. Previous studies have examined the instability of a more dense Bingham fluid positioned above a less dense Bingham fluid in an inclined pipe. The (axial) buoyancy–yield stress balance in (3.1) also arises, but the numerical prefactor depends on the assumed interface configuration. The factor $1/2$ in (3.1) is also that found for concentric interfaces (yield-stress fluid on the outer wall), as in Frigaard & Scherzer (Reference Frigaard and Scherzer1998), while a stratified interface has prefactor $0.6089$ , determined in Frigaard & Scherzer (Reference Frigaard and Scherzer2000). These theoretical predictions, however, tend to be conservative relative to experimental findings (Frigaard & Crawshaw Reference Frigaard and Crawshaw1999), as they must account for the least stable interface, which is not always found experimentally (or numerically).

In the context of off-bottom cement plug placement in oil and gas wells, this leads to a dilemma, if trying to ensure static stability. We can take an expression such as (3.1) as an industrial design rule, but the required yield stresses may make the cement unpumpable! In Charabin & Frigaard (Reference Charabin and Frigaard2024), the estimate $Y\gt Y_{c,exp} = 0.2$ was used, reflecting the observed stability in the experiments. Here we see that many experiments in figure 5 are still static below $Y = \cos \beta$ ; see the filled blue and yellow circles on the right-hand side of the figure. This is inevitable. There is no lower limit on $Y$ that will exclude every possible interface: the horizontal interface in a vertical pipe generates no shear stress when static. Further, many slurries will mix when pumped into a well (Ghazal & Karimfazli Reference Ghazal and Karimfazli2021, Reference Ghazal and Karimfazli2022), meaning that the assumed sharp change $\Delta \hat {\rho }$ across the interface is likely conservative. The practical implication is not that the yield-stress criterion is useless, but rather that it must be applied as an operational risk envelope rather than a deterministic threshold. For example, taking $Y \gt 0.05 \cos \beta$ or $0.1 \cos \beta$ as a design rule might give a better practical prediction of the boundary of flow, but we can see that a number of failures (flows) would result, based on our experiments. Ultimately, this choice depends on the consequences of failed placement.

In figure 5, the empirical transition between the flow and no-flow regimes appears nearly horizontal, particularly at higher longitudinal buoyancy stresses. This indicates that the transition is governed by a critical yield-stress threshold ( $\hat {\tau }_y$ ) rather than the linear theoretical bounds. The theoretical lines therefore represent one-sided upper bounds for stability. At large density differences, the shear stress generated by the initial interface configuration, whether imposed by the gate valve in the experiments or by the initial curved front in the simulations, does not necessarily attain the theoretical maximum. Consequently, fluids with a yield stress exceeding the stress generated by the specific interface geometry remain static. This interpretation is supported by the numerical simulations (star symbols), which also exhibit static behaviour below the theoretical bounds at large buoyancy stresses.

To assess the practical relevance of the three-minute observation window used to define the no-flow boundary, it is useful to compare the characteristic time scale of buoyancy-driven motion in the laboratory and field:

(3.2) \begin{align} \hat {t} \sim \frac {\hat {\mu }\hat {L}}{\Delta \hat {\rho } \hat {g} \hat {D}^{2}}. \end{align}

For a representative field-scale plug (e.g. $\hat {D}=24.4\ \mathrm{cm}$ and $\hat {L}=20\ \mathrm{m}$ ), the strong dependence on $\hat {D}^{2}$ largely offsets the higher viscosity of field cements. As a result, the time scale for buoyancy-driven failure is expected to be of the same order as, or shorter than, that observed in our laboratory experiments. Therefore, if a macroscopic buoyant failure is to occur, it is expected to initiate within the first few minutes after placement. Conversely, if the fluid remains stable during this critical period, subsequent hydration, gel-strength development and thickening increase the effective yield stress and further suppress motion. The three-minute observation window therefore provides a practically relevant criterion for assessing long-term plug integrity.

Figure 6.

Experimental snapshots of the exchange flow at the different regimes. (a) Fingering regime for sample V, $\chi =2.05$ , $Y=0.082$ and $\beta =0^{\circ }$ ; from left to right $t=[390, 940, 1364, 2070, 2834, 3366 ]$ . (b) Slumping regime for sample IX, $\chi =6.98$ , $Y=0.022$ and $\beta =30^{\circ }$ ; from left to right $t=[297, 416, 559, 702, 880, 1011]$ . In all panels, the dimensionless field of view is $2 \times 150$ . The solid black bar(s) in each image represent(s) the pipe support, here and elsewhere.

3.2. General flow behaviours

Let us begin with describing the general behaviour observed. As the more dense fluid is released once the gate valve is opened, the interface evolves in the middle of the pipe, as observed in the upper half of the apparatus. In our experiments, we have observed two distinct flow regimes: fingering and slumping flows.

1. (i) Fingering regime: as shown in figure 6(a), in which only the top half of the pipe is illustrated (above the gate valve), the less dense Newtonian fluid penetrates into the viscoplastic fluid as a central finger, eventually reaching the top end of the pipe, while the more dense fluid descends along the wall into the lower half of the pipe. The dynamics of the viscoplastic flows in the lower half of the pipe is discussed later in § 3.6. This flow regime is observed when the pipe is oriented vertically.

2. (ii) Slumping regime: figure 6(b) shows a buoyant front that advances towards the upper side of the pipe upon opening the gate valve, with the denser viscoplastic fluid slumping downwards below the rising Newtonian front. The front shape appears to advance in a quasi-steady fashion, while instability and mixing are seen far from the front. As the exchange flow continues, the front reaches the top end of the pipe, while the counterflow of the viscoplastic fluid continues towards the lower end of the pipe. The slumping regime occurs when the pipe is deviated from vertical.

Figure 7.

(a,b) Spatiotemporal diagram of the normalised depth-averaged concentration field, $\bar {\bar {c}}(z,t)$ , corresponding to figures 6(a) and 6(b). The colour bar represents the concentration value of the less dense fluid here and everywhere else. The arrow in (a) marks the path of the finger tip. The solid white line in each panel represents distorted pixels due to the presence of the pipe support, here and elsewhere.

Figure 7 shows the average concentration field of the less dense Newtonian fluid in the experiments corresponding to figure 6. The concentration, $\overline {\overline {c}}(z,t)$ , represents the cross-sectionally averaged concentration of the less dense fluid as a function of axial position, $z$ , and time, $t$ . In our coordinate system, $z$ is the axial direction, $y$ is the depth axis aligned with the optical path of the camera and $x$ is the transverse axis across the visible width of the pipe (see figure 1). For a pipe of dimensionless radius $1$ , this true cross-sectional average of the 3-D concentration field $c(x,y,z,t)$ is defined by the double integral over the circular domain:

(3.3) \begin{equation} \overline {\overline {c}}(z,t) = \frac {1}{\pi } \int _{-1}^{1} \left ( \int _{-\sqrt {1-x^2}}^{\sqrt {1-x^2}} c(x,y,z,t)\,\text{d}y \right ) \text{d}x. \end{equation}

Experimentally, this double-averaging is achieved in two steps. The inner integral along the $y$ axis is inherently performed physically by the camera; the captured 2-D pixel intensities represent the depth-integrated light absorption. Subsequently, our image processing code performs the outer integral by averaging these depth-integrated 2-D intensities across the transverse $x$ axis.

This diagram can be used to visualise the flow development, the less dense fluid front evolution and the degree of macroscopic streamwise (axial) mixing between the fluids throughout the flow domain during the exchange flow, as lateral concentration variations are inherently averaged out in this representation.

In figure 7(a), the arrow indicates the less dense Newtonian finger, advancing at constant speed. In this regime, at later times behind the finger, the concentration of the less dense fluid increases in the flow domain ( $t\gt 2100$ ), reflecting the progression of the exchange flow lower down. In contrast, figure 7(b) shows light-coloured streaks behind the finger tip, corresponding to segments of the viscoplastic fluid that descend along the pipe in wave-like fashion, initiating 30–40 radii behind the front. After $t\gt 1200$ , the concentration field remains essentially unchanged, indicating that the exchange flow has reached completion.

Figure 8.

Numerical simulation results for $\chi =7.98$ , $Y=0.022$ and $\beta = 0^{\circ }$ . Snapshots of concentration profile (a) and velocity contour in the z direction (b). Cross-sectional concentration (c) and velocity profiles in the z direction (d) at $z=8$ with increasing time (same as times in a) from left to right. The dimensionless field of view in (a,b) is $2 \times 120$ . The colour bar in (a,c) represents the concentration value and the black regions mark the numerically computed yield surfaces of the viscoplastic fluid. The colour bar in (b,d) denotes the magnitude of the velocity in the z direction. The dashed red rectangles mark the snapshot selected for magnification. The enlarged views clearly show the yielded and unyielded regions along the interface in (a) and the velocity vectors in (b).

Figure 8(a) provides simulation snapshots of the concentration profile for the fingering regime, showing the front development. We have also superimposed regions of unyielded viscoplastic fluid on the snapshots using black colour. For the regularised model these are defined as where the computed stress is below the yield stress, i.e. but the flow still has a small strain rate. As can be seen, far from the finger front, the whole viscoplastic fluid is unyielded. The yielded region of the viscoplastic fluid starts forming near the interface (regions coloured red), where shear stresses are highest. Interestingly, we have also identified a very thin layer of unyielded fluid between the yielded viscoplastic fluid attached to the pipe wall and the finger interface (see magnified snapshot), highlighting the interface deformation and, in particular, the yielded regions surrounding the interface, revealing the localised shear zones on the pipe wall. In the last snapshot, unyielded segments of the viscoplastic fluid are observed far behind the finger, falling down the pipe, as also observed experimentally.

The colourmap in figure 8(b) shows the axial component of velocity (z direction) in the central $x{-}z$ plane at $y = 0$ . The first notable observation is that the maximum velocity occurs in the finger core, away from the fluid–fluid interface and some distance behind the advancing finger front. This 2-D view also reveals that the viscoplastic fluid exhibits a lower velocity magnitude compared with the Newtonian fluid. Since the net flow is zero, the viscoplastic phase occupies a larger cross-sectional area. Focusing on the velocity vectors of the magnified snapshot in figure 8(b), most of the vectors remain aligned in the direction of the flow, i.e. $\pm z$ , except near the finger front, where the velocity must decelerate within the finger and the fluid ahead must move around the sides of the advancing finger.

In figure 8(c), simulations display the colourmap of the more dense viscoplastic fluid concentration at the same times as in figure 8(a), all at $z=8$ radii above the gate valve. Unyielded regions are again shown in black. Initially, at $t=210$ , as the finger width is growing, we can see a small unyielded portion of the viscoplastic fluid. However, in the second snapshot, at $t=287$ , a circular finger is surrounded by unyielded viscoplastic fluid. Subsequently, the finger’s cross-section becomes distorted due to instabilities and mixing induced by the downward motion of the more dense fluid.

Comparing figures 8(a) and 8(c), we see that the initial central finger that advances in an approximately steady and stable fashion is associated with an unyielded sheath of more dense fluid at the interface, over the part of the finger that is uniform. This feature is known to prevent the growth of linear interfacial instabilities (Frigaard Reference Frigaard2001; Moyers-Gonzalez et al. Reference Moyers-Gonzalez, Frigaard and Nouar2004), which would otherwise be expected in a countercurrent flow of two viscous fluids. The importance of an unyielded layer at the interface was also identified in Charabin & Frigaard (Reference Charabin and Frigaard2024) using a one-dimensional (1-D) axisymmetric model, relating the loss of the unyielded sheath to less stable (more inertial) exchange flows that we discuss further below. The simulations allow us to see that this sheath is effectively ruptured as we move behind the front. Since the 1-D flow may in theory extend forever, rupturing of the sheath should be attributed to emergence of stresses that are not significant in the 1-D flow. Figure 8(b) perhaps gives the answer. We see that the Newtonian fluid behind the finger front (entering at $t=287$ ) is moving significantly faster (strong yellow colour) than that at the front. The retardation of the Newtonian stream directly behind the tip requires the conversion of kinetic energy into local gradients in dynamic pressure. Significant inertia is critical to this mechanism: in a strictly viscous regime, such kinematic deceleration would be accommodated by viscous dissipation without generating stresses large enough to break the yield surface. Because the Newtonian flow is inertial, these dynamic pressure gradients are substantial and must be balanced by streamwise variations in interfacial stress. These localised inertially driven stresses ultimately exceed the yield stress of the surrounding fluid, rupturing the protective sheath and leading to a form of helical buckling of the stream and slow growth of the interfacial instability, all at a distance behind the front (see $t=430$ ).

To shed light on the dynamics of the flow in an inclined pipe, the concentration profiles of the fluids with superimposed unyielded viscoplastic fluid regions obtained by the numerical simulation are illustrated in figure 9(a), corresponding to a flow at $\beta =15^\circ$ . As seen, a very thin viscoplastic layer exists on the upper wall of the pipe (right-hand side of images). In other words, even though from the side view an experiment may appear to have a stratified interface, in actuality this may simply be the side view of an eccentrically advancing finger. Interestingly, a thin unyielded layer is found at the interface between the fluids on the lower side of the finger, where the viscoplastic layer is thicker and presumably a reduced shear stress results. Unlike the unyielded sheath in the vertical pipe, the unyielded layer does not surround the advancing Newtonian finger and appears only some distance behind the front.

Figure 9.

Numerical simulation results for $\chi =7.98$ , $Y=0.022$ , and $\beta = 15^{\circ }$ . Snapshots of concentration profile (a) and velocity component and vectors in the z direction (b). Cross-sectional concentration (c) and velocity profiles in the z direction (d) at $z=8$ with increasing time (same as times in a) from left to right. The dimensionless field of view in (a,b) is $2 \times 120$ . The colour bar in (a,c) represents the concentration value and the black regions mark the numerically computed yield surfaces of the viscoplastic fluid. The colour bar in (b,d) denotes the magnitude of the velocity in the z direction.

The axial velocity component colourmap and velocity vectors are depicted in figure 9(b), where we see that the fluids travel with approximately the same speed in opposite directions. One can also observe that the finger tip remains at a small distance from the upper pipe wall, maintaining the eccentric finger configuration, which is further discussed in § 3.4. The velocity vectors in the flow domain are fairly parallel and uniformly aligned in the direction of the downward and upward flows with small velocity fluctuations.

In figure 9(c), the concentration contour of the less dense fluid is illustrated on a 2-D cross-section of the pipe at $z=8$ , at the same times as figure 9(a). We see that the initial finger shape is approximately circular at the centre of the pipe at early times, while it becomes eccentric towards one side of the pipe at later times. The unyielded viscoplastic fluid is superimposed on the concentration field in black. The residual layer of viscoplastic fluid reduces with time. Noting that the less dense Newtonian fluid is pushed upwards by buoyancy, transverse to the flow direction, the wall film is likely squeezed slowly around into the low side of the pipe as the finger advances, i.e. viscous film drainage. A similar flow structure is found in fluid–fluid displacement flows with significant viscosity differences (Etrati et al. Reference Etrati, Alba and Frigaard2018).

In figure 9(d), the velocity contour of the fluids at the pipe cross-section is illustrated. The results indicate that the maximum velocity occurs at the centre of the less dense (Newtonian) fluid and decreases progressively away from the core. The surrounding viscoplastic fluid exhibits lower velocities but occupies a larger fraction of the cross-sectional area. This distribution highlights the confinement effect of the viscoplastic phase, which restricts lateral spreading of the Newtonian finger while simultaneously dominating the overall flow cross-section.

Figure 10.

Experimental snapshots of exchange flow regimes. (a) Helical finger for sample III, $\chi =0.20$ and $Y=0.34$ at $t=[135, 434, 710, 1064, 1450, 2010]$ . (b) Disconnected finger for sample V, $\chi =1.2$ and $Y=0.082$ at $t=[181, 283, 370, 507, 667, 1051]$ . (c) Slug flow for sample X, $\chi =10.46$ and $Y=0.0063$ at $t=[51, 143, 180, 257, 337, 481]$ . In all panels, $\beta = 0^{\circ }$ and the passage of time is from left to right. The dimensionless field of view is $2 \times 150$ .

3.3. Fingering regimes in vertical pipe

In this section, we delve into the detailed dynamics of the fingering phenomena observed during vertical exchange flows. While the macroscopic classification of these flows into three distinct subregimes, namely helical finger, disconnected finger and slug flow, was recently established by Charabin & Frigaard (Reference Charabin and Frigaard2024), the present study significantly extends that foundational work. Here, we move beyond 1-D phenomenological descriptions to provide a rigorous, quantitative analysis of the 3-D internal mechanics and stability thresholds of these regimes. Specifically, we utilise our 3-D numerical simulations to resolve the asymmetric rupturing of the unyielded viscoplastic sheath, and we introduce new empirical metrics to quantify the maximum stable penetration length of the fingers, their dynamic elongation post-detachment and the wavelength of the resulting 3-D corkscrew instabilities. The three subregimes, associated with decreasing $Y$ and increasing inertia as the system moves away from the no-flow boundary, are detailed below:

1. (i) Helical finger: This subregime occurs at relatively high yield stress and low density differences ( $Y$ just below onset), as illustrated in figure 10(a). At short times, the finger grows approximately centrally. Some distance behind the front we first begin to see deformation of the interface in the form of a helical wave; see also our discussion of figure 8. Eventually, the disturbance grows and the finger separates from the rest of the Newtonian fluid. Behind the separated finger, a mixing region develops that expands along the pipe, often preserving a helical structure.

2. (ii) Disconnected finger: Illustrated in figure 10(b). This subregime occurs at moderate yield stress of the more dense fluid and moderate density difference, where the finger detaches from the Newtonian fluid much closer to the gate valve compared with the helical fingering subregime. Unlike the latter, the detachment appears complete: no mixed zone of helical structures is observed connecting to the finger and the detached finger rises more rapidly.

3. (iii) Slug flow: Illustrated in figure 10(c). This subregime occurs at small yield stresses and larger density differences (i.e. when the buoyancy force is relatively large), in which a small finger or ‘slug’ separates just above the gate valve and rises faster than in the other subregimes. Nearer to the gate valve a mixed region is formed. Sometimes further slugs may be formed. In this flow regime, the viscoplastic fluid remains on the pipe wall for a longer time. The presence of the mixed region substantially reduces the buoyancy-driven exchange flow between the fluids. However, it does not entirely halt the flow; exchange becomes intermittent and weaker rather than ceasing altogether.

Figure 11.

(a–c) Spatiotemporal diagrams of the normalised depth-averaged concentration field, $\bar {\bar {c}}(z,t)$ , corresponding to figure 10(a–c), respectively. In (a), the slope of the solid line represents the front velocity ( $V_{\!f}$ ) of the finger and the region enclosed by the ellipse indicates where the finger separates and the concentration values approach zero. In (b), $\hat l_{f,i}$ and $\hat l_{f,e}$ are shown by two arrows. In (c), the arrow marks the upward flow of the mixed region behind the separated finger. The star symbols mark the exact pinch-off coordinate.

Spatiotemporal diagrams associated with figure 10 are presented in figure 11. In figure 11(a), a sharp boundary between different comparative concentrations of the two fluids represents a stable interface at the tip of the finger while moving towards the top end of the pipe. At later times, within the central portion of the upper half of the pipe ( $25\lt z\lt 75$ ), a mixed region is observed, likely due to shear and velocity gradients in this region. The slope of the solid black line in this figure represents the front velocity of the less dense Newtonian fluid, which is quantified in § 3.5. The region enclosed by the ellipse indicates the location where the finger detaches, and the Newtonian fluid concentration approaches zero. Closer to the gate valve, we observe downwards streaks associated with destabilisation and detachment of the wall layer, discussed later.

In figure 11(b), the lengths of disconnected finger at the onset of separation, $\hat l_{f,i}$ , and when it reaches the end of the apparatus, $\hat l_{f,e}$ , are marked by the arrows. A white region is seen behind the finger at longer times, which represents the separation of the finger from the rest of the fluid mixture. Again downwards streaks are observed lower down. In the slug finger regime (figure 11 c) the separation is clear. The path of the small initial slug is observed as a streak line starting at $t=0$ , within which the concentration of the less dense fluid remains below $\approx$ 30 % everywhere, i.e. effectively suggesting the slug occupies less than 55 % of the pipe radius. The mixed region in figure 11(c) moves upwards slowly as indicated by an arrow. Although slower than the slug, spreading of the mixed regime is fast compared with the other regimes. The nearly linear trajectory of the mixed-region front in figure 11(c) indicates an advective rise with approximately constant velocity, rather than a diffusive spreading.

Figure 12(a) illustrates the time of the finger detachment from the fluid, denoted as $t_{f,i}$ , as a function of the buoyancy number ( $\chi$ ) for various flow regimes. This separation time is determined automatically from our image processing data by identifying the earliest time step at which the cross-sectionally averaged concentration, $\overline {\overline {c}}(z,t)$ , drops below a near-zero threshold ( $\approx \kern-1pt 0.05$ ) at an axial location between the gate valve and the finger tip, signalling a complete pinch-off (marked by red star in figure 11). It is evident that the finger detachment occurs at longer times in the helical finger regime, whereas detachment is considerably quicker in the slug flow regime. The corresponding finger lengths at the onset of separation, $l_{f,i}$ , were quantified for each regime and are shown in figure 12(b) as a function of $\chi$ . Consistent with the behaviour of $t_{f,i}$ , the finger length in the helical finger regime is greater than in other flow regimes. Although we have classified the flows using the three subregimes, as also in Charabin & Frigaard (Reference Charabin and Frigaard2024), it is evident that the transition between regimes as $Y$ decreases and $\chi$ increases is broadly speaking defined by an increase in (buoyancy-driven) inertial stresses compared with viscous (including yield) stresses. Some of the classifications may also be experiment-dependent, such as whether the eventual disconnection of all helical fingers is plausible in a sufficiently long pipe.

Figure 12.

(a) Finger separation time versus $\chi$ . (b) Finger length at separation for different $\chi$ . The colour bar shows $Y$ . Different symbols mark helical finger ( $\blacklozenge$ ), disconnected finger ( $\blacktriangleright$ ) and slug ( $\bullet$ ) flow regimes. Here and in subsequent figures, error bars denote the maximum relative experimental uncertainty.

Focusing on the finger regimes, when the less dense Newtonian fluid penetrates into the viscoplastic medium, the process unfolds through a distinct sequence of stages. As shown in figure 13(a), the less dense Newtonian fluid initially advances into the yield-stress medium as a finger of nearly constant width. For sufficiently short times, this interface remains coherent and the finger tip stable. Beyond a critical penetration length (equivalently, distance behind the tip), the interface destabilises, eventually leading to a mixed region with entrainment and irregular interfacial structures.

To quantify this transition, we define the maximum stable penetration length, $l_s$ , as the distance travelled by the finger before the onset of visible interfacial instabilities. This length was systematically measured across a range of density contrasts and viscoplastic properties, and the results were plotted against the dimensionless buoyancy number, $\chi$ , as illustrated in figure 13(b). Each point represents a different experimental condition, with the colour indicating the yield number, $Y$ , according to the colour bar. The data reveal that $l_s$ decreases with increasing $\chi$ , demonstrating that higher buoyancy forces reduce the distance over which the finger remains stable. The black dash-dotted line indicates a trend capturing the overall relationship between $l_s$ and $\chi$ as $l_s=(37 \pm 5) \chi ^{-0.5}$ .

Figure 13.

(a) Typical fingering regime experiment. (b) Stable finger length, $l_s$ , versus $\chi$ . The dash-dotted line marks $l_s=(37 \pm 5) \chi ^{-0.5}$ . (c) Time to onset of interface instability, $t_s$ , versus $\chi$ . The dash-dotted line represents $t_s=(1951 \pm 320)\chi ^{-0.5}$ . In (b,c), the colour and size of the symbols represent $Y$ .

Figure 13(c) presents the corresponding time at which the maximum stable length is reached, $t_s$ , also plotted against $\chi$ and colour-coded by $Y$ . The results show a similar trend: $t_s$ decreases as $\chi$ increases, indicating that destabilisation occurs faster under stronger buoyancy forces or weaker viscous (viscoplastic fluid) resistance. The dash-dotted line again represents a trend capturing the overall trend between $t_s$ and $\chi$ as $t_s=(1951 \pm 320)\chi ^{-0.5}$ . The analysis reveals a quantitative criterion for the onset of interfacial instability: both the maximum stable length and the corresponding time decrease systematically with increasing buoyancy number.

While a formal stability derivation is highly complex for this 3-D yield-stress flow, this inverse square-root scaling for $l_s$ and $t_s$ , reminiscent of typical viscous shear-layer growth (Moyers-Gonzalez et al. Reference Moyers-Gonzalez, Frigaard and Nouar2004), effectively bounds the experimental data, providing a practical predictive estimate for the maximum stable length. The practical benefit of providing these specific scaling fits, along with their standard deviations to account for natural experimental scatter and measurement uncertainty (rather than any numerical error, as the simulations are strictly converged), is that they offer a quantifiable guideline for estimating the operational stability window, in both time and distance, before the onset of detrimental interfacial mixing.

To quantify the finger length during upward propagation, we have measured its length, extracted from the time-series images, when it is separated from the bulk of the fluid as well as when it reaches the top end of the pipe, as marked, respectively, by $\hat {l}_{f,i}$ and $\hat {l}_{f,e}$ in figures 11(b) and 14(a). Figure 14(b) shows the variation of $\hat {l}_{f,i}$ and $\hat {l}_{f,e}$ with $Y$ , indicated by circle and diamond symbols, respectively.

Figure 14.

(a) Illustration of finger length at the separation point, $\hat l_{f,i}$ , and at the top end of the pipe, $\hat l_{f,e}$ . (b) Plots of $l_{f,i}$ and $l_{f,e}$ for different $Y$ . Symbol colour and size denote $\chi$ .

The first observation is that the finger length increases during ascent, i.e. $\hat {l}_{f,e} \gt \hat {l}_{f,i}$ for each constant $Y$ . As the Newtonian finger passes through the viscoplastic fluid, we have seen that there is yielding only immediately in front. The viscoplastic fluid is displaced as the finger rises; it travels down along the sides of the finger, yielding at the pipe wall, often retaining an unyielded sheath of fluid around the finger. Thus, radial expansion is restricted. If the volume of the Newtonian finger is fixed after detachment, any reduction in cross-sectional area must result in elongation.

The second observation is that the length of the Newtonian finger (both $\hat {l}_{f,i}$ and $\hat {l}_{f,e}$ ) increases with increasing $Y$ , i.e. increasing yield stress and decreasing buoyancy force, which can be explained by a simple scaling argument. Consider that the finger generates an axial pressure head across its own length, given by

(3.4) \begin{align} \Delta \hat p \sim \Delta \hat \rho \hat g \hat l_{\!f} , \end{align}

where $\Delta \hat \rho$ is the density difference between the finger and the surrounding viscoplastic fluid. This pressure head produces interfacial shear stresses in the surrounding fluid. At the onset of yielding, the yield stress of the viscoplastic fluid balances the stress generated by the buoyant head, i.e. $\hat \tau _y \;\approx \; \Delta \hat \rho \, \hat g \, \hat l_{\!f}$ . Rearranging gives the characteristic length of the finger as

(3.5) \begin{align} \hat l_{\!f} \approx \frac {\hat \tau _y}{ \Delta \hat \rho \hat g}, \end{align}

or in dimensionless variables $l_{\!f} \propto Y$ . This scaling argument shows that increasing yield stress leads to a longer finger, while decreasing buoyancy also results in an increased finger length. However, as seen in figure 14(b), while the data confirm the qualitative trend that finger length increases with $Y$ , the actual dependence is much weaker than linear. This discrepancy arises because the finger is highly dynamic; viscous stresses and interfacial mixing actively carry much of the buoyant load, rendering a purely static yield balance insufficient for quantitative prediction.

Figure 15.

(a) Corkscrew pattern behind the finger for two experiments: $\chi =3.26$ (left) and $\chi =1.78$ (right). (b) Length of corkscrew, $l_c$ , versus $\chi$ . The dash-dotted line marks $l_c=5.72-0.27\chi$ . The colour and size of the symbols represent $Y$ .

Since the finger volume, $\hat {V}=\; \pi {\hat r_{\!f}}^2 \hat l_{\!f}$ , is conserved, we find that the cross-sectional radius ( $\hat r_{\!f}$ ) decreases with increasing length:

(3.6) \begin{align} {\hat r_{\!f}}^2 = \frac {\hat {V}}{\pi \hat l_{\!f} } \propto \frac {\Delta \hat \rho \hat g}{\hat \tau _y}, \end{align}

i.e. $r_{\!f} \propto Y^{-1/2}$ dimensionlessly. Thus, higher yield stress or lower buoyancy causes the finger to become narrower in cross-section and longer in the axial direction, consistent with the experimental results.

Generally, we do observe that the finger is lighter at the exit than initially. As our measurement relies on scaled image pixel values, lighter image can either mean smaller finger radius or possibly a change in concentration within the finger. The former of these might e.g. result from yield stress confinement, and we do see a reduced radius in the images. The latter could arise from slow entrainment of viscoplastic fluid at the interface, into the finger. Consistent with the simulations in figure 8, which show dispersive mixing within the Newtonian finger, we cannot dismiss this possibility, the result of which is that the observed detached finger volume may not be strictly constant. This, of course, also leads to an increase in length.

We have also observed rotating periodically buckled structures, corkscrew patterns, formed behind the finger, as shown in figure 15(a). The emergence of corkscrew-shaped patterns during the upward penetration of the Newtonian finger into the viscoplastic fluid can be attributed to the asymmetric yielding of the surrounding viscoplastic fluid, as depicted in figure 8(c). As the focused low-viscosity Newtonian fluid finger rises, the local buoyancy force, driven by the density difference, generates strong interfacial shear stresses on the surrounding viscoplastic medium. Where these buoyant stresses locally exceed the yield stress, they create a thin, yielded layer that acts as a lubricating sheath around the finger tip. However, the yielding does not occur uniformly around the finger circumference. Small perturbations, arising from geometric imperfections, flow instabilities or other asymmetries, are likely exacerbated by the nonlinear yielding process of the viscoplastic fluid. As evidenced by the asymmetric unyielded zones in our numerical simulations (e.g. figure 8 c), a slight initial lateral perturbation concentrates the interfacial shear stress on one side of the finger. Because of the fluid’s yield-stress and shear-thinning nature, this side yields more readily and experiences a drastic drop in local effective viscosity, producing an imbalance in the viscous resistance and pressure field that generates a net lateral force and torque on the finger. This leads to a helical or corkscrew-like motion of the advancing interface. At early times, when the driving pressure gradient and local shear rates are high, this asymmetric yielding and associated rotational motion are more pronounced (see the velocity variation across the finger in figure 8 b). As the flow evolves, the yielded region around the finger thickens and becomes more symmetric, resulting in a more uniform velocity distribution across the finger and a gradual attenuation of the corkscrew pattern. The phenomenon thus reflects the dynamic coupling between interfacial shear, localised yielding and pressure redistribution in a confined viscoplastic–Newtonian exchange flow.

To characterise this pattern, we have systematically measured the corkscrew wavelength, $\hat l_c$ . As shown in figure 14, the Newtonian fluid undergoes continuous longitudinal elongation under the confining shear of the surrounding viscoplastic fluid as it rises. Consequently, the helical pattern inscribed on the interface is kinematically stretched, leading to a visibly larger apparent wavelength at axial locations further up the pipe (behind the finger tip). To isolate the fundamental instability at its onset, before it is distorted by this macroscopic elongation, we have applied a strict measurement protocol, in which $\bar {l}_c$ is defined as the average length of the first two consecutive, fully formed wavelengths immediately above the gate valve. This metric was then averaged over several frames during steady propagation, and the resulting dimensionless values are plotted versus the buoyancy number ( $\chi$ ) in figure 15(b).

As illustrated in figure 15(b), with increasing $\chi$ , i.e. increasing the density difference or decreasing the viscous stress, the corkscrew wavelength decreases. Such a pattern has not been seen for $\chi \gt$ 4. Looking at the range of $\chi$ in figure 12, we see that this pattern happens exclusively in the helical fingering regime, suggesting its evolution from the initial helical interfacial instabilities found as the descending viscoplastic layer destabilises. The amplitude is evidently limited by the pipe radius.

Corkscrew patterns have been observed in previous multiphase flow studies (e.g. Bai et al. Reference Bai, Chen and Joseph1992), and helical modes arise naturally when destabilising axisymmetric flows. Alba et al. (Reference Alba, Taghavi, de Bruyn and Frigaard2013) found corkscrew patterns in displacement flows of Carbopol: a similar flow except with a non-zero net flow. From a practical viewpoint, the corkscrew pattern enhances local mixing and may also accelerate contamination of the cement slurry locally, underscoring the importance of controlling density contrasts and rheology to suppress such secondary instabilities. Viewed differently, if the exchange flow cannot be arrested, corkscrew and similar structures in mixed regions effectively reduce the width of channels through which the two fluids flow and exchange. Thus, as long as the mixed regions are confined (and predictable), the operational impact may be managed, i.e. pump more cement.

3.4. Slumping regime in inclined pipe

This section examines the slumping regime in inclined pipes. A buoyant exchange flow develops with the dense viscoplastic fluid slumping beneath the rising Newtonian fluid and moving down the lower wall of the pipe. Figure 16 shows flow snapshots at different inclination angles from vertical. As can be observed in figure 16(a), initially, a quasi-steady interface of nearly constant height with limited mixing is formed. At later times, behind the front, the Newtonian fluid’s height reduces at some points as wavelike instabilities grow. The combination of strong shear across the interface and the longitudinal gravitational component in the inclined pipe leads to the formation of these periodic wave-like structures. Morphologically, these resemble the long-wave instabilities and interfacial roll waves recently detailed in sloping Newtonian exchange flows (Zhu et al. Reference Zhu, Atoufi, Lefauve, Kerswell and Linden2024). However, unlike those fundamentally inertial wave phenomena, the growth and propagation of the structures in this viscoplastic system remain heavily damped and kinematically controlled by the high effective viscosity and yield stress of the heavy phase. This shear-driven interfacial evolution ultimately leads to localised mixing between the fluids. Eventually, the advancing front reaches the pipe end. Increasing deviation from vertical (figure 16 b) leads to the creation of a more stable Newtonian fluid layer in comparison with figure 16(a). Intuitively, this happens due to weakening longitudinal buoyancy forces and an increased stabilising transverse buoyancy at the interface, in a more inclined pipe.

Figure 16.

Experimental snapshots of the slumping flow regime for sample IX, $\chi =9.73$ and $Y=0.022$ at different inclination angles: (a) $\beta =15^{\circ }$ at $t=[58, 139, 216, 420, 604, 710]$ and (b) $\beta =30^{\circ }$ at $t=[34, 120, 205, 390, 505, 602]$ .

Figure 17.

(a,b) Spatiotemporal diagram of the normalised depth-averaged concentration field, $\bar {\bar {c}}(z,t)$ , corresponding to figures 16(a) and 16(b), respectively.

Figure 17 presents the spatiotemporal diagrams corresponding to figure 16. The developed mixing region behind the front is more pronounced in figure 17(a) for $\beta = 15^{\circ }$ . In both subfigures, distinct streak lines can be observed, indicating the downward motion of the viscoplastic fluid. Increasing the inclination to $\beta = 30^{\circ }$ results in greater accumulation of the less dense fluid in the upper section of the pipe ( $z \gt 100$ ). Furthermore, the figure demonstrates that reducing $\beta$ , thereby enhancing axial buoyancy stresses, promotes stronger mixing between the fluids. Interestingly, the front speed is slightly larger at $\beta =30^\circ$ .

Figure 18.

(a) Illustration of an eccentric finger and residual layer of viscoplastic fluid on the pipe wall marked by dashed green lines; the solid and dashed black lines represent the pipe centre and finger centre, respectively. (b) Eccentricity of the fingers, $e$ , versus inclination angle, $\beta$ . Symbol colour and size represent the dimensionless finger radius.

Although at first glance in the experiments the slumping flows appear stratified, the advancing finger tip remains at a small distance from the upper pipe wall. There is a thin viscoplastic fluid layer on the wall, which can extend a significant distance behind the advancing tip. This feature is sometimes difficult to observe in the experiments due to the pipe curvature, but it is evident in the numerical results. See figure 18(a) for a comparison with numerical results and also our earlier example in figure 9. In the cross-sections plotted in figure 9, we see that the finger effectively advances initially with an eccentrically positioned circular shape. Later, as the finger advances, it becomes more stratified, adapting to the shape of the pipe and squeezing the wall layer to slowly drain. To further quantify the initial deviation of the finger from the pipe centreline, we introduce a finger eccentricity parameter defined as

(3.7) \begin{equation} e = \frac {\Delta \hat r}{\hat R}, \end{equation}

where $\Delta \hat r$ is the radial displacement of the finger tip from the pipe centre and $\hat R$ is the pipe radius. Eccentricity was extracted from numerical simulations. To avoid sensitivity to small local fluctuations, we report $e$ measured at $z=0.2$ and averaged over a short time window around the instant shown in figure 18(a).

Figure 18(b) presents the variation of the eccentricity parameter, $e$ , with the pipe inclination angle, $\beta$ . The colour bar represents the dimensionless finger radius. The results show that $e$ increases monotonically with $\beta$ , indicating that the finger tip deviates progressively farther from the pipe centreline as the pipe is inclined. At small inclinations ( $\beta \leqslant 15^\circ$ ), the finger remains nearly axisymmetric ( $e \approx 0$ ), and the dimensionless finger radius varies within a narrow range. For larger inclinations ( $\beta \gt 15^\circ$ ), the finger becomes increasingly eccentric ( $e \gt 0.1$ ), while the radius decreases progressively, indicating narrower and more confined finger structures. This trend is consistent with buoyancy-driven segregation flows as inclination increases, gravity pulling the less dense Newtonian fluid preferentially towards the upper side of the pipe, reducing the distance of the finger tip from the pipe wall (reducing $e$ ) and causing the finger to rise eccentrically with a more confined cross-section.

A limitation of the experimental flow visualisation is noted: the side-view cameras provide depth-averaged concentration fields that obscure the 3-D eccentricity of the interface. As a result, 2-D projections may overestimate apparent mixing or misrepresent the true finger width in the slumping regime. This limitation highlights the need for complementary 3-D simulations. As shown in figures 9 and 18, CFD reveals the internal flow structure and captures the strong eccentricity and transverse interface displacement not accessible experimentally.

Figure 19.

Time-dependent variation in interface height ( $h$ ) of less dense fluid slumping in an inclined pipe over $t = [30, 48, 80, 105]$ (s) for sample IX, $\beta =30^\circ$ , $\chi =9.75$ and $Y=0.022$ , versus $z-V_{\!f} t$ .

The interface height $h$ ( $z,t$ ) in figure 19 was obtained from experimental measurements. Because our optical set-up captures a 2-D projection of the flow, $h$ represents the apparent projected thickness of the light Newtonian fluid rather than a localised 3-D measurement. For each side-view image, $h$ was defined as the local thickness of the less dense Newtonian fluid, measured from the upper wall of the pipe to the interface position, which is identified by tracking the boundary where the path-integrated (depth-averaged) Newtonian-fluid concentration drops to 0.05 (5 %). The profiles shown therefore represent the spatial variation of $h$ ( $z,t$ ) along the pipe axis, from the gate valve ( $z=0$ ) to the less-dense-fluid front at different times.

Analysing interface height profiles ( $h$ ) of the slumping flow also provides insights into interfacial front dynamics. Figure 19 illustrates the evolution of the Newtonian fluid height $h$ throughout a typical experiment, spanning from gate valve ( $z=0$ ) to the less-dense-fluid front for different times. As can be seen, at all times, the interface height oscillates remarkably along the pipe far from the front. However, close to the front the oscillations are minimum. The observed behaviour, apart from the eccentric finger cross-section, is reminiscent of the destabilisation of the helical finger flows earlier.

We now analyse whether the interface height profiles exhibit self-similarity, indicating the dominant physical flow mechanism and simplifying the interfacial flow into a universal model. To determine the dominant physical mechanism governing the interfacial flow, the interface height profiles were evaluated across several reference frames, including the stationary laboratory frame ( $z$ ), an advective frame ( $z/t$ ) and a diffusive frame ( $z/\sqrt {t}$ ). None of these frames collapsed the data, indicating that the interface is not governed by macroscopic diffusion or simple linear stretching. Instead, as shown in figure 19, only the reference frame moving with the advancing tip ( $z-V_{\!f} t$ ) achieves convergence to self-similar curves at almost all times; this indicates a steady travelling Newtonian fluid front as we will also discuss for figure 20.

As shown in figure 19, the interface near the leading front evolves into a steady travelling-wave structure. Visually, this steady-profile transition bears a morphological resemblance to the finite-amplitude travelling waves observed in viscoelastic channel flows (Buza et al. Reference Buza, Beneitez, Page and Kerswell2022). However, unlike elasto-inertial phenomena, the evolution of the interfacial structures here is strictly kinematically controlled and damped by the large effective viscosity and yield stress of the viscoplastic phase. The balance between the buoyant driving force and the stabilising yield stress allows this sharp monoclinal-like front to propagate without continually steepening or dispersing.

Figure 20.

Time-dependent position of the less-dense-fluid front along the pipe axis ( $z$ ). (a) Effect of inclination angle for sample IX, $\chi =9.75$ and $Y=0.022$ ; the arrow shows the increasing trend of $\beta =0^\circ,\ 15^\circ,\ 30^\circ$ . (b) Effect of buoyancy numbers ( $\chi$ ) for $\beta =0^\circ$ ; the arrow shows the increasing trend of $\chi =0.20, 3.56, 6.74$ corresponding to $Y=0.50, 0.013, 0.0003$ . In both panels, the darker line represents the higher value of the quantity under study.

Also nicely evident from figure 19 is the growth in amplitude of the interfacial perturbations with distance behind the advancing finger front. This suggests a linear interfacial instability. Because the front advances at a steady velocity into undisturbed fluid, the axial distance behind the tip corresponds directly to the advective time over which the interface has experienced the underlying shear flow. At the leading front, the Newtonian finger continuously encounters undisturbed viscoplastic fluid, meaning the temporal evolution of the instability initiates from a stable state at the tip. Due to the exceptionally high effective viscosity of Carbopol, the temporal growth rate of the interfacial instability is extremely small. As a result, perturbations require a prolonged period of sustained shear, translating to a substantial advective distance, to amplify to an observable magnitude. In particular, figure 19 appears to show that this sustained instability becomes visible at an approximately constant distance of around 40 radii behind the front. While this physical mechanism explains the spatial delay in wave growth, the complexity of the eccentric 3-D interface geometry makes a formal stability analysis highly challenging.

3.5. Front dynamics

In this section, we elaborate on the dynamics of the Newtonian fluid front. First, we have quantified the time-dependent position of the front from the gate valve until it reaches the top of the pipe. Figure 20(a) plots the time-dependent position of the Newtonian fluid front along the pipe for three different values of $\beta$ . The first observation is that, with increasing $\beta$ , the front reaches the pipe end more quickly. Interestingly, during the early stage of the flow ( $z \lt 20$ ), all curves follow a similar trend; beyond this point, the front advances faster at larger $\beta$ . This may be attributed to enhanced transverse mixing at lower $\beta$ , which reduces the effective density contrast of the mixture and consequently diminishes the longitudinal buoyancy force. Also, in the inclined pipe, the fluids are more segregated, which leads to reduced mixing and consequently a faster exchange flow. In figure 20(a), the steeper slopes observed at lower $\beta$ correspond to slower front velocities compared with those at higher  $\beta$ .

Figure 20(b) shows that increasing the buoyancy number enhances the front velocity. For the lowest $\chi$ , the front initially advances slowly $z\lt 30$ , which is due to the development and adjustment of the front at lower buoyancy forces. It can be seen that the front position variation close to the gate valve is not similar at different conditions, which can be due to the dominance of the buoyancy and viscous stress at the beginning of the flow. In general, figure 20(b) shows that the front velocity exhibits stronger time dependence at lower $\chi$ , whereas in most other cases it appears nearly time-independent, likely due to the negligible influence of inertia.

Figure 21.

(a) Dimensional average front velocity of the less dense fluid, $\hat {\bar {V}}_{\!f,L}$ , versus $\hat \mu _H$ . The colour of symbols shows $\Delta \hat \rho = \hat \rho _H - \hat \rho _L$ . (b) Dimensionless average front velocity, ${\bar {V}}_{\!f,L}$ , versus $\chi \cos \beta$ . The line marks $\bar V_{\!f}=0.01 (\chi \cos \beta )^{1.45}$ . The colour of the symbols represents the values of $Y$ . In both panels, different symbol shapes represent $\beta =0^\circ$ ( $\bullet$ ), $\beta =15^\circ$ ( $\blacktriangle$ ) and $\beta =30^\circ$ ( $\blacklozenge$ ).

Figure 21(a) presents the average front velocity of the less dense fluid, $\hat {\bar {V}}_{\!f,L}$ , defined as the velocity averaged over the entire section of the pipe above the gate valve (i.e. the full spatial domain of $z$ shown in figure 20), prior to the arrival of the less-dense-fluid front at the pipe’s upper end, as a critical parameter governing the exchange flow dynamics. The data are plotted against the viscosity of the surrounding viscoplastic fluid, $\hat {\mu }_H$ , measured in our experiments. In this figure, the symbol size and colour encode the values of $\beta$ and the density difference, $\Delta \hat {\rho } = \hat {\rho }_H - \hat {\rho }_L$ , respectively. The results indicate that $\hat {\bar {V}}_{\!f,L}$ generally decreases with increasing $\hat {\mu }_H$ and $\Delta \hat {\rho }$ , while higher front velocities are associated with lower $\hat {\mu }_H$ and larger $\beta$ .

While our current observations show the front speed increasing with the deviation angle ( $\beta$ ), this variation is ultimately non-monotonic. Extensive studies on miscible exchange flows (Séon et al. Reference Séon, Hulin, Salin, Perrin and Hinch2005; Alba et al. Reference Alba, Taghavi and Frigaard2012) demonstrate that front velocity reaches a global maximum at intermediate inclinations (typically $60^\circ \le \beta \le 80^\circ$ ) before decreasing sharply as the pipe approaches horizontal ( $\beta \to \pi /2$ ). In the present study, these extreme deviations were not investigated; our parameter space is intentionally restricted to smaller angles to reflect the primary industrial application of this work, off-bottom plug placement in well abandonment operations, which predominantly involves near-vertical wellbore configurations.

To examine these trends more systematically, the variation of $\hat {\bar {V}}_{\!f,L}$ is recast in dimensionless form in figure 21(b), where the dimensionless front velocity, $\bar {V}_{\!f,L}$ , is plotted as a function of $\chi \cos \beta$ . In this representation, symbol size and colour correspond to $\beta$ and $Y$ , respectively. Overall, $\bar {V}_{\!f,L}$ increases with $\chi \cos \beta$ , with the highest velocities observed at larger values of $\chi$ and $\beta$ , and smaller values of $Y$ . The data are well described by the empirical correlation, $\bar {V}_{\!f,L}=0.01 (\chi \cos \beta )^{1.45}$ , as depicted by the dotted line in the figure. While the exact constants are empirical, this scaling can be physically rationalised by comparison with the inertial limit. For purely Newtonian exchange flows, Séon et al. (Reference Séon, Hulin, Salin, Perrin and Hinch2005) found that the dimensionless front velocity plateaus at an $O(1)$ constant ( $\approx 0.7$ ). In our system, the exceptionally small prefactor ( $0.01$ ) reflects the severe retardation of the flow caused by the viscoplastic fluid’s yield stress and high effective viscosity. However, as the buoyant driving force increases, the viscoplastic resistance breaks down nonlinearly due to fluid shear-thinning and the shrinking of unyielded plug regions. This simultaneous reduction in flow resistance allows the velocity to recover rapidly towards the inertial Newtonian limit, resulting in the super-linear exponent ( $1.45$ ).

A broad spectrum of correlations has been developed to predict the front velocity of miscible buoyant displacement flows across a variety of fluid systems, including Newtonian, shear-thinning and viscoplastic fluids (Sher & Woods Reference Sher and Woods2017; Akbari & Taghavi Reference Akbari and Taghavi2020; Eslami et al. Reference Eslami, Akbari and Taghavi2022; Akbari & Taghavi Reference Akbari and Taghavi2023; Kazemi et al. Reference Kazemi, Akbari, Vidal and Taghavi2024). These correlations consistently demonstrate that the characteristic velocity in gravity-driven flows scales with the buoyancy-inertial velocity, $\hat {V}_i = \sqrt {2 At \hat g \hat R}$ . This scaling provides a useful approximation of front propagation. Building on this framework, Séon et al. (Reference Séon, Hulin, Salin, Perrin and Hinch2005) showed that front velocities $V_{\!f}$ , of Newtonian–Newtonian exchange flow in a pipe, increased with $\beta$ to a plateau: $\hat V_{\!f}=0.7 \sqrt {2 At \hat g \hat R}$ , for $65^{\circ }\lt \beta \lt 82^{\circ }$ .

To understand the internal mechanics driving these front velocities, figure 22 presents the axial velocity profiles across the pipe diameter at $z=80$ , obtained purely from our numerical simulations under different flow conditions. In figure 22(a), the effect of pipe inclination is shown. For the vertical configuration ( $\beta =0^\circ$ , blue line), the upward-moving Newtonian finger remains centred, exhibiting a Poiseuille-like velocity distribution, while the downward-moving viscoplastic fluid (negative $U_z$ ) displays a plug-like profile (flattened plug region identified by arrows) across the pipe cross-section except close to the walls and centre of the pipe where it is yielded due to the higher applied shear and interfacial stresses. When the inclination increases to $\beta =15^\circ$ (red line), the velocity distribution becomes asymmetric: the Newtonian phase accelerates preferentially along one side of the pipe while the viscoplastic phase is displaced non-uniformly on the opposite side. At higher inclination ( $\beta =30^\circ$ , green line), this asymmetry intensifies, with the Newtonian upward flow confined near one pipe wall and the flattened plug region of the viscoplastic downward motion concentrated towards the centre of the pipe.

Figure 22.

Axial velocity profiles across the pipe diameter at $z=80$ for (a) different inclination angles $\beta =0^\circ$ (blue), $\beta =15^\circ$ (red) and $\beta ={30}^\circ$ (green) for $\chi =9.75$ and $Y=0.022$ , (b) different yield stresses quantified by $Y=[0.003, 0.02, 0.036]$ at $\beta =15^\circ$ and (c) different density differences quantified by $\chi =[9.72, 13.73, 18.47]$ at $\beta =15^\circ$ . The colour intensity of the lines in (b,c) represents the values of $Y$ and $\chi$ , respectively; i.e. a darker line represents a higher $Y$ in (b) and a greater $\chi$ in (c). Arrows show the extent of plug regions.

Figure 22(b) illustrates the influence of yield stress at $\beta =15^\circ$ . Lower yield-stress values (lighter lines) lead to a more fluid-like behaviour of the viscoplastic phase, resulting in stronger downward shear-driven velocities. As the yield stress increases (darker lines), the viscoplastic fluid resists deformation, and the velocity distribution flattens, characteristic of plug flow. Figure 22(c) highlights the effect of density differences at $\beta =15^\circ$ . A smaller density contrast (lighter lines) yields more balanced counterflow between phases, while larger density differences (darker lines) enhance buoyancy-driven segregation, leading to greater velocity asymmetry and intensified fingering. These findings are consistent with earlier observations on buoyancy-driven instabilities in multiphase displacement flows (Dai et al. Reference Dai, Eslami, Schneider, Liu and Schwering2024).

In figure 22, the location of the zero velocity ( $U_z=0$ ) changes systematically with the flow parameters and therefore serves as an indication of the fluid–fluid interface. In figure 22(a), the profiles at $\beta =0^\circ$ exhibit two $U_z=0$ points located near the lateral walls, whereas at $\beta =15^\circ$ these points migrate towards one side, evidencing interface displacement and flow asymmetry induced by the pipe tilt. In figure 22(b), increasing $Y$ shifts the zero-velocity location(s) closer to the wall: for low $Y$ a single crossing is observed at the middle of the pipe, while at the largest $Y$ an additional zero-velocity point appears adjacent to the wall, consistent with a thicker plug region and a narrowed shear zone. In figure 22(c), increasing $\chi$ initially displaces the zero-velocity point towards the less dense fluid, reflecting buoyancy-driven stratification of the interface. The zero-velocity position then remains unchanged for further increases in $\chi$ . This is not coincidental, but rather indicates that the flow has reached an asymptotic strongly segregated regime. At high $\chi$ , the transverse location of the interface becomes geometrically locked by the global zero-net-volumetric-flux constraint of the closed pipe ( $Q_H + Q_L = 0$ ). Consequently, while the velocity profile extrema (maximum and minimum values) continue to scale up with the buoyant driving force, the interface location is stabilised, and the shear distribution and flow intensity within each fluid layer are still affected by the density contrast.

3.6. Dynamics of viscoplastic flow below the gate valve

Figure 23 presents the dynamics of the more dense viscoplastic fluid released upon opening the gate valve into the lower half of the pipe initially filled with dyed Newtonian fluid. Since the Newtonian phase is dyed for contrast, the invading viscoplastic fluid itself is not directly visible; instead, the flow evolution is visualised through colour-coded concentration snapshots, as shown in figure 23( $a_1$ , $a_2$ ), where lighter colours denote regions of lower Newtonian-fluid concentration. In the vertical configuration ( $\beta =0^\circ$ ), the viscoplastic phase exhibits a dripping-like motion, where buoyancy forces act to detach small viscoplastic fluid chunks that fall intermittently along the pipe. This fragmented descent is characteristic of yield-stress fluids within confined geometries, in which the fluid resists deformation until local stresses exceed the yield threshold, at which point small packets of fluid break up and descend. The spatiotemporal map (figure 23 $b_1$ ) highlights this process, with discontinuous and step-like fronts (light-coloured downward streaks) associated with successive detachments.

Figure 23.

(a $_1$ ,a $_2$ ) Sequence of experimental snapshots of the flow below the gate valve for sample IX, $\chi =9.75$ and $Y=0.022$ ; the dimensionless field of view is $2 \times 150$ . (a $_1$ ) $\beta =0^\circ$ ; from left to right, $t=$ [20:125:1700]. (a $_2$ ) $\beta =30^\circ$ ; from left to right, $t=$ [10:90:1300]. (b $_1$ ,b $_2$ ) The spatiotemporal diagram of the normalised depth-averaged concentration field, $\bar {\bar {c}}({z},t)$ , corresponding to (a $_1$ ,a $_2$ ), respectively.

When the pipe is inclined at $\beta =30^\circ$ , the flow mechanism changes substantially due to the combined effects of gravity and wall interaction. In this case, the viscoplastic fluid preferentially adheres to and descends along the lower wall, forming larger coherent masses rather than small dripping droplets. The yield stress promotes cohesive motion, while the inclined geometry enhances gravitational spreading along the boundary, resulting in a more continuous displacement. It can also be seen that the falling viscoplastic chunks become larger with the passage of time. The spatiotemporal diagram (figure 23 $b_2$ ) supports this interpretation, with smoother, slanted concentration bands that indicate the sustained descent of large fluid portions. At longer times ( $t\gt 900$ ), the bottom of the pipe is full of viscoplastic fluid, indicated by the white region in figure 23 $(b_2)$ . This contrast between vertical and inclined configurations underscores the strong sensitivity of viscoplastic displacement to system orientation: while vertical flows are governed by localised yielding and droplet detachment, inclined flows promote wall-guided cohesive transport, leading to more effective replacement of the less dense phase.

Figure 24 quantifies the settling dynamics of the viscoplastic fluid fragments released below the gate valve. The settling velocity is extracted from the spatiotemporal diagram (e.g. figure 23 $b_1$ , $b_2$ ) by fitting straight lines to the distinct descending streaks and averaging the slopes ( $\Delta z / \Delta t$ ) across multiple representative fluid fragments. This figure presents the dimensionless settling velocity, $V_s$ , plotted against the buoyancy number, $\chi$ . Despite variations in inclination and yield stress, the majority of the data cluster around a nearly constant value of $V_s \approx 0.65$ , as marked by the dotted line. This collapse of the data implies that, once normalised by the buoyancy-driven scaling, the settling velocity becomes largely independent of the system parameters. Such a result highlights the influence of buoyancy in controlling descent, with wall inclination and yield effects primarily modulating the fragmentation process rather than the normalised settling rate itself. Overall, these findings demonstrate that while the absolute settling velocity depends strongly on buoyancy, the dimensionless dynamics exhibits a universal character.

Figure 24.

Dimensionless settling velocity versus longitudinal buoyancy $\chi \cos \beta$ . The size and colour of symbols represent $\beta$ and $Y$ , respectively. The dotted line marks $V_s\approx 0.65$ .