2. Theoretical model with ascending HB fluid

We consider two immiscible and incompressible fluids moving in opposite directions in a circular pipe, with a rheology described by the following general constitutive equation:

(2.1)\begin{equation} \boldsymbol{\mathsf{T}} ={-}p \boldsymbol{\mathsf{I}} + 2 \eta \boldsymbol{\mathsf{D}} ={-}p \boldsymbol{\mathsf{I}} + \boldsymbol{{\tau}}, \end{equation}

where $p$ is the pressure, $\boldsymbol{\mathsf{I}}$ is the unit tensor, $\boldsymbol{\mathsf{D}}=(1/2)(\boldsymbol {\nabla } \boldsymbol {v}+\boldsymbol {\nabla } \boldsymbol {v}^T)$ is the strain-rate tensor, $\boldsymbol {{\tau }}= 2 \eta \boldsymbol{\mathsf{D}}$ is the deviatoric tensor and $\eta$ is the apparent viscosity. In particular, in the present section, the first fluid is taken to be described by the HB model (Herschel & Bulkley Reference Herschel and Bulkley1926), given by

(2.2)\begin{equation} \left. \begin{array}{ll@{}} |\boldsymbol{{\tau}}|\le \tau_y & \textrm{if}\ | \dot{\boldsymbol{\gamma}}|=0, \\ \boldsymbol{{\tau}} = \left(\mu_0|\dot{\boldsymbol{\gamma}}|^{n-1}+ \dfrac{\tau_y}{|{{\dot{\boldsymbol{\gamma}}}}|}\right) \dot{\boldsymbol{\gamma}} & \textrm{if}\ | \dot{\boldsymbol{\gamma}}|\ne 0, \end{array}\right\} \end{equation}

where $\mu _0$ is the consistency index, $n$ is the fluid behaviour index, $\tau _y$ is the yield stress (positive), $|\dot {\boldsymbol {\gamma }}| = \sqrt {(\dot {\boldsymbol {\gamma }}:\dot {\boldsymbol {\gamma }})/2}$ and $|\boldsymbol {{\tau }}|=\sqrt {({{\boldsymbol {\tau :\tau }}})/2}$ are a (positive) measure of the shear rate and of the stress, respectively, with $\dot {\boldsymbol{\gamma}}:\dot {\boldsymbol{\gamma}}=\dot {{\gamma }}_{ij} \dot {{\gamma }}^{ij}$ and ${{\boldsymbol {\tau :\tau }}}=\tau _{ij}\tau ^{ij}$; in the latter expressions, the Einstein notation applies. If $\tau _y=0$, the HB model describes an OdW fluid (Morrell & de Waele Reference Morrell and de Waele1920; Ostwald Reference Ostwald1929), if $n=1$, (2.2) represents a Bingham fluid (Bingham Reference Bingham1922) and, if $\tau _y=0$ and $n=1$, (2.2) models a Newtonian fluid. The assumption of an ascending HB fluid holds for the entire § 2 and is relaxed in § 3.

The second fluid is always taken to be Newtonian, hence $\eta =\mu _d$, where $\mu _d$ is the viscosity independent of the shear rate.

The flow problem is fully described by the continuity equation under the incompressibility constraint

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0, \end{equation}

and by the linear momentum balance equation

(2.4)\begin{equation} \rho \dfrac{\partial {\boldsymbol{v}}}{\partial {t}} + \rho (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{v} ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} + \rho \boldsymbol{b},\end{equation}

where $\boldsymbol {v}$ is the velocity field and $\boldsymbol {b}$ is the body force per unit mass, reducing to gravity only in many applications.

Let us consider a vertical pipe where the inner ascending current, of lower density, has the said HB rheology and the heavier descending current, in contact with the pipe walls, is Newtonian; see figure 1. Both currents fall under the buoyancy–viscosity regime. We further assume that the ascending current, or long drop, has a length scale parallel to the $z$ axis equal to $\mathcal {L}$ and a radial scale equal to the pipe radius $R$, with an aspect ratio $\varepsilon =R/\mathcal {L}\ll 1$. To derive the vertical velocity scale $U$, we equate, under steady-state conditions, the buoyancy ${\rm \Delta} \rho gR$ to the viscous forces at the pipe wall, $\mu _dU/R$, yielding

(2.5)\begin{equation} {\rm \Delta} \rho gR \sim \frac{\mu_d U}{R},\quad \textrm{or}\quad U=\frac{{\rm \Delta} \rho g R^2}{\mu_d}, \end{equation}

where ${\rm \Delta} \rho =\rho _d-\rho _a$, and $\rho _a$ and $\rho _d$ are the density of the lighter HB ascending fluid and of the denser Newtonian descending fluid, respectively. By continuity, the radial velocity scale is $UR/\mathcal {L}\equiv \varepsilon U$, the time scale is $\mathcal {L}/U\equiv (R/U)/\varepsilon$ and the pressure scale is $\mu _dU\mathcal {L}/R^2\equiv \mu _dU/(\varepsilon R)$.

Figure 1. Schematic of a long drop in a circular pipe. The descending fluid (red) moves near the walls with an annular cross-section, the ascending fluid moves near the pipe axis. Here, $z_d$ and $z_a$ are the coordinates of the descending and ascending fronts, $\mathcal {L}=z_a-z_d$ is the length of the drop, $U_d$ and $U_a$ are the front speeds of the two fluids.

Dimensionless relevant parameters are the density ratio $\mathcal {R}\equiv \rho _a/\rho _d<1$, the Archimedes number ${Ar}=\rho _d{\rm \Delta} \rho g R^3/\mu _d^2$, equal to the ratio between buoyancy and viscous forces, the Reynolds number ${Re}=(\rho _dR\sqrt {g R {\rm \Delta} \rho /\rho _d})/\mu _d$, the viscosity contrast

(2.6)\begin{equation} \mathcal{M}=\frac{\mu_d}{\mu_0(U/R)^{n-1}}\equiv \frac{\mu_d}{\mu_0}\left(\frac{{\rm \Delta} \rho gR}{\mu_d}\right)^{1-n},\end{equation}

i.e. the ratio between the scales of the viscous Newtonian and power-law stresses, and the modified Bingham number for HB fluids

(2.7)\begin{equation} {Bm}=\frac{\tau_y}{\mu_0(U/R)^{n}},\end{equation}

the ratio between the yield stress and the scale of the shear component of the tangential stress. Large values of $\mathcal {M}$ indicate a dominant viscosity in the descending Newtonian fluid, while large values of $Bm$ indicate a dominant contribution of the yield stress to the dynamics of the ascending HB fluid.

The mass balance and momentum equations in a cylindrical geometry for the two fluids can be derived in dimensionless form as reported in Appendix A. In particular, if $\mathcal {M}\gg \varepsilon$, (A11)–(A14) reduce at $O (1)$ to

(2.8)\begin{align} \left.\begin{array}{@{}ll@{}} \displaystyle \dfrac{1}{\tilde{r}}\dfrac{\partial {}}{\partial {\tilde{r}}}\left(\tilde{r} \dfrac{\partial {\tilde{u}_d}}{\partial {\tilde{r}}}\right) = \dfrac{\partial {\tilde{p}_d}}{\partial {\tilde{z}}} +\dfrac{1}{1-\mathcal{R}}, & \tilde{r}\in[\tilde{\delta},1],\\ \begin{aligned} &\displaystyle\dfrac{1}{\mathcal{M}}\dfrac{1}{\tilde{r}}\dfrac{\partial {}}{\partial {\tilde{r}}} \left(\tilde{r}\left|\dfrac{\partial {\tilde{u}_a}}{\partial {\tilde{r}}}\right|^{n-1} \dfrac{\partial {\tilde{u}_a}}{\partial {\tilde{r}}}\right)\\ &\quad +\,{sgn}\left(\dfrac{\partial {\tilde{u}_a}}{\partial {\tilde{r}}}\right) \dfrac{{Bm}}{\mathcal{M}}\dfrac{1}{\tilde{r}} = \dfrac{\partial {\tilde{p}_a}}{\partial {\tilde{z}}} + \dfrac{\mathcal{R}}{1-\mathcal{R}},\end{aligned} & \textrm{if}\ |{\tau_{rz}}| > \tau_y,\ \tilde{r}\in[0,\tilde{\delta}],\\ \tilde{u}_a=\textrm{const.}, & \textrm{if}\ |{\tau_{rz}}|\le \tau_y,\ \tilde{r}\in[0,\tilde{\delta}],\\ \displaystyle \dfrac{\partial {\tilde{p}_d}}{\partial {\tilde{r}}}=\dfrac{\partial {\tilde{p}_a}}{\partial {\tilde{r}}}=0, \end{array}\right\} \end{align}

for the descending and the ascending fluid, respectively. The symbol $\widetilde {\cdots }$ indicates a dimensionless variable and $\delta$ is the radius of the internal current.

Neglecting the interface tension, the pressure assumes a unique value in the horizontal cross-section for the descending and ascending fluids, $p_d=p_a=p(z)$, and (2.8) can be written as

(2.9)\begin{equation} \left.\begin{array}{@{}ll@{}} \displaystyle \dfrac{1}{{r}}\dfrac{\textrm{d}}{\textrm{d}r} \left({r}\dfrac{\textrm{d}{u}_d} {\textrm{d}{r}}\right) =\underbrace{\dfrac{\textrm{d}p}{\textrm{d}z} +\dfrac{1}{1-\mathcal{R}}}_{\mathcal{P}}, & {r}\in[{\delta},1],\\ \begin{aligned}&\displaystyle\dfrac{1}{\mathcal{M}}\dfrac{1}{r}\dfrac{\textrm{d}}{\textrm{d}{r}} \left({r}\left|\dfrac{\textrm{d}{u}_a}{\textrm{d}{r}}\right|^{n-1} \dfrac{\textrm{d}{u}_a}{\textrm{d}{r}}\right)\\ &\quad +\,{sgn}\left(\dfrac{\textrm{d}{u}_a}{\textrm{d}{r}}\right) \dfrac{{Bm}}{\mathcal{M}}\dfrac{1}{r} = \underbrace{\dfrac{\textrm{d} p}{\textrm{d}z} +\dfrac{\mathcal{R}}{1-\mathcal{R}}}_{\mathcal{P}-1},\end{aligned} & \textrm{if}\ |{\tau_{rz}}| > \tau_y,\ {r}\in[0,{\delta}],\\ u_a=\textrm{const.}, & \textrm{if}\ |{\tau_{rz}}| \le \tau_y,\ {r}\in [0,{\delta}], \end{array}\right\} \end{equation}

where the tilde has been dropped for simplicity. The boundary conditions are

(2.10)\begin{equation} \left.\begin{array}{c@{}} u_d(1)=0,\\ u_a(0)\ \textrm{is finite},\\ u_d(\delta)=u_a(\delta),\\ \displaystyle\left.\dfrac{\textrm{d}u_d}{\textrm{d} r}\right|_{r=\delta}= \dfrac{1}{\mathcal{M}}\left|\dfrac{\textrm{d}u_a}{\textrm{d} r}\right|^{n-1}_{r=\delta} \left.\dfrac{\textrm{d}u_a}{\textrm{d} r}\right|_{r=\delta}+ \dfrac{{Bm}}{\mathcal{M}}{sgn}\left(\left.\dfrac{\textrm{d}u_a}{\textrm{d} r}\right|_{r=\delta}\right), \end{array}\right\} \end{equation}

representing (i) the no-slip condition at the wall, (ii) the finiteness of velocity at the axis and the continuity of (iii) velocity and (iv) shear stress $\tau _{rz}$ at the interface between the two fluids.

Solving the differential problem yields the following velocity profile:

(2.11)\begin{align} \left.\begin{array}{ll@{}} u_d(r)=\dfrac{\mathcal{P}}{4}(r^2-1)-\dfrac{\delta ^2}{2}\ln{r}, & {r}\in[{\delta},1],\\ \begin{aligned} u_a(r) & =\dfrac{n}{n+1}\dfrac{[ \mathcal{M} (1-\mathcal{P})\delta-2 Bm ]^{{(n+1)}/{n}}-[\mathcal{M} (1-\mathcal{P}) r-2 Bm ]^{{(n+1)}/{n}}}{ \mathcal{M} (1-\mathcal{P})2^{1/n}}\\ & \quad +\dfrac{\mathcal{P}}{4}(\delta^2-1)-\dfrac{\delta^2}{2}\ln \delta, \end{aligned} & {r}\in[r_y,{\delta}],\\ u_{a({plug})}\equiv u_a(r_y), & {r}\in[0,r_y], \end{array}\right\} \end{align}

where $\mathcal {P}=\textrm {d}p/\textrm {d}z+1/(1-\mathcal {R})$ is the driving force of the ascending fluid and $r_y=2{Bm}/[\mathcal {M}(1-\mathcal {P})]$ is the radius of the axial plug. Imposing a zero net flux at the generic cross-section of the pipe by setting

(2.12)\begin{equation} \int_0^\delta 2{\rm \pi} r u_a\,\textrm{d}r+\int_\delta^1 2{\rm \pi} r u_d\,\textrm{d}r=0, \end{equation}

we obtain an equation involving $\mathcal {M}, Bm, \mathcal {P}, \delta$ and $n$, to be solved numerically in order to evaluate $\mathcal {P}$.

For the case of a OdW fluid (a HB fluid with null yield stress, ${Bm}=0$) the ascending and descending velocities reduce to

(2.13)\begin{equation} \left. \begin{array}{ll} \displaystyle u_d(r)=\dfrac{\mathcal{P}}{4}(r^2-1)-\dfrac{\delta ^2}{2}\ln{r}, & {r}\in[{\delta},1],\\ \begin{aligned} \displaystyle u_a(r) & =\dfrac{n}{n+1}\left[\dfrac{\mathcal{M} (1-\mathcal{P})}{2}\right]^{{1}/{n}}(\delta^{{(n+1)}/{n}} -r^{{(n+1)}/{n}}) \\ & \quad +\displaystyle\dfrac{\mathcal{P}}{4}(\delta^2-1)- \dfrac{\delta^2}{2}\ln \delta, \end{aligned} & {r}\in[0,{\delta}], \end{array}\right\} \end{equation}

and the zero net flux condition in the generic cross-section, (2.12), becomes

(2.14)\begin{equation} \frac{(1-\delta ^4) \mathcal{P}}{8}- \dfrac{n \delta ^{(3n+1)/n}}{(3 n+1) }\left[ \dfrac{\mathcal{M} (1-\mathcal{P})}{2}\right]^{{1}/{n}}- \frac{ \delta ^2 (1-\delta ^2)}{4}=0.\end{equation}

Equation (2.14) admits an analytical solution for $n=1, 1/2, 1/3, 1/4$ (and for $n=2, 3, 4$, a shear-thickening fluid). For $n=1$ it results in

(2.15)\begin{equation} \mathcal{P}=\delta ^2 \frac{(2 \delta ^2-\delta ^2 \mathcal{M} -2)}{\delta ^4-\delta ^4 \mathcal{M}-1}, \end{equation}

coincident with the solution derived by Picchi et al. (Reference Picchi, Suckale and Battiato2020). For $n=1/2$ the closed-form solution is

(2.16)\begin{equation} \mathcal{P}=1+\dfrac{{5} (1-\delta ^4)}{4 \delta ^5 \mathcal{M}^2}-\dfrac{\sqrt{5}(1-\delta ^2)}{4 \delta ^5 \mathcal{M}^2}{ \sqrt{5 \delta ^4+10 \delta ^2+8 \delta ^5 \mathcal{M}^2+5}}. \end{equation}

The analytical solutions for $n=1/3, 1/4$ are cumbersome and are not shown; a numerical solution to (2.14) can be derived for any value of the fluid behaviour index.

Figure 2(a) shows the overall driving term $\mathcal {P}$ as a function of the long drop radius $\delta$ and figure 2(b) does the same for the pressure gradient; different values of the viscosity ratio $\mathcal {M}$, density ratio $\mathcal {R}$ and flow behaviour index $n$ are considered. On the one hand, the driving term increases with the drop radius, more rapidly for lower $\delta$ values and larger viscosity ratios $\mathcal {M}$, until a quasi-linear increase is reached for unit $\mathcal {M}$; more shear-thinning fluids determine a more rapid increase of $\mathcal {P}$ for larger viscosity ratios, while the influence of rheology is modest for lower viscosity ratios. On the other hand, the pressure gradient increases with the drop radius and also, significantly so, with the density ratio $\mathcal {R}$, while it remains largely unaffected by the fluid rheology. Figure 3(a) shows the velocity profiles computed for a Newtonian and a shear-thinning OdW fluid, for different values of $\mathcal {M}$ and for $\delta =0.6$. Figure 3(b) does the same for the case of an ascending HB fluid, with the presence of a modest plug as a consequence of the non-zero value of ${Bm}$ selected; the velocity profile is almost unaffected by the presence of the plug. A backflow of the inner current is observed in both cases as a consequence of the drag of the descending outer current. An increasing viscosity contrast $\mathcal {M}$ between the internal and the external fluid determines an increasing backflow in the ascending fluid; shear-thinning fluids show a flattened velocity distribution near the axis, more so for low values of $\mathcal {M}$.

Figure 2. An OdW ascending fluid ($Bm = 0$) and Newtonian descending fluid. (a) Dimensionless $\mathcal {P}$ as a function of the long drop radius $\delta$ and of the viscosity ratio $\mathcal {M}$ for a fluid behaviour index of the ascending fluid $n=1,1/2,1/3$; (b) dimensionless pressure gradient as a function of $\delta$ and of the density ratio $\mathcal {R}$ for $\mathcal {M}=10$. The inset shows the domain $\delta -\mathcal {M}$ where the $\mathcal {P}$ for Newtonian fluid exceeds $\mathcal {P}$ for a shear-thinning fluid with $n=1/2$ or $1/3$: the two curves cannot be distinguished.

Figure 3. Velocity profiles for a generic configuration with $\delta =0.6$ for different values of the viscosity ratio $\mathcal {M}$, and for flow behaviour indexes $n=1,0.7$, (a) for OdW fluids (${Bm}=0$) and (b) for HB fluids with ${Bm}=0.05$. The thicker part of the curves near the axis indicates the plug.

Figure 4 shows the radius of the plug for different properties of the HB fluids constituting the internal ascending current. The dashed grey area represents the fluid domain subject to shear for a fixed $Bm$ value: the domain becomes larger for small values of $Bm$. The plug radius decreases for increasing $Bm$ and for more shear-thinning fluids.

Figure 4. Radius of the plug $r_y$ as a function of the radius $\delta$ of the inner ascending current for different values of its properties $n$ and $Bm$ and for $\mathcal {M}=10$. The grey area represents the shearing zone for a HB fluid with $Bm=1$.

The shear stress at the interface is shown in figure 5 as a function of the interface position $\delta$ for different values of $Bm$ and $n$ of the ascending inner fluid and a fixed viscosity ratio. The (negative) shear stress decreases with increasing $\delta$, reaches a minimum and then increases to zero for $\delta =1$. Conversely, for a given $\delta$ the magnitude of the tangential stress increases with $Bm$ and is greater for shear-thinning than for Newtonian fluids.

Figure 5. Tangential stress at $r=\delta$, the interface between the inner and the outer current, as a function of the fluid behaviour index $n$ and of the Bingham number $Bm$ of the inner ascending fluid, for $\mathcal {M}=10$.

Figures 6(a) and 6(b) show the velocity and the tangential stress profiles for an ascending non-Newtonian fluid with zero (OdW) and non-zero (HB) yield stress, respectively. A backflow occurs in both cases, with the non-Newtonian fluid dragged down by the Newtonian one at their interface, while a plug is evident for the HB fluid. The velocity profiles associated with the HB fluid are flatter due to the presence of the plug. The tangential stress shows a linear behaviour, attaining the minimum value at the drop radius, the maximum value at the outer wall and a null value where the downward Newtonian velocity is larger in absolute value.

Figure 6. Velocity profiles (blue hatched) and tangential stress (green hatched) radial distributions for a Newtonian descending fluid and (a) an OdW ascending shear-thinning fluid (${Bm}=0$) and (b) a HB ascending fluid (${Bm}=0.2$). The parameter values are $n=0.5, \mathcal {M}=10, \delta =0.6$.

The average ascending and descending speed of the long drop is obtained by averaging the velocity profiles (2.11) along the cross-section

(2.17a,b)\begin{equation} U_a=\dfrac{1}{{\rm \pi}\delta^2}\int_0^\delta 2{\rm \pi} ru_a(r)\,\textrm{d}r, \quad U_d=\dfrac{1}{{\rm \pi}(1-\delta^2)}\int_\delta^12{\rm \pi} r u_d(r) \,\textrm{d}r. \end{equation}

Since the zero net flux condition, (2.12), can be written as $U_aA_a+U_dA_d=0$, this gives

(2.18)\begin{equation} U_a={-}\dfrac{A_d}{A_a}U_d=\dfrac{\delta^2-1}{\delta^2}U_d. \end{equation}

Figure 7 shows the speed of the ascending and descending fronts as a function of the long drop radius $\delta$ and of the viscosity contrast $\mathcal {M}$ for an ascending fluid that is Newtonian ($n=1$) and shear-thinning OdW fluid (${Bm}=0$) with $n=1/2, 1/3$. It is seen that the speed of the ascending fluid depends on the long drop radius in a non-monotonic fashion, first increasing from zero, reaching a maximum and then decreasing towards zero as $\delta$ goes from zero to unity. The ascent speed strongly increases for a large viscosity contrast $\mathcal {M}$ until $\delta <0.8$, while it is modestly affected by the fluid rheology. The downward speed of the descending Newtonian fluid shows a similar behaviour, except the influence of the rheology of the ascending fluid is relatively more impactful. Further, the $n$-index modulates the speed, and it can happen that a current with $n<1$ shows a higher ascent speed than one with $n=1$. To exemplify this case, figure 8 shows the $\delta -\mathcal {M}$ shaded domain where an OdW ascending fluid with $n=1/2$ or $1/3$ generates a long drop faster than a HB ascending fluid with $n=1$ and $Bm = 0.0, 0.1, 1.0$. The domain is almost rectangular with $\delta$ ranging from $0.05$ to $0.55$ for $\mathcal {M}>80$, while its width narrows for $\mathcal {M}<80$, until at $\mathcal {M}\approx 10.4$ for both $n=1/2$ and $1/3$, the required condition is satisfied only for a single value, $\delta =0.32$. Increasing $Bm$ results in a shrinking of the domain.

Figure 7. Long drop average speed for an OdW ascending fluid (${Bm}=0$). (a) Ascent speed and (b) descent speed as a function of the long drop radius $\delta$ and the viscosity ratio $\mathcal {M}$ for $n=1,1/2,1/3$.

Figure 8. Ranges $\delta$–$\mathcal {M}$ where the ascent speed for an OdW fluid with $n=1/2$ and $n=1/3$ is larger than the ascent speed for a HB fluid with $n=1$ and ${Bm}=0, 0.1, 1$.

The limit condition $\mathcal {M}\rightarrow \infty$ is equivalent to a descending Newtonian fluid having a viscosity much larger than the viscosity (for a Newtonian ascending fluid) or the apparent viscosity (for a HB ascending fluid) of the inner current. The asymptotic speeds are

(2.19)\begin{equation} \left. \begin{array}{l@{}} U_d(\delta,\mathcal{M}\rightarrow\infty,{Bm},n)={-} \dfrac{\delta^4(3-4\ln\delta)+1-4\delta^2}{8(1-\delta^2)},\\ U_a(\delta,\mathcal{M}\rightarrow\infty,{Bm},n)= \dfrac{\delta^4(3-4\ln\delta)+1-4\delta^2}{8\delta^2}, \end{array}\right\} \end{equation}

and do not depend on ${Bm}$ or the fluid behaviour index $n$. This result was expected, since in (2.9), the limit $\mathcal {M}\to \infty$, with finite $Bm$, leads to the condition $\mathcal {P}=1$ with the homogeneous boundary condition $\textrm {d}u_d/\textrm {d}r=0$ at $r=\delta$, which allows integration of the velocity profile of the Newtonian descending current without involving any rheological parameter of the ascending fluid. When $\mathcal {M}\to \infty$ and $Bm$ is of the same order as $\mathcal {M}$, the boundary conditions for the descending fluid reduces to $\textrm {d}u_d/\textrm {d}r=- Bm /\mathcal {M}$ at $r=\delta$, and the asymptotic speeds are

(2.20)\begin{align} \left. \begin{array}{l@{}} U_d(\delta,\mathcal{M}\rightarrow\infty,{Bm}=O(\mathcal{M}),n) ={-}\dfrac{\delta^4(3-4\ln\delta)+1-4\delta^2}{8(1-\delta^2)}+ \dfrac{Bm}{\mathcal{M}}\dfrac{1-\delta^2}{8\delta},\\ U_a(\delta,\mathcal{M}\rightarrow\infty,{Bm}=O(\mathcal{M}),n) =\dfrac{\delta^4(3-4\ln\delta)+1-4\delta^2}{8\delta^2}- \dfrac{Bm}{\mathcal{M}}\dfrac{(1-\delta^2)^2}{8\delta^3}, \end{array}\right\} \end{align}

with the descending current slowed by the tangential stress at the interface of the ascending internal plug. Equations (2.20) are valid only for the $\delta$ values that make the average speed for the descending and ascending current negative and positive, respectively.

The other limit of interest $\mathcal {M}\rightarrow 0$ corresponds to a rigid plug, and the speeds are

(2.21)\begin{align} \left. \begin{array}{l@{}} U_d(\delta,\mathcal{M}\rightarrow 0,{Bm}=0,n)={-} \dfrac{\delta^4(\delta^2-1-\ln\delta-\delta^2\ln\delta)}{2(1-\delta^4)},\\ U_a(\delta,\mathcal{M}\rightarrow 0,{Bm}=0,n)= \dfrac{\delta^2(\delta^2-1-\ln\delta-\delta^2\ln\delta)}{2(\delta^2+1)}, \end{array}\right\} \end{align}

again independent on the fluid behaviour index $n$. Note that, for $\mathcal {M} \to 0$, the limiting radius of the plug tends to infinity unless $Bm = 0$ since $r_y=2 Bm /(\mathcal {M}(1-\mathcal {P}))$ with $\mathcal {P}< 1$, so the latter limit condition is only valid for an OdW fluid with $Bm=0$.

For $\mathcal {M}=O(1)$ and $Bm$ sufficiently large, the ascending current is all plug, with $\mathcal {P}\to 2\delta ^2/(1+\delta ^2)$. The condition $r_y=\delta$ is reached for ${Bm}=\mathcal {M}\delta (1-\delta ^2)/(2+2\delta ^2)$, with a maximum value ${Bm}\approx 0.15\mathcal {M}$, so for $Bm$ not particularly large. The asymptotic speeds are again those in (2.21).

A more in-depth analysis of the asymptotic ascending speed for $\mathcal {M} \to 0$ and $\mathcal {M} \to \infty$, and of the corresponding flow rate, is reported in § 5.

3. The Ellis model for the inner ascending current

In the previous section, the HB model adopted for the ascending fluid reduces to the OdW model for $Bm=0$. Now the OdW model typically shows an inconsistency: the apparent viscosity tends to infinity for shear rates tending to zero (see e.g. Myers Reference Myers2005). In order to eliminate this inconsistency and to adapt the rheological model to the effective experimental behaviour of shear-thinning fluids, we adopt the three-parameter Ellis model, which has been extensively and successfully applied to describe the flow in complex geometries (Al-Behadili et al. Reference Al-Behadili, Sellier, Hewett, Nokes and Moyers-Gonzalez2019; Ali et al. Reference Ali, Hussain, Ullah and Bég2019; Celli, Barletta & Brandão Reference Celli, Barletta and Brandão2021; Ciriello et al. Reference Ciriello, Lenci, Longo and Di Federico2021; Picchi et al. Reference Picchi, Ullmann, Brauner and Poesio2021). The Ellis model relates the apparent viscosity to the stress tensor as

(3.1)\begin{equation} \eta=\dfrac{\eta_0}{1+\left(\dfrac{{{{\boldsymbol{\tau: \tau}}}}}{2\tau_0^2}\right)^{(\alpha-1)/2}}, \end{equation}

where $\tau _0$ is the shear stress corresponding to an apparent viscosity $\eta _0/2$ and $\alpha \ge 1$ for a shear-thinning fluid is an indicial parameter. The apparent viscosity tends to 0 for $\tau _0\to 0$ and $\alpha >1$ and tends to $\eta _0/2$ for $\alpha \to 1$. For intense shear stress $\tau \gg \tau _0$ the apparent viscosity reduces to

(3.2)\begin{equation} \eta=\eta_0\left(\dfrac{\tau}{\tau_0}\right)^{1-\alpha}\quad \alpha\ne 1, \end{equation}

and the fluid behaves as an OdW power law with $\alpha =1/n$ and a consistency factor $\mu _0=\eta _0^n\tau _0^{1-n}$. For $\alpha =1$ the fluid is Newtonian with $\eta =\eta _0/2$.

With the same approximations adopted for an ascending HB fluid, the balance of linear momentum in the vertical direction reads

(3.3)\begin{equation} \dfrac{1}{\mathcal{M'}}\dfrac{\partial \tilde{u}_a}{\partial \tilde{r}}=\left(\underbrace{\dfrac{\partial {\tilde{p}_a}}{\partial {\tilde{z}}} + \dfrac{\mathcal{R}}{1-\mathcal{R}}}_{\mathcal{P}-1}\right) \left[1+\beta\left|\underbrace{\dfrac{\partial {\tilde{p}_a}}{\partial {\tilde{z}}} +\dfrac{\mathcal{R}}{1-\mathcal{R}}}_{\mathcal{P}-1}\right|^{\alpha-1} \left(\dfrac{\tilde{r}}{2}\right)^{\alpha-1}\right] \dfrac{\tilde{r}}{2},\quad \tilde{r}\in[0,\tilde{\delta}], \end{equation}

where a second viscosity ratio

(3.4)\begin{equation} \mathcal{M'}=\dfrac{\mu_d}{\eta_0},\end{equation}

was introduced, and the parameter

(3.5)\begin{equation} \beta=\left(\dfrac{{\rm \Delta} \rho\,g R}{\tau_0}\right)^{\alpha-1},\end{equation}

measures the relative importance of buoyancy and the typical shear stress of an Ellis fluid. The boundary conditions are

(3.6)\begin{equation} \left. \begin{array}{c@{}} u_d(1)=0,\\ u_a(0)\ \textrm{is finite},\\ u_d(\delta)=u_a(\delta),\\ \left.\dfrac{\textrm{d}u_a}{\textrm{d} r}\right|_{r=\delta} =\mathcal{M}'\left[1+\beta\left|\left.\dfrac{\textrm{d}u_d}{\textrm{d} r}\right|_{r=\delta}\right|^{\alpha-1}\right] \left.\dfrac{\textrm{d}u_d}{\textrm{d} r}\right|_{r=\delta}, \end{array}\right\} \end{equation}

with the same meaning as (2.10) and where the tilde has been dropped.

Solving the differential problem, we calculate the following velocity profiles:

(3.7)\begin{equation} \left. \begin{array}{ll@{}} u_d(r)=\dfrac{\mathcal{P}}{4}(r^2-1)-\dfrac{\delta ^2}{2}\ln{r} , & {r}\in[{\delta},1],\\ \begin{aligned} u_a(r) & =\mathcal{M}'(1-\mathcal{P})\dfrac{\delta^2-r^2}{4}+\mathcal{M}' \beta(1-\mathcal{P})^\alpha\dfrac{\delta^{\alpha+1}-r^{\alpha+1}}{2^\alpha (\alpha+1)} \\ & \quad+\dfrac{\mathcal{P}}{4}(\delta^2-1)- \dfrac{\delta^2}{2}\ln \delta,\end{aligned} & {r}\in[0,{\delta}],\\ \end{array}\right\} \end{equation}

where we have assumed $\mathcal {P}<1$ for physical consistency. For $\tau _{rz}\gg \tau _0$ (3.7) simplify as follows:

(3.8)\begin{equation} \left. \begin{array}{ll@{}} u_d(r)=\dfrac{\mathcal{P}}{4}(r^2-1)-\dfrac{\delta ^2}{2}\ln{r} , & {r}\in[{\delta},1],\\ \begin{aligned} u_a(r) & \approx\mathcal{M}'\beta(1-\mathcal{P})^\alpha \dfrac{\delta^{\alpha+1}-r^{\alpha+1}}{2^\alpha(\alpha+1)} \\ & \quad+\dfrac{\mathcal{P}}{4}(\delta^2-1)- \dfrac{\delta^2}{2}\ln \delta, \end{aligned} & {r}\in[0,{\delta}], \end{array}\right\} \end{equation}

which upon substituting $\alpha =1/n$ and $\mu _0=\eta _0^n\tau _0^{1-n}$ collapse to (2.13) valid for an OdW fluid.

Figure 9 shows the velocity profiles computed for an ascending Ellis fluid, for different values of $\mathcal {M}',\beta,\alpha$ and for $\delta =0.6$. The profiles are qualitatively similar to those for an OdW fluid, with a backflow of the inner current due to the drag of the outer current; the backflow is accentuated as the buoyancy dominates over the typical shear stress of an Ellis fluid (larger $\beta$). More shear-thinning fluids (larger $\alpha$) exhibit flatter velocity profiles and hence less backflow. An increasing viscosity contrast $\mathcal {M}'$ reduces the differences between velocity profiles for different values of $\beta$ and $\alpha$.

Figure 9. Ellis ascending fluid and Newtonian descending fluid. Velocity profiles for a generic configuration with $\delta =0.6$, different values of the ratio $\beta, \alpha =2,3$ and (a) $\mathcal {M}'=1$ and (b) $\mathcal {M}'=10$.

The average speeds of the descending outer current and of the ascending inner current are

(3.9)\begin{equation} \left. \begin{array}{c@{}} U_d=\dfrac{1}{8(1-\delta^2)}[2\delta^2(1-\delta^2+2\delta^2\ln\delta)- (1-\delta^2)^2\mathcal{P}],\\ U_a=\dfrac{ \mathcal{M}' \beta\delta ^{\alpha +1} (1- \mathcal{P})^{\alpha }}{2^{\alpha }(\alpha +3)}+\dfrac{1}{8} \mathcal{M}'\delta ^2 (1-\mathcal{P})-\dfrac{1}{4} [2 \delta ^2 \ln\delta+ \mathcal{P}(1-\delta^2)], \end{array}\right\} \end{equation}

and the corresponding flow rates are $Q_d={\rm \pi} (1-\delta ^2)U_d, Q_a={\rm \pi} \delta ^2U_a$.

Figure 10 shows the average speed as a function of the long drop radius $\delta$, of the ratio of scale stresses $\beta$, of the viscosity ratio $\mathcal {M}'$ for Ellis shear-thinning fluids with $\alpha =2, 3$ and $\beta =10,100,1000,10\,000$. As for OdW and HB fluids, the ascent speed peaks at an intermediate value of the long drop radius; further, it is strongly impacted by the $\beta$ value, and significantly decreases with increasing shear-thinning behaviour, especially in the low range of $\delta$. A similar behaviour is shared by the Newtonian descent speed, except the absolute values of the speed are lower.

Figure 10. Long drop speeds for an Ellis ascending fluid and a Newtonian descending fluid. (a) Ascent speed and (b) descent speed as a function of the long drop radius $\delta$ and the ratio $\beta$ for $\alpha =2,3, \mathcal {M}'=10$.

By imposing $Q_a+Q_d=0$, an equation in $\mathcal {P},\mathcal {M}',\beta,\alpha,\delta$ is obtained, which may be used to calculate the driving term $\mathcal {P}$ for given values of the other parameters. The equation admits analytical solutions for $\alpha =1,2,3,4$, and requires a numerical solution for other values of $\alpha$.

The asymptotic velocities for $\mathcal {M}'\to \infty$ and for $\mathcal {M}'\to 0$ collapse to (2.19) and to (2.21), respectively, and are independent of $\beta$.

A more in-depth analysis of the asymptotic ascending speed for $\mathcal {M}' \to 0$ and $\mathcal {M}' \to \infty$, and of the corresponding flow rate, is reported in § 5.

4. The energy balance for a Newtonian descending fluid and a HB or Ellis ascending fluid

The present section formulates the energy balance for the non-Newtonian ascending Taylor drop; the balance is valid after enough time has elapsed from its formation, allowing the drop to elongate so that the validity of the lubrication approximation previously employed is ensured. The ascending drop gains energy in potential form, loses energy due to viscous dissipation and stores some energy as kinetic energy of the two counter-current fluids. Neglecting the work performed by the forces at the interfaces, the balance is between kinetic energy $\mathcal {E}_k$, potential energy $\mathcal {E}_p$ and viscous dissipation $\varPhi$ according to (Picchi et al. Reference Picchi, Suckale and Battiato2020)

(4.1)\begin{equation} \dfrac{\textrm{d}\mathcal{E}_k}{\textrm{d}t}+ \dfrac{\textrm{d}\mathcal{E}_p}{\textrm{d}t}+\varPhi=0,\end{equation}

written in dimensional form. The corresponding scales are $\mathcal {E}_k \sim \rho _d U^2 \mathcal {L}R^2, \mathcal {E}_p\sim (\rho _d-\rho _a) g\mathcal {L}^2R^2$ and $\varPhi \sim \mu _dU^2\mathcal {L}$, hence (4.1) can be written in dimensionless variables as

(4.2)\begin{equation} \dfrac{1}{{Ri}(1-\mathcal{R})}\dfrac{\textrm{d} \tilde{\mathcal{E}}_k}{\textrm{d}\tilde{t}}+\dfrac{1}{1-\mathcal{R}} \dfrac{\textrm{d}\tilde{\mathcal{E}}_p}{\textrm{d}\tilde{t}} +\tilde{\varPhi}=0,\quad {Ri}=1/[({1-\mathcal{R}})\varepsilon {Ar}], \end{equation}

where ${Ri}$ is the Richardson number expressing the ratio between potential and kinetic energy. For the special case ${Ri}\to \infty$ considered hereinafter the kinetic energy variation is negligible and (4.2) reduces to

(4.3)\begin{equation} \dfrac{\textrm{d}\tilde{\mathcal{E}}_p}{\textrm{d}\tilde{t}}={-}(1-\mathcal{R})\tilde{\varPhi}.\end{equation}

The potential energy in (4.1) is computed by integration over the drop region and is equal to

(4.4)\begin{align} \mathcal{E}_p&=\rho_dg{\rm \pi}\left[(R^2-\delta^2)\int_{z_d}^{z_a}z\,\textrm{d} z+R^2\int_{z_a}^{L}z\,\textrm{d} z\right]\nonumber\\ &\quad +\rho_ag{\rm \pi}\left[\delta^2\int_{z_d}^{z_a}z\,\textrm{d} z+R^2\int^{z_d}_{{-}L}z\,\textrm{d} z\right]. \end{align}

In dimensionless form the potential energy results in

(4.5)\begin{equation} \tilde{\mathcal{E}}_p=\tfrac{1}{2}{\rm \pi}[(1-\tilde{\delta}^2) \tilde{U}_d^2-\tilde{\delta}^2\tilde{U}_a^2] \tilde{t}^{\,2}+\textrm{const.}, \end{equation}

or equivalently

(4.6)\begin{equation} \mathcal{E}_p=\dfrac{2\delta^2-1}{2{\rm \pi}\delta^2(1-\delta^2)} Q_a^2t^2+\textrm{const.}, \end{equation}

where $Q_a$ is the flow rate of the ascending current, also defined as the transport number $Te$ (see Huppert & Hallworth Reference Huppert and Hallworth2007), and the tilde symbol is omitted for simplicity.

The dissipation rate is computed as

(4.7)\begin{equation} \varPhi=\int_{\mathcal{V}_d}({{\boldsymbol{\tau:}\boldsymbol{\mathsf{D}}}})_d\, \textrm{d}\mathcal{V}+\int_{\mathcal{V}_a}({{\boldsymbol{\tau:}\boldsymbol{\mathsf{D}}}})_a \,\textrm{d}\mathcal{V}, \end{equation}

where $\mathcal {V}_{a,d}$ refers to the volume of ascending lighter/descending denser fluid.

In a cylindrical geometry and under a HB rheology for the internal current one has

(4.8)\begin{equation} \varPhi\equiv\varPhi_d+\varPhi_a=2{\rm \pi}\mathcal{L}(t)\left\{\mu_d \int_{\delta}^R \left(\dfrac{\textrm{d}u_d}{\textrm{d} r}\right)^2 \,\textrm{d}r +\int_{0}^\delta \left[\mu_0\left|\dfrac{\textrm{d} u_a}{\textrm{d} r}\right|^{n+1}+\tau_y\left|\dfrac{\textrm{d} u_a}{\textrm{d} r}\right|\right]\,\textrm{d}r\right\}, \end{equation}

or in dimensionless form

(4.9)\begin{equation} \tilde{\varPhi}\equiv\tilde{\varPhi}_d+\tilde{\varPhi}_a= 2{\rm \pi}\tilde{\mathcal{L}}(\tilde{t})\left\{\int_{\tilde{\delta}}^1 \left(\dfrac{\textrm{d}\tilde{u}_d}{\textrm{d} \tilde{r}} \right)^2\textrm{d}\tilde{r} +\dfrac{1}{\mathcal{M}} \int_{0}^{\tilde{\delta}} \left[\left|\dfrac{\textrm{d} \tilde{u}_a}{\textrm{d} \tilde{r}}\right|^{n+1}+{Bm} \left|\dfrac{\textrm{d}\tilde{u}_a}{\textrm{d} \tilde{r}} \right|\right]\,\textrm{d}\tilde{r}\right\}.\end{equation}

Upon integration using (2.11), (4.9) transforms into

(4.10)\begin{equation} \left. \begin{array}{c} \varPhi_d = 2{\rm \pi}(U_a-U_d)\left\{ \dfrac{\rm \pi}{8}\mathcal{P} [(1-\delta^4)\mathcal{P}-4\delta^2(1-\delta^2)]- \dfrac{\rm \pi}{2}\delta^4\ln\delta\right\}t, \\ \begin{aligned} \varPhi_a & = 2{\rm \pi}(U_a-U_d)\left\{\dfrac{{\rm \pi} n[\delta \mathcal{M}(1- \mathcal{P})-2 Bm ]^{(2n+1)/n} [2n Bm +\delta \mathcal{M} (2 n+1)(1-\mathcal{P})]}{ 2^{1/n}\mathcal{M}^3 (6 n^2+5 n+1) (1-\mathcal{P})^2}\right.\\ & \quad+\left.\dfrac{ {\rm \pi} Bm [ Bm +\delta \mathcal{M}(1-\mathcal{P})] [2 Bm -\delta \mathcal{M}(1-\mathcal{P})]^2}{3 \mathcal{M}^3 (1-\mathcal{P})^2}\right\}t, \end{aligned} \end{array}\right\} \end{equation}

where the symbol $\widetilde {\cdots }$ is omitted for simplicity, and where $\mathcal {L}(t)=(U_a-U_d)t$ is the time varying length of the drop. This suggests an alternative formulation of (4.9) as

(4.11)\begin{equation} \varPhi=Q_a(U_a-U_d)t\equiv \dfrac{Q_a^2}{{\rm \pi}\delta^2(1-\delta^2)}t, \end{equation}

clarifying that, in dimensionless form, the dissipation rate per unit length is equal to the flow rate. Figure 11 depicts the dissipation rate per unit length $\varPhi /\mathcal {L}$ for two different viscosity ratios and several combinations of ascending fluids. The dissipation rate is zero at $\delta =0,1$ and peaks at an intermediate value of $\delta$ in all cases. A higher total dissipation rate is attained for a larger viscosity ratio, while shear-thinning effects decrease or increase the total dissipation rate with respect to the Newtonian case $n=1$ depending on the value of the viscosity ratio. The appearance of a non-zero yield stress induces a lower total dissipation rate for many combinations of viscosity ratio and fluid behaviour index, while for other combinations the total dissipation rate is practically unchanged. There exists a complex interplay between shear rate and shear stress, so it cannot be immediately established which of the two terms is more affected by the presence of yield stress (non-zero $Bm$) and by the value of the viscosity ratio.

Figure 11. Dissipation rate per unit length $\varPhi /\mathcal {L}$ as a function of $\delta$ for $\mathcal {M}=1,10^3, n=0.5,1.0$ and ${Bm}=0.0,0.1, 1.0$. The continuous curves refer to an OdW ascending fluid ($Bm = 0.0$), the dashed curves to $Bm = 0.1$ and the dot-dashed curves to $Bm = 1.0$ The red curves refer to the ascending HB or OdW current, the blue curves to the descending Newtonian current, the black curves are the total dissipation. For $\mathcal {M}=1$ only curves for $Bm = 0.1$ can be drawn, since for $Bm = 1.0$ it is always $\delta < r_y$.

Looking now at the individual contributions of the ascending and descending currents to the total dissipation rate, it is seen that the contribution of the ascending current depends on the viscosity ratio, fluid behaviour index and modified Bingham number, but tends to be smaller than that of the descending Newtonian current in most cases.

The maximum dissipation rate is attained for intermediate values of $\delta$ in all cases. Figure 12 shows the maximum dissipation rate and the corresponding radius of the ascending internal current, $\delta _{{max}}$, vs $\mathcal {M}$. The maximum dissipation rate increases with the viscosity contrast, while the opposite is true for the corresponding radius. Both results moderately depend on the fluid behaviour index, with values in a range limited by the condition $\delta >r_y$. A modest dependence on the modified Bingham number is evident only for the maximum dissipation rate and not for the corresponding radius.

Figure 12. Maximum dissipation rate (coincident with maximum flow rate) per unit length $\varPhi /\mathcal {L}$ and corresponding radius of the internal ascending current as a function of $\mathcal {M}$. The continuous curves refer to an OdW ascending fluid ($Bm = 0.0$), the dashed and dot-dashed curves refer to a HB ascending fluid with $Bm = 0.1, 1.0$, respectively.

We have now the expressions of the terms involved in the energy balance, (4.3), and since $\mathcal {E}_p=f(\delta (t),n,\mathcal {M},{Bm},t)$, (4.3) can be written as

(4.12)\begin{equation} \dfrac{\partial \mathcal{E}_p}{\partial \delta}\dfrac{\textrm{d} \delta}{\textrm{d}t}+\dfrac{\partial \mathcal{E}_p}{\partial t}={-}(1-\mathcal{R})\varPhi, \end{equation}

or

(4.13)\begin{align} \dfrac{\textrm{d}\delta}{\textrm{d}t}=\dfrac{-(1-\mathcal{R})\varPhi- \dfrac{\partial \mathcal{E}_p}{\partial t}}{\dfrac{\partial \mathcal{E}_p}{\partial \delta}}\equiv{-}\dfrac{\delta(1-\delta^2)(\mathcal{R}-2\delta^2)}{-(1- 2\delta^2+2\delta^4)+\delta(1-3\delta^2+2\delta^4)\dfrac{\partial \ln Q_a}{\partial \delta}}\dfrac{1}{t},\end{align}

which describes the time evolution of the radius of the ascending current, impacted by rheology via the flow rate $Q_a$. Figure 13 shows the evolution of the current by integrating (4.13) with different initial conditions for a given OdW fluid. For $\mathcal {M}=1$, initial conditions within the yellow band evolve into the asymptotic value $\delta _{3\infty }=\sqrt {2\mathcal {R}}/2$; slightly different results are obtained for finite values of $\mathcal {M}$, including $\mathcal {M}\to 0$; the two physically meaningless asymptotes $\delta _{1\infty }=0$ and $\delta _{2\infty }=1$ are reached for initial conditions out of the yellow band. For other $\mathcal {M}$ values, the width of the band varies. For $\mathcal {M}\to \infty$ the asymptotic value $\delta _{1\infty }=0$ is missing, and the domain of attraction of $\delta _{3\infty }$ is larger. The results for a HB ascending fluid (not shown) are very similar, except for the case $\mathcal {M}\to \infty$, which implies again the existence of the asymptotic solution $\delta _{1\infty }$. Note that, for a HB ascending fluid, the model requires that $\delta >r_y$.

Figure 13. Time evolution of the radius of the internal ascending current for $n=0.5, \mathcal {R}=0.8, {Bm} = 0.0$ and (a) $\mathcal {M}=1$, (b) $\mathcal {M}\to \infty$.

Figure 14 depicts the time derivative of the radius of the ascending current, the right-hand side of (4.13), at a given time for $\delta \in [0,1]$ and different values of $\mathcal {M}$, for a shear-thinning OdW or HB fluid. For finite values of $\mathcal {M}$ there are three stable nodes for $\delta _{1\infty }=0, \delta _{2\infty }=1$ and $\delta _{3\infty }=\sqrt {2\mathcal {R}}/2$, albeit the first two are physical meaningless. It is noteworthy that similar results were obtained by Picchi et al. (Reference Picchi, Suckale and Battiato2020) for Newtonian ascending fluids. The domain of attraction of the stable node $\delta _{3\infty }$ depends mainly on $\mathcal {M}$, see figure 15 where the domain of attraction is shown for different values of the fluid behaviour index $n$ and for an OdW or HB fluid. The shear-thinning behaviour favours moderately larger bandwidths, whereas the effect of the yield stress is generally modest.

Figure 14. Evolution of the radius of the internal ascending current at $t=0.1$ for $n=0.5, \mathcal {R}=0.8, {Bm} = 0$, and (a) $\mathcal {M}=0$, (b) $\mathcal {M}=1$, (c) $\mathcal {M}\to \infty$. The grey curves refer to a HB fluid with ${Bm}=0.1$, and can be drawn only in the domain $\delta >r_y$. The latter condition is never satisfied for $\mathcal {M}=0$ (a), is always satisfied for $\mathcal {M}\to \infty$ (b) and is satisfied in a limited range of $\delta$ for $\mathcal {M}=1$ and ${Bm}=0.1$ (c). The red curves refer to a HB fluid with ${Bm}=1$, and can be drawn only in the domain $\delta >r_y$, which is empty for $\mathcal {M}=0, 1$ and is non-empty only for large $\mathcal {M}$.

Figure 15. Domain of attraction of the stable node $\delta _{3\infty }=\sqrt {2\mathcal {R}}/2$ for $n=1,1/2,1/3, \mathcal {R}=0.8, {Bm} = 0$ (OdW fluid). The grey curves, more evident in the enlargement below, refer to a HB fluid with ${Bm} = 0.1$ (thick grey) and to ${Bm} = 1.0$ (thin grey), the blue dashed and dot-dashed curves represent the condition $\delta >r_y$ for ${Bm} = 0.1, 1.0$, respectively. The vertical arrows indicate the direction of evolution over time of the radius of the ascending current towards the stable node $\delta _{3\infty }=\sqrt {2\mathcal {R}}/2$.

If the ascending fluid is modelled according to the Ellis model, the dimensionless dissipation rate is computed as

(4.14)\begin{align} \varPhi\equiv\varPhi_d+\varPhi_a=\mathcal{L}(t)\left[\int_{\delta}^1 2{\rm \pi}\left(\dfrac{\textrm{d}u_d}{\textrm{d} r}\right)^2r\, \textrm{d}r+\int_{0}^\delta 2{\rm \pi}\mathcal{M}'(\tau_{rz}^2+ \beta|\tau_{rz}|^{\alpha+1})r\,\textrm{d}r\right], \end{align}

resulting in

(4.15)\begin{equation} \left. \begin{array}{c@{}} \varPhi_d = 2{\rm \pi}(U_a-U_d)\left\{ \dfrac{\rm \pi}{8}\mathcal{P} [(1-\delta^4)\mathcal{P}-4\delta^2(1-\delta^2)]- \dfrac{\rm \pi}{2}\delta^4\ln\delta\right\}t, \\ \varPhi_a= 2{\rm \pi}(U_a-U_d) \left[\dfrac{8 \beta \delta ^{3+ \alpha} }{\alpha +3}| 1-\mathcal{P}|^{\alpha +1}+2^{\alpha } \delta^4 (1-\mathcal{P})^2\right]\dfrac{{\rm \pi}\mathcal{M}'}{2^{\alpha +3}}t. \end{array}\right\} \end{equation}

The dimensionless dissipation rate is again expressed as a function of the flow rate of the ascending current, see (4.11). The time evolution of the radius of the internal current is identical to (4.13), only the expression of $Q_a$ changes as a consequence of the different rheological model. When compared, the results are similar to those obtained for the subcase of power-law fluid (not shown). It is important to point out that the asymptotic value $\delta _{3\infty }$ does not depend on the rheological behaviour of the ascending fluid, but only on $\mathcal {R}$.

5. The asymptotic ascending speed and flow rate

Adopting the HB model or the Ellis model for the inner ascending current, it is seen that its ascent speed has two plateaus for $\mathcal {M}\,(\mathcal {M}')\to 0$ and for $\mathcal {M}\,(\mathcal {M}')\to \infty$, where it becomes only a function of its radius. Figures 16 and 17 show the curves corresponding to zero and non-zero yield stress for different values of the fluid behaviour index and two different ratios of buoyancy to Ellis shear stress $\beta$ for different values of the indicial parameter $\alpha$, respectively. For both rheologies, the asymptotic value is reached more quickly for less shear-thinning behaviour, without yield stress (${Bm}=0)$ for the HB fluid and larger $\beta$ for the Ellis fluid. The insets show the values of $\mathcal {M},\mathcal {M}'$ above which 95 % of the asymptotic ascent speed is reached, as a function of $n,{Bm}$ or $\alpha,\beta$. The asymptotic values of the viscosity ratio are an increasing function of ${Bm}$ and a decreasing function of $\beta$. For $\mathcal {M}>150$ or $\mathcal {M}'>50$, one can assume the asymptote is always reached irrespective of the other parameters (except $\mathcal {R}$ that is included in $\delta _{3\infty }$); by interpolating the numerical values of the asymptotic velocity computed for several combinations $\mathcal {R}-\mathcal {M}$ in the range $0.5<\mathcal {R}<1$ and $\mathcal {M}>150$ (or $\mathcal {M}'>50$), the ascent speed and the corresponding flow rate can be approximated as

(5.1)\begin{equation} \left. \begin{array}{@{}l} U_a\approx 1.34\exp({-}4\mathcal{R})\\ Q_a\equiv Te \approx2.1\mathcal{R}\exp({-}4\mathcal{R}) \end{array}\quad \textrm{for}\ \mathcal{M}>150(\mathcal{M}'>50) \ \textrm{and}\ 0.5<\mathcal{R}<1,\right\}\end{equation}

where a negative exponential function of the density ratio $\mathcal {R}$ was adopted for interpolation. Equation (5.1) becomes in dimensional variables

(5.2)\begin{align} \left. \begin{array}{@{}l} {U}_a\approx 1.34\exp({-}4\mathcal{R})\dfrac{{\rm \Delta} \rho\,gR^2}{\mu_d}\\ {Q}_a\approx 2.1\mathcal{R}\exp({-}4\mathcal{R}) \dfrac{{\rm \Delta} \rho\,gR^4}{\mu_d}\end{array} \quad \textrm{for}\ \mathcal{M}>150\ \textrm{or} \mathcal{M}'>50\ \textrm{and}\ 0.5<\mathcal{R}<1. \right\}\end{align}

The transport number $Te$ in (5.1) is consistent with the transport number computed by adopting (4.7) in Picchi et al. (Reference Picchi, Suckale and Battiato2020), also reproducing the results by Wallis (Reference Wallis1969) and Brauner & Ullmann (Reference Brauner and Ullmann2004) for Taylor bubbles in the viscosity-controlled regime; it is generally smaller than $Te =0.125\ \textrm {for}\ \mathcal {M}\gg 1$ in Huppert & Hallworth (Reference Huppert and Hallworth2007), where it was computed with reference to the condition of maximum flow rate.

Figure 16. The HB model for the ascending fluid, ascent speed $U_a$ as a function of $\mathcal {M}$ for fluid behaviour index $n=0.2, 0.4, 0.6, 0.8, 1.0$. Bold curves show the ascent speed corresponding to the asymptotic core radius $\delta _{3\infty }=\sqrt {2}/2$ and ${Bm}=0$; dashed and dot-dashed curves refer to ${Bm}=0.10, 1.0$, respectively. The inset shows the values of $\mathcal {M}$ as a function of $n$ above which 95 % of the asymptotic ascent speed is reached.

Figure 17. As in figure 16 but with the Ellis model describing the internal ascending fluid: ascent speed $U_a$ as a function of $\mathcal {M}'$ for an indicial parameter $\alpha =1,1.5,2,2.5,3$. Bold curves refer to a ratio of buoyancy to Ellis shear stress $\beta =100$, dashed curves to $\beta =10$.

6. The experiments

6.2. Experimental set-up and protocol

To validate the model and verify the assumptions and limitations, we conducted a series of tests in the Hydraulic Laboratory of the University of Parma. The experimental apparatus consists of a polymethyl methacrylate (PMMA) pipe with an internal diameter of $14.1 \pm 0.2\ \textrm {mm}$ (calculated using a volumetric method) and a length of 150 cm, see figure 18(a). The pipe is closed at both ends and, in the middle, has a slide gate that allows the two 75 cm portions of the pipe to be separated or connected. In a cross-section 5 cm from the gate, a flat-walled cell of PMMA was fitted around the pipe, to be filled with glycerol, which has a refractive index 1.474 at $20\,^{\circ }\textrm {C}$ almost identical to that of PMMA; the refractive index of aqueous glycerol solutions increases linearly with increase in glycerol concentration and can be expressed as (Takamura, Fischer & Morrow Reference Takamura, Fischer and Morrow2012) $n_r=1.333w_w+1.474w_{gl}$, where $w_w$ and $w_{gl}$ are the weight fraction of water and glycerol in the mixture. Since the tests were carried out using glycerol as a descending fluid, always adhering to the walls of the pipe in the experimental configuration adopted, the experimental apparatus allows correction of the refraction with a negligible residual distortion. Some tests were conducted with honey as a descending fluid, which has an index of refraction almost identical to glycerol. To check the refractive index correction, a needle with marks engraved every half a millimetre was inserted diametrically into the pipe, see figure 18(b): the marks look equispaced even in the presence of the circular walls of the pipe.

Figure 18. Experimental set-up. (a) Vertical pipe with the USB microscope camera and the video camera for large-scale image analysis; (b) a photo of the pipe, as seen from the USB camera microscope, containing glycerol and inserted in the box filled with glycerol, in order to correct the image distortion. The needle with a series of equispaced marks shows the efficiency of the distortion correction. See also the enlargement.

The verticality of the pipe was controlled with an electronic spirit level accurate to one tenth of a sexagesimal degree; the temperature of the fluids, in thermal equilibrium with the environment, was measured with a thermometer accurate to a tenth of a degree Celsius.

The fluid of the outer current, which is heavier and more viscous than the inner current, is glycerol as anticipated; the fluid of the inner current is water (Newtonian), water and Carboxymethyl cellulose and Xanthan Gum (both shear thinning), water and Carbopol (a polyacrylic acid) neutralized with NaOH (shear thinning with yield stress, non-thixotropic). Rheometric measurements were carried out with an Anton Paar TwinDrive MCR 702 parallel plate rheometer and with a Haake Rotovisco RT10 with Searle coaxial cylinder measuring system, in shear rate control; density measurements were performed with a pycnometer and an Anton Paar DMA 5000 M densimeter.

Images were taken in the refraction corrected section with a microscopy webcam with $640\times 480$ pixels at 25 frames per second (f.p.s.), with a field of view approximately $30\times 25\ \textrm {mm}^2$ and with a resolution $\approx 0.05\ \textrm {mm}\ \textrm {pixel}^{-1}$. The inner current was coloured with aniline dye to facilitate image analysis. The diameter of the inner current was estimated using proprietary software in Matlab that involves manual identification of contours and conversion from pixels to metric coordinates. The operation is repeated in ten cross-sections for ten frames, chosen in the absence of evident perturbations. The interface between inner and outer currents is not sharp, and a film separating the two fluids is evident in the shots; measurements were taken considering (i) the interface between film and inner current, and (ii) the interface between film and outer current.

The average ascent speed of the internal lighter current was measured by image analysis taken from an iPhone at 30 f.p.s., allowing estimation at different times of the position of the upwelling front with respect to an external rule with spacing of 1 mm.

The experiments were conducted by first filling the bottom half of the tube with the lighter fluid, closing the gate and then filling the top half of the tube with the heavier fluid. When the gate is opened, the upward current, the innermost one for all the tests carried out, starts to advance.

Figure 19 shows three snapshots of the elongated drop of water in glycerol.

Figure 19. Experiment 2 (see table 1) with inner ascending and outer descending Newtonian fluids. Here, $\mathcal {R}=0.796, {Ar}=1.11, \mathcal {M}=10^3, {Bm}=0$. Three snapshots at different stages of the current evolution are shown.

Table 1. Parameters of the experiments with a HB model for the ascending fluid. Here, $\rho _{a,d}$ is the density of the ascending/descending fluid, $\mathcal {R}=\rho _a/\rho _d$, $Ar$ is the Archimedes number, $T$ is the temperature of the fluids during the test, $\mu _d$ is the viscosity of the descending fluid, $n,\mu _0,\tau _y$ are the fluid behaviour index, the consistency index and the yield stress of the ascending fluid, $\mathcal {M}$ and $Bm$ are the two dimensionless groups defined in (2.6) and (2.7), ${Re}=\rho _dR\sqrt {{\rm \Delta} \rho g R}/\mu _d$ is the Reynolds number and $\delta _{exp}$ and $U_{a-exp}$ are the experimental values of the core radius and of the ascent speed, respectively. The last column indicates the composition of the inner ascending fluid: W stands for water, Cb stands for Carbopol, A stands for air, XG stands for Xanthan Gum and CM stands for Carboxymethyl cellulose. The experiments with a number followed by a star have also been interpreted with an Ellis model for the inner ascending fluid, see table 2. The descending fluid is glycerol except for experiments 6 and 11, where honey was used. The symbol $\clubsuit$ indicates that a video is available as supplementary material.

For some tests, a velocity profiler (UVP) DOP 5000 by Signal Processing SA, Switzerland, was used, with an 8 MHz carrier frequency probe inclined at $75^\circ$ to the vertical. The fluid was seeded with $\textrm {TiO}_2$ particles, showing a high sonic impedance, so that the probe could measure their velocity at different distances (gates) along the axis of the ultrasonic cone on the basis of the Doppler shift of the echoes. Doppler shift and gate position depends on the speed of sound, which varies as a function of the temperature and of the density of the fluid. A correction to the gate position and velocity was made in order to compensate for the different velocities of the Ultrasounds in the inner and outer currents (see Petrolo & Longo Reference Petrolo and Longo2020). The space resolution was 0.3 mm, the data rate was $\approx 5$ profiles per second, so that the overall accuracy in fluid velocity measurement is $\le 4\,\%$. Figure 20 shows two velocity profiles for tests with a Newtonian and a shear-thinning ascending fluid, respectively, and a third test with a stream of water with known flow rate in order to compare measurements with the theoretical parabolic velocity profile. The overall accuracy is fairly good, although the finite size of the volume of measurement of the UVP smooths the velocity. An improved agreement between experiments and theory is obtained by considering the space average of the theoretical velocity over the thickness of the volume of measurements (dashed curves). In a first attempt, we used the UVP data to estimate the zero-crossing distance of the velocity profile, related to the radius of the internal current. However, the limited spatial resolution of the instrument did not guarantee adequate accuracy, hence we decided to measure the radius optically with the webcam microscope.

Figure 20. Velocity profiles measured with ultrasonic profiler, average values over approximately 20 s. (a) Ascending Newtonian fluid, $n=1, \mathcal {R}=0.80, \mathcal {M}=640, {Ar}=1.88$; (b) ascending shear-thinning fluid, $n=0.55, \mathcal {R}=0.83, \mathcal {M}=5.7, {Ar}=2.63$. (c) Water stream in laminar viscous regime. The blue bold curves are the theoretical values, the red dashed curves are the theoretical values averaged in a volume equal to the volume of measurement of the gates of the velocity profiler, symbols are the experimental values, error bars refer to one standard deviation and are representative of the variability of the sample of the velocity profiles. The interface ($r=\delta$) is located at the intersection between the two branches of ascending and descending fluid velocity (the velocity shows a cusp), based on theoretical profiles (blue continuous curves).

6.4. The results of the experiments

A series of 12 experiments were carried out, 5 with Newtonian (2 with air), 3 with shear-thinning and 4 with yield shear-thinning fluids constituting the inner ascending current. The outer descending fluid was always Newtonian, glycerol or honey. The experiments were not repeated and a single realization for each set of fluids properties is documented and discussed. The main parameters of the experiments are listed in table 1, where the OdW model is used for the shear-thinning ascending fluid without yield stress (experiments 1–8). Table 2 lists the parameters for experiments 1–8 with the inner ascending fluid modelled with the Ellis rheology. Measurements with both ascending and descending Newtonian fluid currents are already available in the literature, see Stevenson & Blake (Reference Stevenson and Blake1998) and Goldsmith & Mason (Reference Goldsmith and Mason1962); this choice was initially adopted in the present apparatus in order to verify its correct functioning, as well as to extend the database.

Table 2. Parameters of the experiments with an Ellis model for the ascending fluid. For caption, see table 1. Here, $\alpha, \eta _0$ and $\tau _0$ are the three parameters of the Ellis model, $\mathcal {M}'$ and $\beta$ are two dimensionless groups defined in (3.4) and (3.5). These experiments have also been interpreted with an OdW model for the inner ascending fluid.

A clear rationale for the appropriateness of the Ellis model instead of the OdW model is evident in figure 21(a), which shows the experimental rheometric measurements for the shear-thinning fluid adopted in experiment 6, see table 2, with data interpolated with the Ellis and the OdW models. The Ellis model better interprets the low shear-rate data, and this affect the velocity profile mainly near the axis, where the shear rate is at a minimum. At high shear rate, the two models are almost coincident.

Figure 21. Rheometric measurements (a) for the shear-thinning fluid in experiment 6, see tables 1 and 2. The Ellis model (bold curves) and the OdW model (dashed curves) are used to interpolate the experimental data $\tau -\dot {\gamma }$ (crosses), and $\eta -\dot {\gamma }$ (open squares). (b) For the HB fluid in experiment 12. Data have been decimated for clearer visualization. The grey symbols refer to experimental points not included in the interpolation process due to their limited accuracy.

Figure 21(b) shows the rheometric data for the HB fluid in experiment 12, where the continuous curve has been obtained by interpolating the experimental points. The flow curves for the HB fluids adopted in experiments 9–11 are shown in Appendix C. Note that the determination of rheometric parameters for HB fluids can be based on different criteria, see Di Federico et al. (Reference Di Federico, Longo, King, Chiapponi, Petrolo and Ciriello2017) where some of the most frequently adopted techniques are compared in the appendix. In the present study, we have adopted the simplest procedure, which involves a relatively high uncertainty. A more detailed analysis would also have required checking for the existence of slip phenomena, which lead to an underestimation of the measured parameters. We leave this analysis to further developments of both the adopted model and the rheometric measurements. To some extent, the high uncertainty in parameter estimation also includes some deliberately neglected phenomena.

The experiments cover a wide range of viscosity contrast $\mathcal {M}=O(1)-O(10^6)$ always with a stable ascending current, while the density ratio $\mathcal {R}\approx 0.8$, except for two tests where the density ratio is much smaller. The Reynolds number is ${Re}=O(0.1)$, the Archimedes number is ${Ar}=O(1)$, the modified Bingham number for HB fluids, when non-zero, is ${Bm}=O(0.1)$ and the ratio between buoyancy and the typical shear stress of an Ellis fluid is $\beta =O(1)- O(10^4)$. Note that experiments with even higher $Bm$ values showed clear instabilities, with a rising current fragmented into numerous discontinuous plugs. The maximum and minimum values of the Richardson number ${Ri}$ in the experiments are ${Ri}_{{max}}\approx 840/\epsilon$ and ${Ri}_{{min}}\cong 4.14/\epsilon$, and in general if ${Ri}$ is requested to be larger than a lower threshold ${Ri}_0$ (say, ${Ri}_0=100$) this is equivalent to requiring that the drop ascent develops for a time larger than a given time, i.e. $t>t_0=[(1-\mathcal {R})R\,{Ar}\,{Ri}_0]/(U_a-U_d)$. In these experiments one has $t_{0{min}}=0.3\ \textrm {s}$ and $t_{0{max}}=24\ \textrm {s}$, corresponding to a value of $\mathcal {L} < 17\ \textrm {cm}$; only in experiment 4 does the value of $\mathcal {L}$ exceed the length of the tube, with the long drop characterized by ${Ri}<{Ri}_0$.

Figure 22 compares the ascent speed for the lock-exchange configuration of the present experiments with our theoretical results for the HB model. The curves refer to the theoretical asymptotic radius $\delta _{3\infty }$ and to the maximum flow rate. The agreement between the experimental data and the theoretical model is generally good, with most of the values thickening around the curves for $\delta _{3\infty }$, except for the experiments with a yield-stress ascending fluid. The approximations are better at low $\mathcal {M}$, and also adequately reproduce the curve knee for $\mathcal {M}<10$. Figure 23 compares the experimental ascent speed with theory based on the Ellis model, doing so only for experiments without yield stress. Again, the experimental agreement is fairly good and better than the case of the same experiments with a fluid described by the OdW model, especially for low values of $\mathcal {M}'$.

Figure 22. Ascent speed $U_a$ as a function of $\mathcal {M}$ for the experiments listed in table 1. Curves refer to the ascent speed corresponding to the asymptotic core radius $\delta _{3\infty }$, symbols are the experiments, with empty circles representing the value if the internal ascending fluid is air; dashed thick curves refer to a HB fluid with large Bingham number (not in the experiments). The descending external fluid is Newtonian. Error bars (almost the same size as the symbols) indicate two standard deviations.

Figure 23. Ascent speed $U_a$ as a function of $\mathcal {M}'$ for the experiments listed in table 2. Curves refer to the ascent speed corresponding to the asymptotic core radius $\delta _{3\infty }$, symbols are the experiments, with empty circles representing the value if the internal ascending fluid is air. The descending external fluid is Newtonian. Error bars (almost the same size as the symbols) indicate two standard deviations.

Figure 24 compares the experimental radius of the ascending internal current with the theoretical value $\delta _{3\infty }$ for the dissipative and non-dissipative cases: also for this quantity, the agreement is good, with a maximum deviation $\le 10\,\%$.

Figure 24. Drop radius in the stable node state $\delta _{3\infty }$ as a function of $\mathcal {M}$. The horizontal curves refer to the non-dissipative regime with $\mathcal {R}=1$ and to the dissipative regime with $\mathcal {R}=0.8$, symbols refer to the experiments, with empty circles for the experiments with air as internal ascending fluid. Error bars indicate two standard deviations.

In some experiments a bubble was present from the beginning on the wavefront of the ascending current, and was sometimes able to distort the geometry of the front to an appreciable extent (see, for example, the video for experiment 12 in the supplementary material available at https://doi.org/10.1017/jfm.2022.676). However, the bubble had no effect on the speed (different from the case where the bubble was injected in a subsequent phase because, in this second configuration, the bubble was able to perturb the current, see the analysis of bistability in § 7).

These data confirm, even for non-Newtonian fluids such as those used in the present experiments, that the stationarity condition and the corresponding radius of the rising inner current do not coincide with the maximum value of the flow rate; instead, they approximate the value $\delta _{3\infty }$. The discrepancies between theory and experiments can be attributed both to measurement uncertainties and to effects unaccounted for by the theoretical model, first and foremost the miscibility of the two currents, which tends to modify the velocity profiles and all related variables. It should be noted that the progressive growth of a coaxial boundary layer of fluid, having characteristics intermediate between the inner and outer fluids, involves first of all a reduction of the buoyancy; it then entails a more or less significant variation of the rheological parameters. Further, the boundary layer is presumably non-homogeneous in the direction of motion, with a minimum thickness near the front of the inner current and a progressive growth towards the tail.

7. Observed bistability

It is noteworthy that core-annular flow in vertical tubes is inherently bistable, as documented by direct numerical simulations (see Suckale et al. Reference Suckale, Qin, Picchi, Keller and Battiato2018): there exist two distinct values of the radius of the internal current that are both stable (at least to perturbations of infinitesimal amplitude). The bistability leads to two distinct regimes, one with thin-core current and the other with thick-core current, straddling the maximum of the flow rate curve.

In some experiments, transitions from one flow regime to another were observed, triggered by disturbances such as small air bubbles leading up to the inner current. Figure 25 shows the experimental front position for two Newtonian experiments in similar conditions. In both cases, the initial thick current suddenly reduces its radius with a corresponding increase in speed. Figure 25(b) also shows a return to the initial thick condition after a segment in the thin configuration. In both cases, the thick radius corresponds to $\delta _{3\infty }=0.632$ and the thin radius is still in the thick branch of the transport parameter curve, see figure 26 depicting the transport parameter and the dissipation function for the two fluids. The bistability observed in these experiments appears to be different from that documented, for example, by Suckale et al. (Reference Suckale, Qin, Picchi, Keller and Battiato2018), since the transition occurs between two configurations, both of which have a core radius greater than the radius at which the maximum flux corresponds. It appears evident that, in addition to the effects of the boundary conditions invoked by Suckale et al. (Reference Suckale, Qin, Picchi, Keller and Battiato2018), also the presence of perturbations of finite amplitude of various nature can favour the transition between two apparently stable regimes. In particular, a metastability is evident in the experiment illustrated in panel (${b}$) of figures 25 and 26, in which the disappearance of the air bubble at the top leads to a return to the starting configuration. This means that the predictability of the flow regime is even more impaired: not only is the principle of maximum dissipation invalid, which would lead to a current of maximum flow, but also bistability does not seem to be preferentially associated with the boundary conditions. Further direct or numerical experiments, accompanied perhaps by a more refined theoretical approach, are needed to discriminate between the different contributions to the onset of bistability.

Figure 25. Front position for two experiments with Newtonian ascending and descending fluids, where bistability occurs. Blue dots are the experimental results, the solid curves are the linear interpolation.

Figure 26. Transport parameter $Te$ and dissipation functions for two Newtonian experiments where bistability occurs, see figure 25, (a) with transition from A to B, and (b) with transition from A to B and vice versa. The bold curve is the transport parameter (coincident with the flow rate and with the total dissipation), the dashed and dot-dashed curves show the dissipation in the descending and ascending fluid, respectively. The arrows indicate the direction of the transition.