1. Introduction
This paper considers the entrainment of a viscoplastic fluid into a Newtonian fluid, due to the passage of bubbles across the stratified interface between the two layers. The lower layer represents an industrial suspension within which bubbles are continually created via methanogenesis, held in place by the fluid yield stress. On growing to a critical size the bubbles are released, rise and further expand. Aside from the gas release, the entrainment of suspended particles affects the turbidity of the upper layer. Turbidity, defined as the cloudiness or opaqueness of water, serves as a critical indicator of the overall health and dynamics of aquatic ecosystems. It provides insights into the presence and distribution of suspended particles, the penetration of sunlight and the consequent absorption of solar energy and the biological productivity of the aquatic ecosystem. In the industrial setting of a tailings pond, excessive turbidity, entraining pollutants and sediments from below delay land and water reclamation management.
As an example, Base Mine Lake (BML) is a constructed water body in Alberta, Canada. The lake consists of water with an average depth of 9 m, covering a layer of fluid fine tailings (FFT) of up to 45 m deep. FFT includes elements such as dissolved organic matter, suspended organic and inorganic particles that can collectively contribute to turbidity. Examination of suspended solids concentrations in the water cap above the FFT suggested oscillating bottom currents beneath wind-waves as a potential contributor to FFT erosion (Lawrence, Ward & MacKinnon Reference Lawrence, Ward and MacKinnon1991). However, the water cap's substantial depth in BML renders these currents, which attenuate exponentially with increasing depth, too weak to erode the FFT. Echo soundings conducted beneath the ice at BML have revealed the ascent of gas bubbles through the water column, as depicted in figure 4 of Lawrence, Tedford & Pieters (Reference Lawrence, Tedford and Pieters2016). These observations, combined with water turbidity data measured during winter, as reported in Tedford et al. (Reference Tedford, Halferdahl, Pieters and Lawrence2019), give rise to the hypothesis that these bubbles are likely facilitating the transportation of some solids from the FFT into the water cap, while also potentially inducing circulation and mixing within the water cap. The primary objective of the present study is to examine the entrainment of fluid from the lower layer to the upper layer, by rising bubbles.
BML is one of many gas-emitting lakes in North Eastern Alberta, Canada, resulting from the oil sands industry. As reviewed by Small et al. (Reference Small, Cho, Hashisho and Ulrich2015), these lakes cover some 175 km
$^2$ and contain a wide range of residual solids, bitumen and diluents. The lower layer in which bubbles are produced (FFT layer) has been characterised as viscoplastic (Derakhshandeh Reference Derakhshandeh2016). In such fluids the material flows only if the imposed stress surpasses the yield stress (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014). In the context of bubbles rising in a viscoplastic fluid, one might expect that stresses arise from both interfacial tension and buoyancy effects, and are resisted by the yield stress of the material. Finding the shape and position of the yield surfaces, i.e. the boundaries between the yielded and unyielded regions, is non-trivial. This complicates understanding bubble migration through viscoplastic fluids, and further the entrainment process when one layer features yield stress.
During the last few decades, many theoretical and computational studies have been performed focused at determining the yield surfaces around moving objects in viscoplastic fluids (Tsamopoulos et al. Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008; Sikorski, Tabuteau & de Bruyn Reference Sikorski, Tabuteau and de Bruyn2009; Dimakopoulos, Pavlidis & Tsamopoulos Reference Dimakopoulos, Pavlidis and Tsamopoulos2013; Lopez, Naccache & de Souza Mendes Reference Lopez, Naccache and de Souza Mendes2018; Chaparian & Frigaard Reference Chaparian and Frigaard2021; Pourzahedi et al. Reference Pourzahedi, Chaparian, Roustaei and Frigaard2022; Daneshi & Frigaard Reference Daneshi and Frigaard2023; Kordalis et al. Reference Kordalis, Pema, Androulakis, Dimakopoulos and Tsamopoulos2023). Theoretical and computational studies amongst these are often based upon a model of ideal (or simple) viscoplastic fluids, as defined in Frigaard (Reference Frigaard2019). Experimentally, deviation from ideal viscoplastic behaviour is always present. Experimental observations find fore-and-aft asymmetries in moving bubbles, even at low Reynolds numbers, characterised by an inverted teardrop shape of the bubble as it rises and a negative wake at the rear of the bubble (Sikorski et al. Reference Sikorski, Tabuteau and de Bruyn2009; Mougin, Magnin & Piau Reference Mougin, Magnin and Piau2012; Lopez et al. Reference Lopez, Naccache and de Souza Mendes2018; Zare & Frigaard Reference Zare and Frigaard2018; Pourzahedi, Zare & Frigaard Reference Pourzahedi, Zare and Frigaard2021; Daneshi & Frigaard Reference Daneshi and Frigaard2023). These features are attributed to viscoelastic behaviour of the fluid around the bubble (Moschopoulos et al. Reference Moschopoulos, Spyridakis, Varchanis, Dimakopoulos and Tsamopoulos2021). These same viscoelastic stresses result in successive released bubbles following the same pathway (Lopez et al. Reference Lopez, Naccache and de Souza Mendes2018). These damaged or weakened pathways appear relevant to bubble migration in large reservoirs. Zare, Daneshi & Frigaard (Reference Zare, Daneshi and Frigaard2021) showed that the shape and trajectory of bubbles are significantly influenced when they move within a distance (not greater than their yielded zone) from a weaker region. This influence arises from alterations in residual stress distributions.
In the context of entrainment, we are, in particular, interested in the flow around bubbles near viscoplastic–Newtonian interfaces. We conducted a series of experiments to investigate the behaviour of bubbles rising in a two-layer system consisting of viscoplastic fluid below water, as detailed in Zhao et al. (Reference Zhao, Tedford, Zare, Frigaard and Lawrence2022). Our findings revealed that as the bubble enters the water layer, it entrains viscoplastic fluids, which accumulate above the interface in a conical shape. In addition, we observed the formation of water conduits extending downward into the lower layer, when multiple bubbles passed through. As the bubbles ascended through these conduits, they displaced water, and once they exited, water returned to fill the conduits. This phenomenon holds significant environmental significance, as it can facilitate the movement of contaminants to and from submerged sediments. To the best of the authors’ knowledge, there is no other research that explores the ascent of bubbles through a rheologically non-uniform viscoplastic medium.
The concept that a moving particle permanently displaces the surrounding fluid was quantified by Darwin (Reference Darwin1953). Early investigations have been based on this classical concept, known as Darwin's drift. In these studies, two fluid layers have the same physical properties, coloured differently, and the entrainment of fluid by body translation is predicted in Stokes regime (Eames & Duursma Reference Eames and Duursma1997; Bush & Eames Reference Bush and Eames1998). Numerous experimental and numerical investigations have been conducted to explore the process of bubbles passing through a Newtonian liquid–liquid interface since then (Mercier et al. Reference Mercier, Da Cunha, Teixeira and Scofield1974; Greene, Chen & Conlin Reference Greene, Chen and Conlin1988, Reference Greene, Chen and Conlin1991; Reiter & Schwerdtfeger Reference Reiter and Schwerdtfeger1992b; Manga & Stone Reference Manga and Stone1995; Kemiha et al. Reference Kemiha, Olmos, Fei, Poncin and Li2007; Dietrich et al. Reference Dietrich, Poncin, Pheulpin and Li2008; Uemura, Ueda & Iguchi Reference Uemura, Ueda and Iguchi2010; Bonhomme et al. Reference Bonhomme, Magnaudet, Duval and Piar2012; Natsui et al. Reference Natsui, Takai, Kumagai, Kikuchi and Suzuki2014; Emery, Raghupathi & Kandlikar Reference Emery, Raghupathi and Kandlikar2018). Greene et al. (Reference Greene, Chen and Conlin1988, Reference Greene, Chen and Conlin1991) established a crude criterion for bubble penetration of the liquid–liquid interface by comparing buoyancy and interfacial tension between the liquid layers. Their work identified the conditions under which a bubble can be trapped at the interface. They noted that the entrainment volume can be eliminated by increasing the density ratio between the two liquids. Interestingly, they observed that interfacial tension had a relatively modest impact on the entrainment volume but played a role in initiating entrainment. Furthermore, they found that a decrease in the viscosity of the lower liquid resulted in a significant increase in the entrainment volume.
Reiter & Schwerdtfeger (Reference Reiter and Schwerdtfeger1992a,Reference Reiter and Schwerdtfegerb) conducted a study in which they recorded the duration of the bubble's presence at the interface (residence time), the vertical extent of the column formed beneath the bubble and the properties of drops formed in the upper phase. Shopov & Minev (Reference Shopov and Minev1992) conducted transient numerical simulations at low and intermediate Reynolds number, to study bubbles crossing the interface between immiscible fluids. Specifically focusing on bubble and interface deformations as well as film drainage dynamics, without considering the tail pinch-off phase. Their findings revealed that at very low Weber and Reynolds numbers, bubbles could take on a prolate shape (elongated in the direction of motion) when passing through the interface, particularly in cases that the upper layer is less viscous, a result supported by Manga & Stone (Reference Manga and Stone1995). In contrast, at higher Weber and Reynolds numbers, inertial forces and interaction with the interface cause bubbles to adopt oblate shapes, often forming a concavity at the rear and spherical cap shapes during the crossing.
Dietrich et al. (Reference Dietrich, Poncin, Pheulpin and Li2008) experimentally investigated interface crossing in inertial regimes, establishing a relationship between the crossing time of the interface by a bubble and the ratio of terminal velocities in the two continuous phases. They observed the coexistence of fluid motion entrained in the bubble wake and an opposing flow driven by gravity once the tail formed. Bonhomme et al. (Reference Bonhomme, Magnaudet, Duval and Piar2012) considered bubbles crossing the interface between immiscible fluids at inertial conditions, examining a wide range of bubble shape configurations from spherical to toroidal. They provided a comprehensive map of bubble shapes and entrained column geometries based on Archimedes and Bond numbers. Their findings indicated that smaller bubbles could be slowed down or even stopped at the interface, whereas larger spherical cap bubbles, with larger cross-sectional areas, generally crossed the interface more easily. Emery et al. (Reference Emery, Raghupathi and Kandlikar2018) explored the crossing of single bubbles considering the tail pinch-off mode and found that the length of the entrained liquid column was longer when the bubble's velocity remained relatively constant during the interface crossing. The tail often remained connected for an extended period before rupturing, with the liquid shell covering the bubble breaking before the column. They also studied the crossing of a stream of bubbles, providing insights into various flow regimes and the potential formation of clusters.
The majority of previous work on bubble passage through liquid–liquid interface reviewed has focused solely on Newtonian fluids. However, the passage of a single bubble through a viscoplastic–Newtonian fluid has not been studied previously to the best of the authors’ knowledge and very little is known about bubble's motion in rheologically non-uniform fluids. This paper embarks on a comprehensive exploration of the relationship between the elementary physical mechanisms at play in the entrainment of viscoplastic fluid by the rising of a single bubble, using computational methods. The study does not focus specifically on any one aspect but instead intends to provide a broad characterisation and overview of the possible flow regimes. The following sections describe the computational methods, provide technical details and discuss the results in depth, offering valuable insights into the dynamics of bubbles at fluid–fluid interfaces and subsequent entrainment.
2. Problem formulation
The present study is centred on the flow around a single gas bubble crossing an interface between two miscible liquid layers. Our primary focus revolves around entrainment, which refers to an amount of liquid from the lower layer that is transported to the upper layer by the rising bubble. Figure 1 illustrates the flow geometry examined herein.
Figure 1. Schematic of the flow geometry and coordinate system.
Our simulations are computed in an axisymmetric domain (
$[0,100 \hat {R}_b] \times [-15\hat {R}_b, 50\hat {R}_b])$ in 
$(\hat r,\hat z)$ coordinates. Here 
$\hat {R}_b$ denotes the radius of the initially spherical bubble. The liquids are initially static separated with a horizontal interface at 
$\hat z=0$. The bubble is initially placed in the lower part of the flow domain at 
$\hat z=-5$.
The flow surrounding a rising bubble is governed by the usual momentum and mass balances, and we make them dimensionless through the following approach. We scale all lengths with the equivalent radius, 
$\hat {R}_b$, of a spherical bubble of the same volume. We scale velocities with 
$\hat U$, obtained by balancing buoyancy and viscous forces, i.e. 
$\hat U=\hat {\rho }_1\hat {g}\hat {R}_b^2/\hat {\mu }_1$, where 
$\hat {g}$ is the gravitational acceleration and 
$\hat {\mu }_1$ is the plastic viscosity of the surrounding viscoplastic fluid. Pressure and other stresses are scaled with 
$\hat {\rho }_1\hat {g}\hat {R}_b$. Index 
$k$ varies for the lower layer fluid, upper layer fluid and gas phase, respectively, 
$k =1,2,3$.
Various dimensionless groups emerge. The Archimedes number (
$Ar=\hat {\rho }_1^2\hat {R}_b^3\hat {g}/\hat {\mu }_{1}^2$), represents the ratio of gravitational stresses to the effects of plastic viscous stresses. The density ratio (
$\rho _k$), represents the ratio of the density of fluids to the density of the lower layer, note that air density is negligible 
$\rho _3 \ll 1$, and 
$\rho _1 = 1$ so 
$\rho \equiv \rho _2$. The yield number is 
$Y=\hat {\tau }_{Y}\hat {R}_b/\hat {\rho }_1\hat {g}$, representing the competition between yield and buoyancy stresses. The viscosity ratio is 
$m_k=\hat {\mu }_k / \hat {\mu }_1$; note that 
$m_1=1$ and since 
$Ar$ is defined using the viscosity of the lower layer fluid, the bubble viscosity ratio
$m_3$ varies with the selected 
$Ar$. Nevertheless, for all 
$m_3$ studied we have 
$m_3 \ll 0.002$. Hence, 
$m \equiv m_2$. The Bond number represents the ratio of gravitational forces to interfacial tension effects. Here we assume the same interfacial tension between both liquids and gas, i.e. 
$\hat {\sigma }_{13} = \hat {\sigma }_{23}$, and the liquids are miscible 
$\hat {\sigma }_{12} = 0$. Following the methods of Tofighi & Yildiz (Reference Tofighi and Yildiz2013), the Bond number is 
$Bo=\hat {\rho }_1\hat {g}\hat {R}_b^2/\hat {\sigma }_{13}$, between any liquid and gas. This term vanishes between the two liquids.
The scaled momentum and continuity equations are given by
$$\begin{gather} \rho Ar [\boldsymbol{u}_t+\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u}]={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol \tau + \frac{k \boldsymbol n \delta_s}{Bo} + \rho \boldsymbol{e_{g}}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol u=0, \end{gather}$$
where 
$\boldsymbol u,p,\boldsymbol \tau$ denote the velocity, pressure and deviatoric stress, respectively.
A volumetric force formulation has been adopted for the interfacial tension, and this term appears in the right-hand side of the Navier–Stokes equation (2.1), where 
$k$ is the mean curvature, 
$\boldsymbol n$ the normal unit vector and 
$\delta _s$ is a surface Dirac 
$\delta$-function that is non-zero only on the interface. Dimensionally this term is 
$\hat {\boldsymbol {f}}= \hat {\sigma }_s k \boldsymbol n \delta _s$, where 
$\boldsymbol n \delta _s$ is approximated from 
$\boldsymbol \nabla C$, which is the gradient of volume fraction. In the numerical scheme, this term is evaluated using height functions (Popinet Reference Popinet2009).
The upper layer fluid is treated as a Newtonian fluid, whereas the lower layer fluid is considered to be either a Newtonian or a viscoplastic fluid. Where we have Newtonian liquid in our flows, this is modelled by
\begin{equation} \boldsymbol{\tau}(\boldsymbol u)=m {\dot{\boldsymbol{\gamma}}}(\boldsymbol{u}), \end{equation}
where 
$\dot {\gamma }(\boldsymbol {u})$ is the strain rate tensor:
\begin{equation} {\dot{\boldsymbol{\gamma}}}(\boldsymbol{u})= \boldsymbol{\nabla} \boldsymbol u +(\boldsymbol{\nabla} \boldsymbol{u})^{T} . \end{equation}For the viscoplastic lower fluid layer the Bingham model has been used, which has been solved computationally using a commonly used smooth regularisation method proposed by Papanastasiou (Reference Papanastasiou1987):
\begin{equation} \boldsymbol{\tau}(\boldsymbol u)=\Big[1+Y\frac{1-{\rm e}^{{-}N |\dot{\gamma}(\boldsymbol{u})|} }{|{\dot{\gamma}}(\boldsymbol{u})|}\Big]{\dot{\boldsymbol{\gamma}}}(\boldsymbol{u}). \end{equation} Here 
$\hat {\tau }_{Y}$ denotes the yield stress of the Bingham fluid. The regularisation parameter 
$N \gg 1$ controls the closeness of approximation to the exact Bingham fluid model: 
$1/N$ represents a small (scaled) strain rate below which the fluid becomes very viscous, i.e. for 
$N |\dot {\gamma }(\boldsymbol {u})| \ll 1$, we find at leading order (2.5) becomes
\begin{equation} \boldsymbol{\tau}(\boldsymbol u) \sim NY {\dot{\boldsymbol{\gamma}}}(\boldsymbol{u}) . \end{equation}
In dimensional terms the strain rate and 
$1/N$ have been scaled with 
$\hat U /\hat {R}_b = \hat {\rho }_1\hat {g}\hat {R}_b / \hat {\mu }_1$.
Both fluids are considered miscible with density ratios of (
$\rho = 1$) and (
$\rho = 0.7$). Thus, there is no additional interfacial term required in (2.1). The experimental timescale and length-scales are such that the Péclet number between liquids is extremely large and molecular diffusive effects can be neglected, i.e. the miscibility is physically appropriate.
In steady state, the drag force on the bubble is balanced by the buoyancy. The drag coefficient 
$C_d$ compares the drag force with the inertial force on the bubble:
\begin{equation} C_d=\frac{32\hat{R}^3_b\hat{g}}{3\hat{U}^2_b\hat{W}^2}, \end{equation}
where 
$\hat {W}$ denotes the width of the steadily rising bubble.
3. Computational method
The computation has been performed using a volume-of-fluid (VoF) method. For three-phase flows we introduce the volume fraction 
$C_i$ for each phase and the density and viscosity of the mixture are related to the volume fraction of each phase through the linear laws as
\begin{equation} \rho=\sum_1^3 C_i \rho_i, \quad \mu=\sum_1^3 C_i \mu_i. \end{equation}Each phase is advected
\begin{equation} \partial_tC_i + \boldsymbol{\nabla} \boldsymbol{\cdot}(C_i \boldsymbol{u})=0. \end{equation}The advection equation (3.2) is solved using a piecewise-linear geometrical VoF scheme that is generalised for the quad/octree spatial discretisation; for further detail, please refer to Scardovelli & Zaleski (Reference Scardovelli and Zaleski1999).
The quadtree adaptive mesh refinement technique has been used to accurately track the interfaces between the two and three phases. The local mesh refinement densifies the grid in regions of strong spatial variations of the velocity, fluid concentrations and the strain rate 
$|\dot {\gamma }(\boldsymbol {u})|$. The algorithm is implemented in the open-source multiphase flow solver Basilisk (based on the algorithms of the former Gerris), which is specifically suited to such applications (Popinet Reference Popinet2003, Reference Popinet2009; Fuster & Popinet Reference Fuster and Popinet2018). The system of equations is resolved using a time-splitting projection method described in Popinet (Reference Popinet2009). A physically motivated maximum timestep is imposed and further controlled via a CFL constraint. The velocity advection term is estimated using the Bell–Colella–Glaz second-order upwind scheme at the intermediate timestep. The values of viscosity, density and pressure at the intermediate timestep are also used.
Boundary conditions for the simulations are no-slip on the walls at the side and top of the domain. The pressure as well as the 
$z$-derivative of the normal velocity is set to zero at the top of the domain (
$z=10$). Within the range of our dimensionless parameters, the size of yield surfaces developing around the bubbles is expected not to exceed five times the size of the bubble; see Zare et al. (Reference Zare, Daneshi and Frigaard2021). The width of the domain is significantly larger than the width of yielded zone and the height of the domain is tall enough to allow for transients of interest to be captured.
Basilisk employs the interpolation error between a field value at a grid point on grid level 
$n$ and its interpolated value from a coarser grid level 
$n-1$ as a criterion to decide whether to adjust the grid locally. This adjustment involves either merging four squares into a parent square or dividing a square into four subsquares in two dimensions, as appropriate. The grid is refined or coarsened depending on whether this interpolation error is higher or lower than 
$10^{-3}$, 
$10^{-3}$ and 
$10^{-9}$ for velocity, concentrations and strain rate, respectively. The domain is initially divided into 
$128 \times 128$ cells. Each of these may subdivided according to the refinement, to a maximum of 9 quadtree levels. Spurious currents are also minimised by using a height function method to estimate the interface curvature. Basilisk has been widely validated against analytical, numerical and experimental interfacial flows and has been frequently used for Newtonian bubble flow calculations (Balla et al. Reference Balla, Tripathi, Sahu, Karapetsas and Matar2019; Berny et al. Reference Berny, Deike, Séon and Popinet2020; Zhang et al. Reference Zhang, Dabiri, Chen and You2020).
3.1. Validation
We have validated the code by benchmarking our results with the experimental results obtained by Bonhomme et al. (Reference Bonhomme, Magnaudet, Duval and Piar2012). The experiments were based on a single gas bubble crossing an interface between a lower phase made of water or water plus glycerin and an upper, slightly lighter, phase made of silicon oil.
Here we adopt the procedure suggested by Tofighi & Yildiz (Reference Tofighi and Yildiz2013) for modelling a three-phase system within the VoF framework. The interfacial tension between phases 
$i,$ and 
$j$ are denoted by 
$\sigma _{ij}$. These are transformed to phase specific surface tensions using the following relations:
$$\begin{gather} \sigma_1 = 0.5(\sigma_{12}+\sigma_{13}-\sigma_{23}), \end{gather}$$
$$\begin{gather}\sigma_2 = 0.5(\sigma_{12}+\sigma_{23}-\sigma_{13}), \end{gather}$$
$$\begin{gather}\sigma_3 = 0.5(\sigma_{13}+\sigma_{23}-\sigma_{12}), \end{gather}$$
and the 
$\sigma _i$ are used within the VoF method as explained in Tofighi & Yildiz (Reference Tofighi and Yildiz2013). For this validation example, we have used the values of 
$\sigma _{ij}$ from the experiments of Bonhomme et al. (Reference Bonhomme, Magnaudet, Duval and Piar2012). Dimensionless numbers from Bonhomme et al. (Reference Bonhomme, Magnaudet, Duval and Piar2012) 
$(Ar_B, Bo_B, \varLambda, R$) are translated to those used in this study as follows: 
$(Ar=Ar_B^2/8,\ Bo_k=Bo_{k,B}/4,\ m=\varLambda, \rho =1-R)$.
The results are plotted in figure 2(b) against those obtained by Bonhomme et al. (Reference Bonhomme, Magnaudet, Duval and Piar2012) in figure 2(a). The bubble crosses the interface and tows a column of heavy fluid behind it. The shape of bubble before and after interface, and the column of heavy fluid have been found to compare well with the experimental results.
Figure 2. Evolution of shape of a bubble and interface between two fluids as the bubble crosses the interface. Parameters: 
$Ar=226.8$, 
$Bo=7.375$, 
$Y=0$, 
$\rho =0.794$ and 
$m=1.11$. Here 
$\sigma _{23}/\sigma _{12}$ is 0.425 and 
$\sigma _{13}/\sigma _{12}$ is 0.615. Figures show (a) the experimental observation (Bonhomme et al. Reference Bonhomme, Magnaudet, Duval and Piar2012) and (b) our computational results. As the bubble crosses the fluid–fluid interface, it entrains a column of heavy fluid into the upper layer.
3.2. Scope of study
Our study is based on 167 simulations covering a wide range of governing parameters. The values of governing parameters are selected from a wide range: 
$(Y, Bo, \rho, m)= ([0, 0.1, 0.15], [0.1, 1, 5, 10, 50], [0.7, 1], [0.1, 1, 10])$ and various 
$Ar=1, 5, 10, 50, 500$.
To quantify the entrainment, the volume of liquid from the lower layer (
$\,f_1$) that has been transported to the upper layer (
$\,f_2$) is calculated, see figure 3. As the bubble crosses the liquid–liquid interface, it displaces the initially horizontal interface between the liquid layers and it carries liquid from the lower layer to the upper layer either as a towed column or trapped in the recirculating wake. In this context, we have estimated the entrainment volume by
\begin{equation} V_e = \frac{2{\rm \pi} \varSigma_{z>0} \hat r_{ij} \hat{\varDelta}_{ij}^2C_{ij}}{4/3{\rm \pi}\hat{R}_b^3}, \end{equation}
where 
$C_{ij}$ is the volume fraction of 
$f_1$, 
$\hat {\varDelta }$ represents the size of cells encompassing 
$f_1$ located at 
$r_i,z_j$ and 
$z>0$, and 
$r_{ij}$ represents the radial distance of cells located at 
$r_i,z_j$ to the axis of symmetry. This value has been normalised with the volume of rising bubble.
Figure 3. As the bubble crosses the liquid–liquid interface (between 
$f_1$ and 
$f_2$), it displaces the initially horizontal interface and tows a column of 
$f_1$. The mesh has been adapted and refined at the interfaces and inside the bubble. The entrainment volume is essentially the cumulative volume of cells 
$z>0$, encompassing 
$f_1$.
In our code benchmarking, we computed 
$Ar=5$, 
$Bo=5$, 
$Y=0.1$, 
$\rho =1$, 
$m=1$ using various values of maximum level of refinement (
$MLR$), and calculated the volume of entrainment 
$V_e$. The results are presented in figure 4. Due to the automated quadtree mesh refinement, the number of cells (
$NC$) varies at each time step, hence the 
$NC$ here shows the numbers of cells when the bubble centroid reaches 
$z_b \sim 5$. Notably, the liquid–liquid interface shows good alignment for computations with 
$MLR \geq 7$. To assess convergence concerning 
$MLR$, we consider the solution at 
$MLR=9$ as the converged value and calculate the relative error of 
$V_e$. Beyond 
$MLR \geq 8$, the relative error of 
$V_e$ remains relatively stable, approximately at 
$10^{-2}$. For the subsequent computations in this paper, we adhere to 
$MLR=9$.
Figure 4. Interface between the lower and upper liquid layers at different 
$MLR$ for 
$Ar=5$, 
$Bo=5$, 
$Y=0.1$, 
$\rho =1$ and 
$m=1$: 
$MLR={6}$, 
${7}$, 
${8}$, 
${9}$ with red, green, blue and black colours, respectively. The inset graph displays the relative error in the volume of entrainment (
$V_e$) calculated at the time bubble reaches (
$z_b \simeq 5$) for varying numbers of cells (
$NC$) relative to the solution at 
$MLR=9$ and 
$NC=126289$. The 
$NC$ values corresponding to 
$MLR=6,7,8$ and 
$10$ are 6745, 11 425, 21 295 and 294 001, respectively.
4. Parametric effects
In this section we present our results of examining the impact of the governing dimensionless parameters on the entrainment process.
4.1. Effect of bubble's shape
It is of interest to examine how the change in the bubble shape influences the entrainment. We consider two scenarios, (a) both upper and lower liquid layers are Newtonian and (b) the lower layer is a viscoplastic (Bingham) fluid. We begin with a relatively low inertia regime 
$(Ar=5)$, where the liquid layers are isodense and isoviscous i.e. 
$\rho =1$, 
$m=1$, and we study the effect of changing 
$(Bo)$, figure 5.
Figure 5. Effect of variation of bubble's initial shape on the entrainment at 
$Ar=5$, 
$\rho = 1$, 
$m= 1$, for (a) 
$Y=0$ and (b) 
$Y=0.1$; from left to right 
$Bo=0.1$, 10 and 50. The bubble and the evolution of the liquid–liquid interface as bubble centroid reaches five specified heights, 
$z= -1$, 5, 15, 25, 35 with black, red, green, blue and brown colours, respectively, is shown in each panel. Due to break up of the entrained liquid column small droplets form and in cases that the bubbles breaks up at the interface small satellite bubbles form. The former is indicated with D and the latter is indicated with BD.
When the liquid layers are Newtonian, given values of 
$\rho$ and 
$m$, the bubble does indeed rise within a single liquid. The bubble reaches its steady-state shape before reaching the interface between the two liquids and maintains this shape after crossing the interface. By increasing 
$Bo$ from 0.1 to 10, the shape of the bubble undergoes a transformation from a spherical shape to a ellipsoidal cap, with wider profile. Tips also develop at the sides of the bubble and the rear side becomes indented. At 
$Bo=50$, a bubble with a smooth and steady skirt forms. Often, the formation of skirted bubbles is a prelude to the onset of instabilities for a rising bubble. Here we have not considered instabilities directly, but simply advise caution for such shapes. Even when not visibly unstable, there is no guarantee that our assumption of axisymmetry is not preventing asymmetric instabilities.
As depicted in figure 5(a), the entrainment behaviour varies accordingly. At the beginning, the spherical bubble pulls a portion of the lower layer along with it. After travelling a distance of 
$z_b \approx 15$, it eventually separates from the interface. Therefore the main contributing factor to entrainment (at 
$Bo=0.1$) is the displacement of the interface. When the bubble adopts a cap-shaped configuration, its broader profile leads to a more significant displacement of the fluid at the interface compared with when it is spherical. Consequently, a column of liquid from the lower layer accumulates behind the rising bubble. In the case of the indented cap-shaped bubble (
$Bo=10$), this liquid column extends beyond 
$z=35$. The detachment of the bubble from this towed liquid column occurs later, and no entrained liquid remains attached to the bubble. When 
$Bo$ increases to 50, apart from the interface displacement, the bubble entrains liquid from the lower layer within the primary wake, constituting approximately 20 % of the bubble's volume, forming a new compound. This compound then separates from the column and continues its ascent. At high 
$Bo$, similar to low 
$Bo$ values, no liquid column forms.
The results indicate that all bubbles rising in a Newtonian liquid cause displacement of the interface. Cap-shaped bubbles, in particular, not only induce interface displacement but also tow behind a liquid column. The thickness of the liquid column is influenced by the width of the bubbles and, eventually, the bubble detaches from the column. On the other hand, skirted bubbles rise as a separate entity, carrying along a portion of the liquid from the lower layer. As part of the detachment, for Newtonian fluids, we often see a break up of the entrained liquid column into small droplets. This phenomenon is marked with a D in figure 5(a), and subsequent figures.
We now turn our attention to the entrainment process when the lower layer is a Bingham fluid. At 
$Y = 0.1$, corresponding to 
$Bo$ value, the bubble develops different steady-state shapes before reaching the interface between two fluid layers. The steady-state bubble shapes varies from a spherical shape to a prolate shape, as 
$Bo$ increases. For small values of 
$Bo \leqslant 1$, the bubble maintains an almost spherical shape, primarily due to the significant influence of interfacial tension. However, at medium values of 
$Bo$, the bubble forms a prolate shape. In this case, gravitational forces surpass interfacial tension effects, resulting in an decreased rate of strain, subsequently increased effective viscosity, around the equatorial plane of the bubble compared to the poles. Consequently, the bubble begins to deform preferentially in the direction of its poles, taking a prolate shape. As 
$Bo$ increases and the influence of interfacial tension decreases, the prolate bubble becomes more elongated. For large 
$Bo=50$, the bubble retains an indentation. These bubble shapes are in agreement with those predicted earlier by Tsamopoulos et al. (Reference Tsamopoulos, Dimakopoulos, Chatzidai, Karapetsas and Pavlidis2008). At lower 
$Bo$ we see a detachment of the entrained liquid column and breakup into small droplets, as before for the Newtonian fluids.
Remarkably, as the bubble approaches the interface, it undergoes increased elongation, as depicted in figure 6(a). This elongation is primarily driven by an elevated rate of strain occurring near the north pole where the bubble intersects the interface. This elevated strain rate leads to a reduction in the effective viscosity near the pole, and decreased rate of strain around the equatorial plane, subsequently causing further elongation of the bubble. This elongation becomes even more pronounced with higher values of 
$Bo$ and 
$Y$, as illustrated in figure 6(b).
Figure 6. (a) Contours of strain rate on the computational cells as a bubble reaches the interface at 
$Ar=5$, 
$Bo=10$, 
$Y=0.1$, 
$\rho =1$ and 
$m=1$. (b) The profile of a bubble for 
$Ar=5$, 
$\rho =1$ and 
$m=1$ at different values of 
$Bo$, including 0.1, 5, 10 and 50, as well as for 
$Y=0$, 
$0.1$ and 
$0.15$.
We observe two different types of bubble breakup at 
$Bo=10$ and 
$Bo=50$. At a moderate value of 
$Bo$, the bubble experiences breakup as it enters the upper Newtonian fluid, whereas at a high 
$Bo$, the bubble undergoes breakup within the viscoplastic fluid layer, with a small portion of it remaining confined within this layer. Both breakup phenomena occur when viscous shear forces surpass the capillary forces, indicating the presence of a critical capillary number (
$Ca = {Bo}/{Ar}$). The initial breakup of rising bubbles in the viscoplastic fluid happens at 
$Ca_c \approx 5$ (over the set of our computational parameters). Breakup of drops has been seen at the tip when 
$m_b = \hat {\mu }_b/\hat {\mu }_p \ll 1$ and for cases with capillary number above a 
$Ca_c$ (Chu et al. Reference Chu, Finch, Bournival, Ata, Hamlett and Pugh2019). We also see a breakup/detachment process later in the flow as the bubble rises. However, unlike the small 
$Bo$ results the droplets of entrained fluid often contain smaller residual bubbles of the gas, so that the droplets encapsulate bubbles on breakup, particularly at larger 
$Bo$. We mark figures where this phenomenon occurs with BD.
The contribution of entrainment through the displacement of the interface diminishes when there is a viscoplastic fluid in the lower layer. It is mainly due to the fact that as the bubble approaches the interface, buoyancy yields only a partial region of the lower fluid while the rest remains unyielded. It can be seen as the passage of bubble from a narrow yielded column to the upper Newtonian layer. Therefore, a limited radial length of approximately 5 of the interface becomes displaced and distances from 
$z=0$. For 
$Bo =10$ and 
$50$, some Bingham fluid becomes trapped behind the rising bubble. At 
$Bo=10$, an inverted cone of viscoplastic fluid remains attached to the bubble and forms a new configuration. The rising bubble and attached cone break off from the towed liquid column around 
$z=15$, while the bubble remains connected to the corresponding Newtonian towed column. In general, the rising bubble transports less fluid to the upper layer and entrainment decreases significantly.
The variation of the volume of entrainment as the bubble rises, as well as the rise speed of the bubble for both scenarios (i.e. Newtonian and Bingham lower layers) are shown in figure 7. The entrainment volume (
$V_e$) vs the bubble's position, as shown in figure 7(a), demonstrates that in the case of Newtonian lower layers, spherical-shaped bubbles transport a smaller amount of liquid to the upper layer when compared with cap-shaped or skirted bubbles. This can be attributed to two factors: (i) the displacement at the interface increases proportionally with the area of the bubble; (ii) the towed column becomes thicker for the same reason. As a result, cap-shaped or skirted bubbles are more efficient in bringing liquid to the upper layer. In contrast, for Bingham lower layers, spherical-shaped bubbles transport a larger quantity of liquid to the upper layer compared with cap-shaped or skirted bubbles. This is because, as 
$Bo$ increases, the width of the bubble's profile decreases in the lower viscoplastic fluid, leading to a thin entrained column and a reduction in entrainment.
Figure 7. Evolution of (a) the volume of entrainment 
$V_e$, and (b) the rise speed (
$U_b$) of bubble as a function of the dimensionless bubble position for the cases considered in figure 5.
Surprisingly, we do not see a noticeable difference in the entrainment for 
$Bo=5$, 10 or 50. By comparing the shape of bubbles and their rising velocity, one might notice that despite the variety in the shapes, the width of bubbles and their rising velocity are quite similar, see figure 7(b). It implies that the effect of entrainment on drag forces imposed on the bubble is minimal. Despite that, the rising velocity of the spherical bubble at 
$Y=0.1$ increases gradually and then reaches a steady-state value similar to the Newtonian counterpart figure 7(b). However, the rising velocity at 
$Bo=10$ and 
$50$ undergoes an acceleration followed by a deceleration. The speed of rising bubble accelerates as it enters the Newtonian upper layer, but due to attached Bingham fluid it undergoes a deceleration.
4.2. Effect of inertia
We investigated the effect of inertia on the transport of both Newtonian and viscoplastic fluids at moderate and high values of Bond numbers, 
$Bo=5$ and 
$50$.
The Newtonian results for 
$Bo=5$ are presented in figure 8(a). As we discussed in § 4.1, as the bubble approaches the initially horizontal Newtonian–Newtonian interface, the interface becomes displaced. This displacement occurs within a specific radial domain and dissipates at a critical radius denoted 
$\lambda$, see figure 5. By increasing 
$Ar$, 
$\lambda$ decreases 
$\sim$10, 9.7, 8 and 4. Here 
$Ar$ can be interpreted as a Reynolds number, with characteristic velocity 
$\hat U$ from the balance of buoyant and viscous stresses. As 
$Ar$ (i.e. 
$\hat U$) increases the ratio of advective to viscous timescales decreases and disturbance of the interface behind the bubble (i.e. 
$\lambda$) decreases.
Figure 8. Effect of increasing 
$Ar$ on the entrainment at moderate 
$Bo$: 
$Bo=5$ and 
$\rho =1$, 
$m=1$ for (a) 
$Y=0$ and (b) 
$Y=0.1$; from left to right 
$Ar=5,\ 10,\ 50$ and 
$500$. For colours refer to the caption of figure 5.
The bubble develops different shapes before reaching the interface. It undergoes a transition from a spherical cap to an ellipsoidal cap with a wider profile as 
$Ar$ increases. This alteration in shape leads to a wider towed column. Subsequently, the length of the towed column extends as 
$Ar$ increases. For 
$Ar=5$ and 
$10$, no liquid remains attached to the bubble and the bubble leaves off the towed column before reaching top of the computational domain (
$z= 50$). At 
$Ar=50$, some liquid remains initially attached to the bubble; however, it gradually lags behind and ultimately separates from the bubble. In the highly inertial regime, specifically at 
$Ar=500$, the entrained liquid divides into two parts: a towed column and a recirculating torus behind the bubble. The towed column becomes progressively thinner, and the recirculating bodies gradually move further away from the bubble.
The results for 
$Y=0.1$ are depicted in figure 8(b). In this scenario, there is no significant change in 
$\lambda$ with the increase in 
$Ar$, and only a small portion of the interface moves away from 
$z=0$. The interface remains nearly flat as the bubble crosses it and ascends further away. The shape of bubble before reaching the interface varies by 
$Ar$. As 
$Ar$ increases, the aspect ratio decreases, causing the profile to shift from a prolate shape to one more akin to an oblate shape with a flat rear side. The main contribution to the entrainment comes from the Bingham fluid attached to the rising bubble at 
$Ar=50$ and 
$500$. In contrast with the Newtonian counterpart, the inverted cone of viscoplastic fluid remains attached to the bubble and does not decay with time/distance. The fluid in the inverted cone is partly unyielded. At 
$Ar=500$, the rate of strain within the Newtonian fluid surrounding the towed column is higher than that of Bingham fluid inside the column. Consequently, the interface becomes concave, leading to a significantly different flow pattern around the bubble and entrainment, compared to the Newtonian counterpart.
The results for 
$Bo=50$ and 
$Ar=1,\ 5,\ 10$ and 
$50$ are shown in figure 9. When the lower layer is a Newtonian fluid, by increasing 
$Ar$, the profile of bubble becomes wider and unstable skirts form. As a result, more liquid becomes trapped behind the bubble and larger ascending compounds form. At 
$Ar=50$ and 
$Bo=50$, the bubble becomes highly unstable and splits into two parts at 
$z \sim 25$.
Figure 9. Effect of increasing 
$Ar$ on the entrainment at large 
$Bo$: 
$Bo=50$ and 
$\rho =1$, 
$m=1$, for (a) 
$Y=0$ and (b) 
$Y=0.1$; from left to right 
$Ar=1,\ 5,\ 10$ and 
$50$. For colours refer to the caption of figure 5.
The impact of inertia on the displacement of a viscoplastic fluid (
$Y=0.1$) is depicted in figure 9(b). As 
$Ar$ increases, the bubble becomes wider and consequently, the thickness of the entrained column increases. At low and medium values of 
$Ar$, where 
$Ca \geqslant 5$, bubble breaks up before entering the upper layer and some air becomes entrapped within the viscoplastic fluid. The increase in 
$Ar$, causes a greater amount of liquid to trap behind the ascending bubble within the upper Newtonian layer.
4.3. Upper and lower layers with similar properties
So far our results show how a rising bubble induces flow and transports liquid at specific values of 
$Ar$ and 
$Bo$. As discussed in §§ 4.1 and 4.2, the lower layer liquid could be transported by either the displacement of the interface (in the case of spherical shaped bubbles), or the towed column (in the case of cap-shaped bubbles) or the wake behind a rising (skirted-shaped) bubble. When buoyancy forces are dominant, high-
$Bo$ regime, the liquid attached to the cap-shaped bubbles break off from the column and rises up with the bubble. Increasing inertia in each regime amplifies the process and increases the retention time at low and moderate 
$Bo$ values.
4.4. Effect of viscosity
When liquid layers have different viscosities but similar densities, the viscosity ratio (
$m$) plays a significant role in altering the Archimedes number in the upper layer. An increase in 
$m$ results in a decrease of 
$Ar$ by 
${1}/{m^2}$ while the Bond number remains unaffected. We investigated the effect of varying 
$m$ on the rise of bubble and entrainment in moderate 
$Ar$ and 
$Bo$ regime, where bubbles form an indented cap shape. The passage of the bubbles and subsequent entrainment at 
$m=0.1,\ 1$ and 
$10$ are shown in figure 10.
Figure 10. Effect of variation of the viscosity ratio (
$m$) on the entrainment for 
$Ar= 10$, 
$Bo=5$, 
$\rho =1$ and (a) 
$Y=0$ and (b) 
$Y=0.1$; 
$m$ increases from left to right as 0.1, 1 and 10. For colours refer to the caption of figure 5.
The results for the case in which the upper and lower layers are Newtonian are shown in figure 10(a). By increasing 
$m$, 
$\lambda$ increases 
$\sim$8, 14 and 46, which is counterintuitive. By increasing 
$m$, the bubble experiences deceleration and consequently the advective to viscous time scale increases and the disturbance of the interface behind the bubble (i.e. 
$\lambda$) increases.
When the bubble enters the upper layer at 
$m=0.1$, the Archimedes in the upper layer increases to 1000, so we expect that the bubble becomes wider and the entrained column becomes thicker. However, in contrast to these expectations, the increase in Archimedes (to 1000) does not lead to a proportional increase in the entrainment, and the width of the liquid column does not become significantly thicker with respect to 
$Ar=10$ and 
$m=1$. This could be attributed to the high shear stresses that are generated in the upper layer, which cause the towed liquid column to have a thickness comparable to the width of the bubble. On the other hand, when the viscosity ratio is 
$m=10$, as the bubble enters the upper layer, the Archimedes number decreases to 0.1. In this situation, the bubble's inertia becomes less significant than the viscous stresses. Consequently, the bubble deforms and takes on a spherical shape, as expected. The entrainment in this case is caused by a significant displacement of the interface. Due to the initial shape of the bubble and lower viscosity of the lower layer, which results in generation of higher strain rates, the bubble displaces a larger area of the interface and creates an initially thick towed column. However, as the flow regime transitions and the shape of the bubble changes, it separates from the Newtonian–Newtonian interface and the liquid column is not extended to 
$z>15$.
When the lower layer is a Bingham fluid, as illustrated in figure 10(b), as the bubble approaches the interface, its aspect ratio increases. This elongation becomes more pronounced when the Newtonian viscosity is less than the plastic viscosity (
$m=0.1$). In this case, the bubble breaks up when it has nearly completely emerged from the interface. Formation of break-up bubbles is caused by the acceleration upon entering the upper layer as it has a higher Archimedes number. These break-up bubbles become trapped between the two liquid layers. In the case of (
$m=10$), the bubble assumes a more spherical shape upon entering the upper layer. The bubble initially tows a liquid column, and therefore the trailing part of the bubble becomes slightly unstable, giving rise to the formation of small drops. The entrainment of viscoplastic fluids manifests as the towed liquid column for 
$m\leqslant 1$. When 
$m >1$, the towed column is short and truncated around 
$z=5$. In contrast to Newtonian fluids, there is no significant change in 
$\lambda$ with 
$m$ .
4.5. Effect of density variation
By keeping liquid layers isoviscous and reducing the density of the upper layer to 70 % of the lower layer, the Archimedes decreases by a factor of 0.49 and the Bond number decreases by a factor of 0.7. To study these effects, we considered 
$Bo=50$, at 
$Ar=10$, for 
$Y=0$ and 
$0.1$ the results are shown in figure 11.
Figure 11. Effect of density variation on the entrainment for 
$Ar= 10$, 
$Bo=50$ and 
$m=1$: (a) 
$Y=0$; (b) 
$Y=0.1$. The left figures show the entrainment for 
$\rho =1$ and the right figures show the entrainment when 
$\rho =0.7$. For colours refer to the caption of figure 5. (c) Fragmentation of bubble into smaller bubbles as it enters the upper layer. (d) The volume of entrainment (
$V_e$) as a function of bubble position. (e) Entrapment of heavy Bingham fluid as well as a small bubble in the wake of rising bubble.
The reduction in Archimedes and Bond numbers, both suggest the bubble deformation to a more indented spherical cap shape upon entering the upper layer. When the lower layer is Newtonian, this transformation causes the bubble to fragment into smaller bubbles as it departs from the interface, see figure 11(c). Interestingly, the small bubbles generated above the interface are encapsulated within the towed heavy Newtonian liquid and transport some of the surrounding liquid and prevent the towed column from descending.
The results for the case where both upper and lower layers are Newtonian are displayed in figure 11(a). As the bubble enters the lighter fluid, the displacement at the interface and the subsequent entrainment noticeably decreases. The column of the entrained liquid thins and nearly diminishes as the bubble reaches 
$z=5$. In contrast to the isodense liquid layers in which the interface rises until the bubble detaches, here the interface gradually settles back.
The results for 
$Y=0.1$ are shown in figure 11(b). When the lower layer is a viscoplastic fluid, the bubble tows a column of this fluid, but the towed liquid does not change with time; it remains stationary above the interface. The density difference is not sufficient to cause the entrained liquid to settle back after the bubble breaks off. However, this density difference becomes apparent in the thinned towed column. Interestingly, there is some heavy Bingham fluid trapped in the wake of the rising bubble as well as a small bubble (figure 11e). The small bubble has been formed during shape transformation above the liquid–liquid interface,
In both cases, i.e. Newtonian and Bingham lower layer, bubbles tow some liquid as a column, which settles down shortly after the bubbles travels a distance of 
$z \sim 5$, and also entrap some in their wake. Initially, a greater amount of liquid is transported from a Newtonian layer compared with a viscoplastic fluid layer. However, a portion of the Newtonian entrained liquid settles down due to the density difference. Consequently, the entrained heavy Newtonian liquid is less than the entrained heavy viscoplastic liquid; see figure 11(d).
The results of high-inertia regime with 
$Ar=500$ are shown in figure 12. When the bubble rises from a heavy Newtonian layer to a light Newtonian layer, it assumes a wide shape both below and above the interface. As a result, the bubble does not undergo fragmentation as it enters the lighter fluid, and the heavier towed column splashes onto the liquid–liquid interface; see figure 12(a). Interestingly, the interface connecting the bubble to the towed column changes from convex to concave. The negative buoyancy force causes the strain rates in the towed column near the equatorial plane around the bubble to become lower than those in the surrounding fluid. This results in a change in the concavity of the interface and absence of the recirculating regions behind the rising bubble. The entrained heavy liquid, visible as the towed column and wake behind the bubble, gradually settle down.
Figure 12. Effect of density variation on the entrainment for 
$Ar= 500$, 
$Bo=5$ and 
$m=1$: (a) 
$Y=0$; (b) 
$Y=0.1$. The left figures show the entrainment for 
$\rho =1$ and the right figures show the entrainment when 
$\rho =0.7$. For colours refer to the caption of figure 5.
The results for the case with a Bingham fluid in the lower layer are shown in figure 12(b). The main difference between the entrainment at 
$\rho =0.7$ and at 
$\rho =1$ is in the entrapment of liquid behind the rising bubble in the upper layer. At 
$\rho =1$, an inverted cone of Bingham fluid forms behind the rising bubble, while at 
$\rho =0.7$, the liquid behind separates from the rising bubble and falls down.
4.6. Effect of density and viscosity variation
It is of interest to explore the entrainment in the case that Archimedes number increases as the bubble crosses the interface, but the Bond number decreases. The results obtained for 
$Ar=10$, 
$Bo=50$, 
$\rho =0.7$ at 
$m=1$ and 
$m=0.1$ are shown in figure 13.
Figure 13. Effect of density and viscosity variation on the entrainment for 
$Ar= 10$, 
$Bo=50$ and 
$\rho =0.7$: (a) 
$Y=0$; (b) 
$Y=0.1$. The left figures show the entrainment for 
$m=1$ and the right figures show the entrainment when 
$m=0.1$. For colours refer to the caption of figure 5.
In both cases of the Newtonian and Bingham lower layers, as the bubble rises and encounters the upper layer, it adopts a wider profile and develops an unstable skirted shape. The bubble's profile accommodates transportation of more liquid to the upper layer. However, notably, the entrained liquid does not remain within the wake of the skirted bubble and instead descends. This presents an interesting contrast with the case of 
$\rho =0.7$ and 
$m=1$. Furthermore, the bubble breaks up in the lower Bingham fluid before reaching the interface, 
$Ca=5$, and a small portion of bubble remains within the lower layer.
5. Regimes
As the bubble crosses the interface, the interface changes its shape and a portion of lower layer fluid is entrained into the upper layer fluid. During the ascent of the bubble in the upper layer, the entrainment volume changes by time. We have identified the following four different characteristic behaviours. (i) Plateau regime. The rising bubble breaks off from the liquid–liquid interface, either with or without carrying a portion of entrained fluid in its trailing wake. Upon the separation of the bubble from the entrained fluid, it does not return to the lower layer. Consequently the entrained volume reaches a plateau and remains relatively unchanged over time. (ii) Descending–steady-state regime. Initially, the bubble transports some fluid, but the fluid separates from the liquid–liquid interface. It may do so either with or without carrying some fluid in its trailing wake. In either case, a portion of the entrained fluid settles down. Eventually, the entrained volume reaches a plateau and does not change over time. (iii) Descending regime. If the reduction in 
$V_e$, explained in (ii), continues, it is called descending. (iv) Ascending regime. The rising bubble remains connected to the liquid–liquid interface as it rises up, so the entrained volume is ascending, i.e. 
$V_e$ keeps increasing.
The variation of entrainment with bubble position for these four distinct entrainment regimes are shown in figure 14(a) and representative examples of the evolution of the interfaces are shown in figure 14(b), for each case.
Figure 14. (a) The variation of entrainment with bubble position for four distinct entrainment regimes. (b) Representative examples of four entrainment regimes. For colours refer to the caption of figure 5.
It is of interest to identify key parameters controlling the entrainment regimes. The entrainment regimes for 
$Y=0$ and 
$0.1$ and 
$\rho =1$ and 
$0.7$ are investigated. The results are presented within the 
$Ar$–
$Bo$ domain, with the size of the symbol indicating the value of 
$m$.
The entrainment regimes for 
$\rho =1$ and 
$Y=0$ are shown in figure 15(a). For 
$\rho =1$, there is no driving force for entrained fluid to settle. In this case, plateau and ascending regimes occur. At low 
$Bo$, the dominant regime is the plateau regime. The bubble assumes a spherical shape, and the entrainment happens only in the form of displacement of the interface. At moderate 
$Bo$ and 
$m=0.1$ and 
$1$, the gravitational forces are stronger than interfacial tension, and the bubble remains connected to the liquid–liquid interface and a towed column forms. There are no significant forces leading to break off of bubble with the column. Thus, we may observe the ascending regime in this case. This is the only case where the ascending regime occurs.
Figure 15. Classification of entrainment regimes over the 
$Ar$–
$Bo$ domain for (a) 
$\rho =1$ and 
$Y=0$ and (b) 
$\rho =1$ and 
$Y=0.1$. The size of symbols represents the value of 
$m$ at each point, i.e. 
$m={0.1}$, 1, 10 with magenta, black and cyan colours, respectively. Regimes are represented by symbols as follows: 
$\square$ plateau; 
$\diamond$ descending–steady-state; 
$\triangledown$ descending; and 
$\triangle$ ascending.
The entrainment regimes for 
$\rho =1$ and 
$Y=0.1$ are shown in figure 15(b). In this case plateau and descending–steady-state regimes occur. At low 
$Bo$, the dominant regime is the plateau regime. The bubble assumes a spherical shape, and the entrainment is limited to formation of low towed columns above the interface. These transported viscoplastic fluids remain in the upper layer. At moderate to large 
$Bo$, two behaviours are observed depending on the value of 
$m$: descending–steady-state or plateau regimes. If any liquid becomes trapped in the wake behind the rising bubble, the entrainment regime is identified as descending–steady-state (
$m\leq 1$); otherwise, it is denoted as the plateau regime (
${m>1}$).
The entrainment regimes for 
$\rho =0.7$ and 
$Y=0$ are shown in figure 16(a). In this case descending and descending–steady-state regimes occur. Due to buoyancy, the entrained liquid is driven to settle downward. When 
$m \leq 1$, the settlement is faster. At low 
$m$, the entrained liquid experiences less viscous shear stresses, also due to break up of the bubble upon entering the upper layer and the formation of small bubbles, it takes longer to reach a steady state 
$V_e$. We observe the descending regime herein. At large 
$m$, the entrained liquid experiences larger shear stresses in the upper layer. Combined with the negative buoyancy force 
$V_e$ reduces but reaches a plateau, i.e. descending–steady-state.
Figure 16. Classification of entrainment regimes over the 
$Ar$–
$Bo$ domain for (a) 
$\rho =0.7$ and 
$Y=0$ and (b) 
$\rho =0.7$ and 
$Y=0.1$. The size of symbols represents the value of 
$m$ at each point, i.e. 
$m={0.1}$, 1, 10 with magenta, black and cyan colours, respectively. Regimes are represented by symbols as follows: 
$\square$ plateau; 
$\diamond$ descending–steady-state; 
$\triangledown$ descending; and 
$\triangle$ ascending.
The entrainment regimes for 
$\rho =0.7$, 
$Y=0.1$ are shown in figure 16(b). In this case plateau and descending–steady-state regimes occur. The entrained liquid descends and settles down due to negative buoyancy forces that it experiences. The descending–steady-state regime occurs when Bingham liquid has been trapped within the wake behind the rising bubble. The plateau regime occurs when there is no liquid trapped in the wake of the bubble.
The plateau regime is of particular interest since the Bingham liquid that has been pulled up into the lighter Newtonian liquid does not settle. This suggests that the buoyancy stresses are not strong enough to surpass the yield stress of the material, hence it remains in the upper layer.
6. Maximum volume of entrainment 
$V_e$
These four regimes above primarily emphasise how far up from the liquid–liquid interface the bubbles can transport the lower layer fluid. In all regimes except for the ascending regime, the maximum volume of entrained liquid has been transported cumulatively. However, there is no direct link between the entrainment regimes and the maximum volume of entrainment.
It is of interest to identify the set of flow parameters for which most liquid is transported. We have estimated the entrainment volume by integrating over the domain above 
$z=0$, including entrained fluid. This value has been normalised with the volume of the rising bubble. Variation of the volume of entrained fluid as it rises, for 
$\rho =1$ and 
$\rho =0.7$, is shown in figures 17 and 18, respectively. When the upper layer and the lower layer have the same density (
$\rho =1$), the maximum 
$V_e$ occurs when both liquid layers are Newtonian. When the lower layer is a Bingham fluid, 
$V_e$ decreases significantly, such that a bubble transports at most its own volume.
Figure 17. The volume of entrained fluid (
$V_e$) against the position of the bubble, where the fluid layers have similar densities (
$\rho =1$).
Figure 18. The volume of entrained fluid (
$V_e$) against the position of the bubble, where the upper fluid layer has lower density than the lower fluid layer (
$\rho =0.7$).
In the Newtonian case, the amount of liquid that has been transported decreases mainly by reducing 
$m$, i.e. the maximum entrainment observed occurs when 
$m=10$ and it reduces by reducing 
$m$. By increasing 
$m$ from 0.1 to 10, the viscosity of the fluid in the upper layer increases. The higher viscosity causes the column of entrained liquid to break off earlier. However, a larger area of the liquid–liquid interface has been displaced from the horizontal line (
$z=0$). In addition, for 
$m<1$ and by increasing 
$Ar$ and 
$Bo$, a wake forms behind the rising bubble and some liquid becomes trapped there. The results show that the displaced fluid at the interface constitutes the more significant part of the entrained volume. At each 
$m$, intermediate 
$Bo$ values result in more entrainment and increasing 
$Ar$ decreases the entrainment.
In the Bingham case, when 
$m < 1$ the bubble transports more fluid from the lower layer in comparison with 
$m >1$. The maximum value is quickly attained and generally remains constant (plateau regime). At each 
$m$, decreasing 
$Bo$ values and increasing 
$Ar$ results in higher 
$V_e$. For 
$\rho =0.7$, similar to 
$\rho =1$, the maximum 
$V_e$ happens for the Newtonian case with 
$m=10$. Although the density difference causes an appreciable reduction in the displacement of the interface, it is still comparable to other values of 
$m$. Lastly, 
$V_e$ does not diminish by having a Bingham fluid in the lower layer and a light Newtonian fluid in the upper layer.
7. Summary
The dynamics of single bubbles rising from a viscoplastic fluid layer into a Newtonian fluid layer, and the displacement of fluid from the lower layer to the upper layer, have been studied computationally. The effect of density and (effective) viscosity differences between the fluid layers on the entrainment process has been investigated, and the results have been compared with the Newtonian counterpart. The bubble is initially positioned in a viscoplastic fluid below a Newtonian fluid. As the bubble rises up it can undergo breakup within the viscoplastic fluid layer, with a small portion of it remaining confined within this layer. This phenomenon is found to occur for 
$Ca \geqslant 5$. Our results show that bubbles take on a prolate shape when passing through the viscoplastic–Newtonian interface at 
$Bo >1$, and it is more pronounced at 
$m\leqslant 1$.
Bubbles transport material from the lower layer to the upper layer through various mechanisms: displacing the liquid–liquid interface, forming a towed column of the lower layer liquid and/or entrapping it within its wake. In cases where the lower layer is a Bingham fluid, the yield stress has been found to be sufficiently high to prevent the horizontal interface from being significantly displaced from its initial position as the bubble crosses it. This reduces the entrained volume. In general, the maximum 
$V_e$ is found when the lower layer is Newtonian.
As the bubble ascends through the upper layer, the entrainment volume undergoes changes over time, resulting in four distinct entrainment regimes: plateau, descending–steady-state, descending and ascending regimes. On a larger scale, these regimes may imply that rising bubbles could create mixing regions (of fluids from two layers) at different depths. In the plateau regime, the mixing region is situated directly above the liquid–liquid interface and has a thickness of around 5 radii. A small fraction of the entrained liquid is likely to remain attached to the bubble and travels to the free surface. In the ascending regimes, bubbles have the capacity to transport the entrained material over greater distances and form thicker mixing layers above the liquid–liquid interface. In the case of the Bingham lower layer, the practical implication on a large scale is that the entrained fluid could either remain above the interface, forming localised mixing regions in the plateau regime, or be transported higher and even reach the free surface.
It is found that the density difference of 
$70\,\%$ is not enough to prevent the entrainment of a viscoplastic fluid, 
$Y=0.1$. Although the entrainment decreases in comparison with 
$\rho =1$, there is still a possibility for formation of localised mixing regions of upper and lower fluids due to entrainment. Moreover, the entrained heavy viscoplastic fluid does not recede over time, i.e. 
$V_e$ for the case of 
$\rho =0.7$ and 
$Y=0.1$ does not change after the bubble tail pinches off. It is because the denser viscoplastic fluid can support itself while static.
Acknowledgements
We thank the reviewers for their critical and interesting comments.
Funding
This work was supported by the NSERC discovery (M.Z. and I.F., grant number RGPIN-2020-04471) and the NSERC grant (M.Z. and G.L., grant number AWD-001563). The computations described in this paper were carried out using the resources of the Digital Research Alliance of Canada. This support is gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
