On the rising and sinking motion of bouncing oil drops in strongly stratified liquids

Abstract

When an immiscible oil drop is immersed in a stably stratified ethanol–water mixture, the Marangoni flow on the surface of the drop can experience an oscillatory instability, so that the drop undergoes a transition from levitating to bouncing. The onset of the instability and its mechanisms have been studied previously (Li et al., Phys. Rev. Lett., vol. 126, 2021, 124502; Li et al., J. Fluid Mech., vol. 932, 2022, A11), yet the bouncing motion of the drop itself, which is a completely different problem, has not yet been investigated. Here we study how the bouncing characteristics (jumping height, rising and sinking time) depend on the control parameters (drop radius, stratification strength, drop viscosity). We first record experimentally the bouncing trajectories of drops of different viscosities in different stratifications. Then a simplified dynamical analysis is performed to get the scaling relations of the jumping height and the rising and sinking times. The rising and sinking time scales are found to depend on the drag coefficient $C_D^S$ of the drop in the stratified liquid, which is determined empirically for the current parameter space (Zhang et al., J. Fluid Mech., vol. 875, 2019, 622–656). For low-viscosity (5 cSt) oil drops, the results on the drag coefficient match those from the literature (Yick et al., J. Fluid Mech., vol. 632, 2009, pp. 49–68; Candelier et al., J. Fluid Mech., vol. 749, 2014, pp. 184–200). For high-viscosity (100 cSt) oil drops, the parameter space had not been explored and the drag coefficients are not readily available. Numerical simulations are therefore performed to provide external verification for the drag coefficients, which well match with the experimental results.

1. Introduction

When a drop is placed in a stably stratified liquid with a concentration gradient, a Marangoni flow is generated on the surface of the drop. The interplay between this Marangoni flow and gravity will make the drop levitate or bounce continuously (Li et al. 2019, 2021). The transition from levitating to bouncing is caused by an oscillatory Marangoni instability (Li et al. 2021), which has two different mechanisms depending on the drop viscosity (Li, Meijer & Lohse 2022). However, once the drop starts to bounce, the motion of the drop acts back on the velocity field, and the flow field around the drop becomes completely different, namely unstationary rather than stationary. Thus to calculate the motion of the drop itself becomes a completely new problem, as compared to calculating the onset of this instability. Though the bouncing cycle has been described briefly and qualitatively before (Li et al. 2019), many questions regarding the properties of the bouncing trajectory remain unanswered.

The problem is actually that of a drop moving inside a density stratification with the presence of Marangoni flow. This problem is very complicated since the Marangoni advection is coupled with the drop motion and diffusion.

Before diving further into this topic, it is beneficial to first review briefly the relatively simpler and more basic problem of a solid particle moving through a stratified medium, which is common in natural environments and of great interest in various scientific fields. One example is marine snow, which plays a central role in the marine carbon cycle, and understanding its delayed vertical motion due to stratification is essential for bio-geochemical processes (Prairie et al. 2013, 2015). Another example is the motion of aerosols in the stratified atmosphere, which is of significant importance to the Earth's climate system, since they scatter and absorb a considerable amount of radiation (Jacobson 1999; Huneeus, Chevallier & Boucher 2012). A third example is the formation and sinking of ice crystals in the stably stratified layering of the saturated salt water in the Dead Sea in winter (Burns & Meiburg 2012, 2015; Sutherland, Barrett & Gingras 2015), which is further complicated by the growth of the crystals during the sinking process. The complicated physics involved when a spherical particle moves through a (sharply) stratified medium drew the attention of fluid dynamicists, and the problem has since then been studied analytically (Zvirin & Chadwick 1975; Candelier, Mehaddi & Vauquelin 2014; Mehaddi, Candelier & Mehlig 2018), numerically (Torres et al. 2000; Doostmohammadi, Dabiri & Ardekani 2014; Lee, Fouxon & Lee 2019; Zhang, Mercier & Magnaudet 2019), experimentally (Srdić-Mitrović, Mohamed & Fernando 1999; Abaid et al. 2004; Hanazaki, Kashimoto & Okamura 2009; Camassa et al. 2022) or by a combination of the above (Yick et al. 2009; Camassa et al. 2010).

Coming back to the problem of a drop moving inside stratified liquids with the presence of Marangoni flow: apart from the motion of the particle, now a Marangoni flow is coupled as well. There has been some research on this topic. For example, Blanchette & Shapiro (2012) studied the motion of a denser drop inside a sharply stratified liquid, where the length scale of the stratification is smaller than the drop size. It was found that the settling drop could bounce up due to Marangoni flow. Mandel et al. (2020) studied the rising of lighter drops in a two-layer density stratification. Here, the length scale of the stratification is much larger than the drop size, but the potential Marangoni force is analysed only as a side effect, because in that case it is weak (at most comparable to the drop's buoyancy) and only assists the rising of the lighter drop, which will rise anyway due to buoyancy. Other research either ignored the density gradient of the surrounding medium (Young, Goldstein & Block 1959; Leven & Newman 1976; Chen & Stebe 1996) or did not consider the Marangoni effect (Bayareh et al. 2013; Shaik & Ardekani 2020). It is also worth mentioning that the effects of surfactants on drop motion have also been studied (Levich 1962; Leven & Newman 1976; Chen & Stebe 1996; Martin & Blanchette 2017). For a more extensive summary regarding the motion in stratified liquids, the reader is referred to the reviews by Magnaudet & Mercier (2020) and More & Ardekani (2022). For a review on further physicochemical hydrodynamical phenomena of drops, we refer to Lohse & Zhang (2020).

For the problem considered here, the Marangoni force plays a major role in determining the speed and direction of the drop, and the length scale of the stratification is much larger than the size of the drop. In order to understand how in this case the drop's trajectory changes with the physical properties, such as the drop viscosity, the drop size and the concentration gradient, experiments are performed on low (5 cSt) and high (100 cSt) viscosity silicone oil drops in various stratified ethanol–water mixtures. A simplified dynamical model is developed to help to understand the bouncing trajectories. The drag coefficient of the drop $C_D^S$ in the stratified fluid is found to be the key parameter to determine the rising and sinking time scales. We take its dependence on the parameters from literature. For the low-viscosity drops, the drag coefficient has been suggested by Yick et al. (2009) and Candelier et al. (2014). For the high-viscosity oil drops, there is no existing research, so numerical simulations are performed to provide independent verifications for the drag coefficient. The drag coefficients thus obtained are found to agree well with the experimental results.

The paper is organized as follows. In § 2, the experimental procedure and the numerical method are described. In § 3, the general characteristics of the bouncing trajectory are described. After that, in § 4, a simplified dynamical analysis is provided to derive predictions for the minimum and maximum bouncing positions, as well as the dominant time scale for the rising and sinking motion of the drop, in which the drag coefficient $C_D^S$ is found to be the key parameter. In § 5, we compare the theoretical predictions regarding the minimum and maximum bouncing positions to our experimental observations, finding good agreement. In § 6, we summarize briefly the main results discussed in the literature regarding the drag coefficient on spherical objects in stratified media. The results of the rising and sinking time scales of low and high viscosities are discussed in §§ 7 and 8, respectively. The paper ends with conclusions and an outlook in § 9.

2. Experimental procedure and numerical methods

2.1. Experimental procedure

A sketch of the experimental set-up is shown in figure 1(a). A cubic glass container (Hellma, 704.001-OG, Germany) with inner horizontal extension $L= 30$ mm contains the linearly stratified ethanol–water mixture. The mixture is prepared using a ‘double-syringe’ method (Li et al. 2022), which is a slightly modified version of the double-bucket method (Oster 1965). To avoid bubble formation during mixing, both ethanol (Boom B.V., 100 %(v/v), technical grade, the Netherlands) and Milli-Q water are degassed in a desiccator at $\sim$2000 Pa for 20 min before making the mixture. Two layers of uniform ethanol concentration are located at the bottom (weight fraction $w_{b}$) and at the top (weight fraction $w_{t}$), between which the ethanol weight fraction $w_{e}$ increases linearly; see figure 1(b). Immediately after the mixture is prepared, $w_{e}(y)$ is measured by laser deflection (Lin et al. 2013; Li et al. 2019). The two uniform layers $w_{b}$ and $w_{t}$ are used to increase the accuracy of this method. The density of the mixture $\rho (y)$ is calculated from $w_{e}(y)$ using an empirical equation (Khattab et al. 2012), and the height at which the density of the mixture $\rho (w_{e}')$ matches the density of the oil $\rho '$ is set as $y = 0$; see figure 1(c). The values of $w_{b}$ and $w_{t}$, as well as the height of the stratified layer, are varied to change the stratification strength $\mathrm {d}w_{e}/\mathrm {d}y$.

Figure 1. ($a$) Sketch of the experimental set-up. Using a modified version of the double-bucket method, a stable and linearly stratified ethanol–water mixture is generated in the middle of the container. Two liquid layers of uniform concentration are injected at the top ($w_{t}$) and at the bottom ($w_{b}$). The cubic glass container has an inner horizontal extension $L= 30$ mm. Silicone oil drops of varying radii $R$ and viscosities $\nu '$ are released from the top. ($b$) Ethanol weight fraction of the mixture $w_{e}$ as a function of height. ($c$) Density of the mixture $\rho$ as a function of height. The density of the mixture matches that of the drop $\rho '$ at $y=0$. This is called the density-matched position, where $\rho (w_{e}')=\rho '$. The height of the drop $h$ is measured with respect to this density-matched position.

Drops are released from the top layer using a 1 $\mathrm {\mu }$l syringe (Hamilton, KH7001) through an attached needle, whose outer diameter is 0.515 mm. They are released one at a time to ensure that only a single drop is present in the container. The properties of the different silicone oils (Sigma-Aldrich, Germany) are reported in table 1. Due to their robustness to surface contaminations (Young et al. 1959), which alternatively would alter the interfacial surface properties and hence the motion of the drop, silicone oil droplets form the ideal candidate for our study. The interfacial tensions $\sigma$ between both silicone oils and several ethanol–water mixtures are measured on a goniometer (OCA 15Pro, DataPhysics, Germany) by using the pendant-drop method; see Appendix A. The drop is illuminated by a collimated LED (Thorlabs, MWWHL4), and its motion is captured by a side-view camera (Nikon D850) connected to a long-working-distance lens (Thorlabs, MVL12X12Z plus 0.25X lens attachment). All images are recorded at 30 frames per second. After the drop has completed its third bouncing cycle, it is carefully taken out of the mixture using a second thin needle. In all the cases, the third bouncing cycle is used to study its dynamics. After this, another slightly smaller drop is released, whose motion is now captured. We repeat this process several times but for no longer than 40 min, as by that time diffusion will affect the linear stratification of the mixture. The stability of the background density gradient as a function of time has been discussed by Li et al. (2019), where it remained stable for more than an hour. If still more data are desired, a new mixture is generated and the entire process is repeated.

Table 1. Properties of the silicone oils used in the experiments.

As the drop moves though the liquid, its motion will cause perturbations in the flow field that extend over a typical length scale. Since this length scale depends on the stratification strength of the mixture (Phillips 1970; Wunsch 1970), it can happen that for weak density gradients it is comparable to or even larger than the actual size of the container. The finite size of the container might therefore affect the motion of the drop as the stratifications become weaker (Li et al. 2021). To exclude such container-size effects, experiments with weaker density stratification, i.e. $\mathrm {d}w_{e}/\mathrm {d}y<60\ {\rm m}^{-1}$, are performed in a larger container (Hellma, 704.003-OG, Germany) with an inner horizontal extension $L= 50$ mm.

2.2. Numerical simulation

Numerical simulations of the whole bouncing cycle are also performed to provide independent verification for the drag coefficients. The simulations utilize the numerical framework that has been developed to simulate the evaporation of multi-component drops (Diddens 2017). Below, we briefly discuss its details.

In order for the model to be compared to the experimental data, it has to account for all the relevant physical mechanisms during the bouncing process. In particular, these are the flow driven by Marangoni and buoyancy effects, and the advection and diffusion of the ethanol–water mixture outside the drop. Mass transfer across the drop's interface is not taken into account because the solubility of ethanol in silicone oil is negligible, or vice versa. The interface of the drop needs to be well-resolved to capture the Marangoni flow. A sharp-interface finite element method has been developed, where the mesh is always conforming with the moving interface. In addition, the mesh is treated as a pseudo-elastic body (Cairncross et al. 2000), so that the bulk nodes follow the motion of the interfacial nodes. With moving mesh nodes, the numerical approach belongs to the class of arbitrary Eulerian–Lagrangian methods, which furthermore require us to consider the nodal movement $\dot {\boldsymbol {R}}$ at the interface. The problem is solved in axisymmetric cylindrical coordinates.

The governing incompressible Navier–Stokes equations are

(2.1)$$\begin{gather} \rho_0 ( \partial _t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ) ={-} \boldsymbol{\nabla} p^{\phi} - \rho^{\phi} (w_{e} ) g\,\boldsymbol{j} + \boldsymbol{\nabla} \boldsymbol{\cdot} [\mu^{\phi} (w_{e} ) ( \boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u} )^{\rm T} ) ], \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$

where $\phi = {d}, {b}$ denotes the phase (i.e. the drop and the bulk liquid), respectively. Using the Boussinesq approximation, the composition-dependent mass density $\rho ^{b} (w_{e} )$ is considered only for the gravitational term in (2.1). For the inertial term, a constant density $\rho _0$ is assumed. The mass fraction $w_{w}^{b}$ for water is determined using the advection–diffusion equation

(2.3)\begin{equation} \partial _t w_{w}^{b} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} w_{w}^{b} = \boldsymbol{\nabla}\boldsymbol{\cdot}(D^{b}(w_{e})\,\boldsymbol{\nabla} w_{w}^{b}), \end{equation}

with $D^{b} (w_{e} )$ the composition-dependent diffusivity in the bulk. The ethanol content is determined from the remainder of the water content. To account for the high $Pe$ numbers due to the low diffusivity, the equations are stabilized by a streamline upwind Petrov–Galerkin method (Brooks & Hughes 1982). At the interface, the boundary conditions are

(2.4)$$\begin{gather} (\boldsymbol{u} - \dot{\boldsymbol{R}} ) \boldsymbol{\cdot} \boldsymbol{n} = 0, \end{gather}$$

(2.5)$$\begin{gather}\boldsymbol{\tau}^{d} \boldsymbol{\cdot} \boldsymbol{n} - \boldsymbol{\tau}^{b} \boldsymbol{\cdot} \boldsymbol{n} = \sigma (w_{e} )\,\kappa \boldsymbol{n} + \boldsymbol{\nabla}_{S} \sigma (w_{e} ). \end{gather}$$

Here, $\boldsymbol {n}$ is the unit interface normal pointing from the drop to the bulk liquid, $\boldsymbol {\tau }^{\phi }$ are the stress tensors in both phases, $\sigma (w_{e} )$ is the composition dependent surface tension, $\kappa$ is the curvature, and $\boldsymbol {\nabla }_{S}$ is the surface gradient. All composition-dependent physical properties of the bulk liquid are implemented by a best fit of their corresponding values on the ethanol weight fraction; see Appendix A. The above equations are implemented in the finite element framework OOMP-LIB (Heil & Hazel 2006) on triangular mixed first-order Lagrange elements for the composition, and conventional Taylor–Hood elements for the velocity and pressure. To prevent the mesh from deforming significantly, it is reconstructed whenever the mesh quality falls below a specific threshold. Identical simulations are repeated with different mesh and domain sizes to ensure that the obtained results are not affected significantly.

Results of the numerical simulations and their comparison to the experimental observations are discussed in more detail in § 8. The sections that follow first focus on the experiments.

3. General characteristics of the bouncing cycle

When an immiscible drop is placed in a stably stratified ethanol–water mixture, the ethanol concentration gradient leads to a surface tension gradient on the surface of the drop, which points downwards, and a downward Marangoni flow is generated as a consequence. When the Marangoni flow is stable, the drop levitates at a fixed height; otherwise, for larger drops, the Marangoni flow is oscillatory due to an oscillatory instability (Li et al. 2022), and the drop bounces continuously. Whereas our previous studies looked into the onset of the bouncing instability, what mechanism triggers it and how it depends on the viscosity of the oil (Li et al. 2019, 2021, 2022), here the dynamics of the bouncing cycle itself is studied. Figure 2($a$) shows the trajectories of two 5 cSt drops of different radii in a stable stratification with $\mathrm {d}w_\mathrm {e}/\mathrm {d}y=55\ {\rm m}^{-1}$. As can be seen, the drops bounce periodically. Before the background gradient is changed due to diffusion and the bouncing of the drop, each period of the bouncing trajectory is identical. Figure 2($d$) zooms in on one period of the trajectory of the $R=188\ \mathrm {\mu } {\rm m}$ drop. The drop first sinks towards its density-matched position (before 110 s) due to gravity. During this period, an entrained liquid layer with almost uniform concentration is dragged down with the drop, and the Marangoni flow on the drop is very weak (Li et al. 2019). Also, because of the enhanced drag caused by the stable stratification, the sinking velocity decreases exponentially (Zvirin & Chadwick 1975; Li et al. 2019). Later, this entrained layer breaks due to diffusion, and the Marangoni flow velocity increases exponentially (Li et al. 2019), leading to the sudden upward jump of the drop (after $\sim$110 s). Because of this exponential behaviour, the Marangoni flow velocity reaches its maximum value in a short time period. Later, as the drop moves up, the buoyancy force on it increases, so it decelerates, until a time when the drop's upward velocity is almost zero. At this time, the drop is not moving and the strong Marangoni flow mixes the liquid around the drop, thus greatly weakening the Marangoni flow itself, which will finally make the drop sink. During sinking, the drop first accelerates for a short period, then decelerates towards its density-matched position again. The drop repeats this bouncing cycle thereafter. For most of the time, the drop is decelerating, so that the rising trajectory curves upwards and the sinking trajectory curves downwards. This is part of the reason why the trajectory is asymmetric. The other reason is that the average rising velocity of the drop is larger than the sinking velocity. For drops of higher viscosity, the rising velocity is smaller, thus their trajectories are less asymmetric; see figures 2($a$–$c$).

Figure 2. ($a$,$b$) Experimentally obtained trajectories of 5 cSt silicone oil drops of different radii in two different linearly stratified ethanol–water mixtures, with $\mathrm {d}w_{e}/\mathrm {d}y \approx 55\ {\rm m}^{-1}$ and $\mathrm {d}w_{e}/\mathrm {d}y \approx 105\ {\rm m}^{-1}$, respectively. At $t = 0$ s, the tracking procedure is initialized. The density-matched position, i.e. the position at which $\rho = \rho '$, is at $h = 0$ mm. ($c$) Trajectories of two 100 cSt silicone oil drops of different radii inside a linearly stratified ethanol–water mixture, with $\mathrm {d}w_{e}/\mathrm {d}y \approx 60\ {\rm m}^{-1}$. ($d$–$f$) Velocity $\dot {h}$ and height $h$ as functions of time for a single bouncing cycle for the smaller drop inside the corresponding linearly stratified mixture. The definitions of some characteristics of the bouncing cycle are indicated in ($d$), namely the two extrema, $h_{top}$ and $h_{bot}$, as well as the definition of the rising and sinking time intervals, $\tau _{rise}$ and $\tau _{sink}$. ($g$–$i$) Experimentally measured profiles of the mass density $\rho (y)$ and the interfacial surface tension $\sigma (y)$ using the laser deflection technique.

The bouncing cycles are different for drops of different radii, as can be seen from figure 2($a$). They also depend on the degree of stratification and on the drop viscosities; see figures 2($b$,$e$) and 2($c$,$f$), respectively. We denote the highest position of the trajectory as $h_{top}$, and the lowest position as $h_{bot}$, the rising time as $\tau _{rise}$, and the sinking time as $\tau _{sink}$; see figure 2($d$). The bouncing cycles can be characterized by these four quantities. In addition, figures 2($g$–$i$) show the experimentally determined profiles of the mass density $\rho (y)$ and the interfacial surface tension $\sigma (y)$. Approaching the uniform top layer gives rise to the pronounced nonlinearity in the profile of the mass density; see figure 2($h$). Only small drops in very strong stratified liquids will reach this region; see the more detailed analysis in § 5. As mentioned above, the trajectories of the drops are highly asymmetric, where a fast rise is followed by a slow descent. Although both oil viscosities show the same feature, this asymmetry is more dramatic for 5 cSt than for 100 cSt oil drops. In the following sections, we analyse theoretically how the four characteristic quantities of the bouncing cycle ($h_{top}$, $h_{bot}$, $\tau _{rise}$, $\tau _{sink}$) vary with the drop radius $R$, the strength of the stratification $\mathrm {d}w_{e}/\mathrm {d}y$, and the drop viscosity $\mu ^\prime$. The resulting scaling theory is then compared with the experimental results of the 5 cSt drops in §§ 5 and 7, and with the experimental and numerical results of the 100 cSt drops in § 8.

4. Dynamical analysis of the drop motion

The acceleration of the drop is caused by the forces acting on the drop: gravity and buoyancy $F_{B}$, drag force $F_{D}$, and a propulsion $F_{M}$ caused by the Marangoni flow. Thus

(4.1)\begin{equation} m' \ddot{h} = F_{B} + F_{D} + F_{M}, \end{equation}

where $m'$ is the mass of the drop, and $\ddot {h} = {\mathrm {d}^2h}/{\mathrm {d}t^2}$ is the acceleration of the drop. The added mass force and the Basset history force are not taken into account because they are found to be negligible in the parameter space studied here (Yick et al. 2009); see also § 7. They would at most modify prefactors. The Reynolds number is

(4.2)\begin{equation} Re=\frac{\vert\dot{h}\vert\,R}{\nu}, \end{equation}

where $\vert \dot {h}\vert$ is the absolute value of the drop velocity, $\nu =\mu /\rho$ is the kinematic viscosity of the mixture at the height of the drop, with $\mu$ and $\rho$ respectively the viscosity and density of the mixture at the height of the drop.

In our case, the Reynolds number is small (see figure 5), thus the drag force can be written as $F_{D}=-{\rm \pi} C^S_D\,Re\,\mu R\dot {h}/2$, where $C^S_D$ is the drag coefficient of the drop in the stratified liquid, which will be discussed extensively in the following sections. The buoyancy force is $F_{B}=-V^\prime g(\rho ^\prime -\rho )$, where $V^\prime$ is the volume of the drop. The propulsion is actually the viscous force caused by the Marangoni flow, $F_{M}=k\mu V_{M}R$, where $\mu$ is the viscosity of the mixture at the position of the drop, $V_{M}$ is the Marangoni flow velocity, and $k$ is a prefactor to be determined.

Because the Marangoni flow is oscillatory for the bouncing drops (Li et al. 2021, 2022), $V_{M}$ is not constant. When $V_{M}$ is strong, the drop can reach its highest position $h_{top}$, at which height the density of the ethanol–water mixture is $\rho _{top}$; when $V_{M}$ is weak, the drop will reach its lowest position $h_{bot}$, at which height the density of the ethanol–water mixture is $\rho _{bot}$. By balancing the Marangoni force with the buoyancy force – i.e. when $F_{M} = -F_{B}$ – we can obtain the density differences when the drop is at its highest and lowest positions:

(4.3)\begin{equation} \rho^\prime-\rho_{{top/bot}}=\frac{k \mu V_{M}R}{V^\prime g}. \end{equation}

From Young et al. (1959), we know that for the case of infinitely large diffusivity, zero density gradient and constant viscosity $\mu$, the Marangoni flow velocity at the equator of the drop is

(4.4)\begin{equation} V_{M}\vert_{equator} ={-}\frac{1}{2}\,\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}}\,\frac{\mathrm{d}w_{e}}{\mathrm{d} y}\,\frac{R}{\mu + \mu'}, \end{equation}

where $\sigma$ is the interfacial tension between the drop and the mixture, and $\mu '$ is the viscosity of the drop. But in our case, the diffusivity is not zero. Marangoni advection tends to smooth the concentration gradient close to the drop (Li et al. 2021), thus weakening the Marangoni flow itself. The Marangoni flow at the highest position $V_{M, top}$ and the lowest position $V_{M, bot}$ can then be written as

(4.5a,b)\begin{align} V_{{M, top}} \vert_{equator}={-}p \boldsymbol{\cdot} \, \frac{1}{2}\,\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}}\,\frac{\mathrm{d}w_{e}}{\mathrm{d} y}\,\frac{R}{\mu + \mu'}, \quad V_{{M, bot}}\vert_{equator} ={-}q \boldsymbol{\cdot} \, \frac{1}{2}\,\frac{\mathrm{d}\sigma}{\mathrm{d}w_{e}}\,\frac{\mathrm{d}w_{e}}{\mathrm{d} y}\,\frac{R}{\mu + \mu'}, \end{align}

where $0< q< p<1$ are two prefactors to be determined. Since $p$ and $q$ represent the influence of Marangoni advection on the concentration field, we do not expect that they can be calculated beforehand. But we do expect them to vary with the viscosity of the drop. A higher drop viscosity normally leads to a weaker advection, thus the concentration field is less distorted by advection, so that the prefactor is larger. This trend has been confirmed by the levitating drops (Li et al. 2022). According to Young et al. (1959), substituting (4.5a,b) and the volume of the drop into (4.3), we obtain

(4.6a,b) \begin{equation} \rho^\prime-\rho_{top}={-}\alpha \boldsymbol{\cdot} \,\frac{3}{2}\,\frac{\mu}{\mu + \mu'}\,\frac{\mathrm{d}\sigma}{\mathrm{d}y}\,\frac{1}{gR},\quad \rho^\prime-\rho_{bot} ={-}\beta \boldsymbol{\cdot} \,\frac{3}{2}\,\frac{\mu}{\mu + \mu'}\,\frac{\mathrm{d}\sigma}{\mathrm{d}y}\,\frac{1}{gR}, \end{equation}

where $\alpha =kp$ and $\beta =kq$. Thus $0<\beta <\alpha < k$ are two prefactors to be determined. Before we compare (4.6a,b) with the experimental values in § 5, we first analyse the governing time scale of the bouncing intervals $\tau _{rise}$ and $\tau _{sink}$.

Because the Reynolds number in our experiments is very small (see figure 5), the relevant time scale is effectively the time scale of the inertia-free system. That is, the acceleration occurs much faster than the force balance and is thus negligible. The time scale is then given by $F_{B}+F_{D}+F_{M}=0$. For a linear gradient, $\rho =\rho ^\prime +h\,\mathrm {d}\rho /\mathrm {d}y$, we thus have $F_{B}=V^\prime gh\,\mathrm {d}\rho /\mathrm {d}y$. This yields

(4.7)\begin{equation} \dot{h} = \frac{b}{a}\,h + \frac{c}{a}, \end{equation}

where

(4.8a–c)\begin{equation} a = \frac{\rm \pi}{2}\,\mu R C^S_{D}\,Re, \quad b =V' g\,\frac{d\rho}{{\rm d}y}, \quad c = k\mu V_{M}. \end{equation}

Not only the drag coefficient $C^S_{D}$ but also the Reynolds number $Re$ depends on the position of the drop. However, according to Zvirin & Chadwick (1975), (4.7) can have an accuracy to first order when evaluating the drag coefficient $C^S_{D}$ and the Reynolds number $Re$ at the instant the drop reaches its peak velocity $\dot {h}_{p}$. Equation (4.7) implies an exponential behaviour (Zvirin & Chadwick 1975; Li et al. 2019) with a governing time scale for the drop to reach its equilibrium position through either rising or sinking, which is

(4.9)\begin{equation} \tau_{rise/sink}\sim\tau_1\approx{-}\frac{a}{b}\sim\frac{\nu_{p} C^S_{D,{p}}\,Re_{p}}{N^2_{p}R^2}, \end{equation}

where

(4.10)\begin{equation} N_{p}=\sqrt{-\frac{g}{\rho}\,\frac{\mathrm{d}\rho}{\mathrm{d}y}} \end{equation}

is the Brunt–Väisälä frequency. Note that (4.9) is consistent with the results of Zvirin & Chadwick (1975). To evaluate the quantities in (4.9), ethanol weight fractions at the positions where the drops reach their peak velocities (during rising and sinking) are used to obtain the corresponding density $\rho _{p}$, viscosity $\mu _{p}$, and interfacial tension $\sigma _{p}$ (see Appendix A for the concentration dependence of these properties). It has become apparent that the drag coefficient $C_D^S$ is of great importance for the overall dynamics. Therefore, we provide more details regarding the drag on spherical objects in stratified liquids in § 6, followed by a discussion on the experimental and numerical results in §§ 7 and 8.

5. The density differences at the maximum and minimum bouncing positions

In this section, we compare the theoretically predicted density differences at the maximum and minimum bouncing positions with the experimentally measured ones. The quantities on the right-hand sides of (4.6a,b) for 5 cSt drops are calculated (excluding the prefactors $\alpha$ and $\beta$) and plotted against the experimentally measured density differences $\rho ^\prime -\rho _{top}$ and $\rho ^\prime -\rho _{bot}$ in figures 3(a,b), respectively. The position-dependent material properties $\mu$ and $\sigma$ are evaluated at either the highest or lowest position, correspondingly. The results for 100 cSt are shown in figure 4.

Figure 3. ($a$) Maximum and ($b$) minimum bouncing height of 5 cSt silicone oil drops of different radii in linearly stratified ethanol–water mixtures with indicated stratification strengths in m$^{-1}$. The deviations from the linear trend at larger weight fraction gradients can be explained by the occurrence of ceiling effects, i.e. the drop approaching the upper region of constant density. The measured prefactors are $\alpha _{5cSt} = 0.33$ and $\beta _{5cSt} = 0.08$.

Figure 4. ($a$) Maximum and ($b$) minimum bouncing position of 100 cSt silicone oil drops of different radii in linearly stratified ethanol–water mixtures with indicated stratification strengths in m$^{-1}$. The deviations from the linear trend at larger weight fraction gradients can be explained by the occurrence of ceiling effects i.e. the drop approaching the upper region of the constant density. The measured prefactors are $\alpha _{100cSt}= 0.39$ and $\beta _{100cSt}= 0.20$.

In both cases, the density differences follow a linear trend initially, as predicted by (4.6a,b), but are accompanied by a considerable amount of scatter. The scatter originates from the nonlinear profiles of both the viscosity $\mu$ and surface tension $\sigma$ as functions of the ethanol weight fraction $w_{e}$; see Appendix A.

From this procedure, the prefactors $\alpha$ and $\beta$ can be determined and are found to be $\alpha _{5cSt}=0.33$, $\beta _{5cSt}=0.08$, $\alpha _{100cSt}=0.39$ and $\beta _{100cSt}=0.20$, respectively. As drops move higher, the density differences – or equivalently, the bouncing heights – saturate, because the drops are reaching the top uniform layer (see figure 1), which forms the ‘ceiling’ of the bouncing trajectory. Not surprisingly, only small drops in large concentration gradients can reach the ceiling.

Comparing the obtained prefactors, it appears that they depend on the viscosity of the oil. To rationalize this observation, we recall the origin of the prefactor $3/2$ (Young et al. 1959) on the right-hand sides of (4.6a,b), where an infinitely large diffusivity was assumed. Given the fact that the Marangoni advection is important in our case, this assumption no longer holds, causing a reduction of this prefactor, hence $0<\beta <\alpha <1$. And since the Marangoni flow is inversely proportional to the viscosity of the oil $\mu '$ (see (4.4)), $\alpha$ and $\beta$ will increase with $\mu '$. Finding the exact functional form of $\alpha$ and $\beta$ on $\mu '$ is beyond the scope of the present paper.

6. Drag on spherical objects in stratified media

As discussed in § 4, the governing time scale of the motion of the drop (4.9) depends on the drag coefficient of the drops in the stratified liquid $C_D^S$. As mentioned in the Introduction, there are only a few studies on the topic of a drop moving inside stratified liquids with the presence of Marangoni flow. In particular, there are no available data on the drag coefficient of a drop moving inside such stratifications. However, we think the drag coefficient of a solid sphere in such stratifications could be used here, since we are interested only in the scaling of the rising and sinking times, and the drag coefficient of drops and solid spheres normally only differ by a prefactor (Hadamard 1911; Rybczynski 1911). This is confirmed later, in § 7. In this section, we will therefore provide a brief overview on the current understanding of drag on solid spheres in stratified liquids.

Analytically, the settling of a spherical body in a stratified liquid has been studied in the past by several authors. The relevant parameters for this problem are the Froude number, defined as

(6.1)\begin{equation} {Fr} = \frac{\vert \dot{h} \vert}{NR}, \end{equation}

the Péclet number, defined as

(6.2)\begin{equation} Pe = \frac{\vert \dot{h} \vert\,R}{D}, \end{equation}

and the Richardson number, defined as

(6.3)\begin{equation} Ri = \frac{Re}{{Fr}^2} = \frac{N^2R^3}{\nu\,\vert \dot{h} \vert}, \end{equation}

where $D$ is the solute diffusivity. For example, Zvirin & Chadwick (1975) derived that in the limit of $Re\ll 1$ and ${Fr}\ll 1$, and under the assumption that advection dominates, i.e. $Pe \rightarrow \infty$, the drag coefficient scales as

(6.4)\begin{equation} C_D^S \sim Ri^{1/3}/Re. \end{equation}

In the opposite limit, when diffusion dominates advection, i.e. when $Pe \ll 1$, Candelier et al. (2014) derived that

(6.5)\begin{equation} C_D^S \sim (Pe\,Ri)^{1/4}/Re. \end{equation}

It has been shown recently by Mehaddi et al. (2018) that both of these scaling relations can be obtained from the very same derivation, where the expressions above are simply limiting cases.

The drag coefficient of a particle settling in a density stratification has also been investigated experimentally and numerically (Srdić-Mitrović et al. 1999; Torres et al. 2000; Yick et al. 2009; Zhang et al. 2019), and different forms of the drag coefficient have been suggested in different parameter ranges. Some of the studied parameter ranges in terms of Froude number versus Reynolds number are summarized in figure 5. The figure also shows the parameter range of the present study: see the two red boxes for 5 and 100 cSt drops, respectively. It can be seen that the parameter range of the 5 cSt drops overlaps almost entirely with that of Yick et al. (2009), where the drag coefficient was determined empirically as

(6.6)\begin{equation} C_D^S \sim \frac{Ri^{0.51 \pm 0.11}}{Re}. \end{equation}

An argument for the discrepancy between this result and the analytical solutions was provided by Zhang et al. (2019). In their work, they defined three different stratification regimes depending on the relative magnitudes of three length scales: the viscosity length scale $l_{\nu }$, the diffusivity length scale $l_{D}$, and the stratification length scale $l_s$. If $l_s \ll l_{D} \ll l_{\nu }$ (Regime 1), then the drag coefficient scales as predicted by Candelier et al. (2014). Otherwise, if $l_{D} \ll l_{s} \ll l_{\nu }$ (Regime 2), then Zvirin & Chadwick (1975) gave the correct prediction. If $l_{D} \ll l_{\nu } \ll l_{s}$ (Regime 3), then $C_D^S \sim ({Fr}\,Re)^{-1}$ (Zhang et al. 2019). They argue that since the experimental and numerical results by Yick et al. (2009) mix two of these asymptotic regimes, the obtained scaling relation is an ad hoc approximation.

Figure 5. Black boxes are parameter spaces studied previously on the settling of spherical objects in linearly stratified liquids (Torres et al. 2000; Yick et al. 2009; Doostmohammadi et al. 2014; Zhang et al. 2019; Mandel et al. 2020). Both numerical and experimental results are included. The red boxes are parameter spaces of the present study in terms of $Re_{p}$ and ${Fr}_{p}$, for 5 and 100 cSt silicone oil drops, respectively. The limit of $Re\ll 1$ and ${Fr}\ll 1$ was studied analytically by Zvirin & Chadwick (1975) for $Pe \rightarrow \infty$, and by Candelier et al. (2014) for the opposite limit, $Pe \ll 1$.

Since we do not know the appropriate scaling relation of the drag coefficient a priori, especially for the 100 cSt oil drops that fall in a parameter regime not yet studied in the available literature (see figure 5), we assume a general expression as

(6.7)\begin{equation} C_D^S \sim Ri^q / Re, \end{equation}

where $q$ is an exponent to be determined. The scaling relation of the dominant time scale, (4.9), can then be rewritten as

(6.8)\begin{equation} \tau_{{rise/sink}} \sim \frac{R}{\vert \dot{h} \vert_{{p,rise/sink}}}\,Ri^{(q-1)}_{p}. \end{equation}

In the next section, (6.8) will be evaluated using the experimentally determined time intervals and drop velocities for 5 and 100 cSt oil drops, respectively. The obtained values of $q$ are then compared to the above-mentioned scaling relations (§ 7) and to our numerical simulations (§ 8).

In the following sections, we discuss first the experimental results for the low-viscosity drops, then the experimental and numerical results for the high-viscosity drops.

7. Results for the low-viscosity drops: 5 cSt

The result after evaluating (6.8) with the experimentally determined time intervals and drop velocities for 5 cSt oil drops during rising and sinking are shown in figure 6($a$,$b$), respectively. The solid lines show the best fits through the experimental data. It follows that since $q-1 = -0.53$ for rising and $q-1 = -0.75$ for sinking, a scaling relation of the drag coefficient during rising and sinking is obtained as

(7.1a,b)\begin{equation} C^S_{D,rise,5cSt} \sim \frac{Ri^{0.47}_{p}}{Re_{p}},\quad C^S_{D,sink,5cSt} \sim \frac{Ri^{0.25}_{p}}{Re_{p}}. \end{equation}

The scaling relations are clearly different. For the rising drop, the determined scaling relation is in good agreement with the empirical result of Yick et al. (2009), which we write as $C_D^S\sim Ri^{0.51\pm 0.11}/Re$; see (6.6). Although the latter was established for a solid particle, here it is shown that the same relation also applies to drops. The Marangoni velocity $V_{M}$ does not influence the governing time scale for rising or sinking (see (4.7)) and will influence the drag coefficient $C_D^S$ only by changing the parameter range (${Fr}, Re$). Thus as long as the parameter range fits, the results of the drag coefficient $C_D^S$ on particles can be applied directly to drops, which is in line with the observation that the studies overlap almost entirely in the parameter space shown in figure 5.

Figure 6. Non-dimensionalized ($a$) rising and ($b$) sinking times of 5 cSt silicone oil drops of different radii in linearly stratified ethanol–water mixtures with indicated stratification strengths as functions of $Ri_{p}$ on a double logarithmic scale. The solid line shows the best fit through the experimental data.

For the sinking drop, these scaling relations do not match. The obtained result seems to be in agreement with the scaling relation as predicted by Candelier et al. (2014), $C_D^S \sim (Pe\, Ri)^{1/4}/Re$, which is the asymptotic limit of the viscous-diffusive regime where diffusion dominates over advection. This is because, as discussed in § 3, during sinking, the drop is surrounded by an entrained shell where diffusion is dominant (Li et al. 2019). The fact that $C^S_{D,{rise}} > C^S_{D,{sink}}$ is the second reason mentioned earlier why the bouncing trajectory is asymmetric.

8. Results for the high-viscosity drops: 100 cSt

Before evaluating (6.8) with the experimentally measured values, it is important to realize that the parameter space spanned by the high-viscosity oil is not covered by any available literature; see figure 5. Although analytical results do exist for limiting cases, as discussed in § 6, they – as well as other similar studies (Lee et al. 2019; Dandekar, Shaik & Ardekani 2020; Shaik & Ardekani 2020) – relied on the crucial assumption $Ri \ll 1$, i.e. weakly stratified liquids, due to the applied mathematical method of asymptotic matching. Since in our case for the 100 cSt drops $1.7 \leq Ri_{p} \leq 27.6$, the liquid is strongly stratified, and it is questionable whether the above scaling relations apply to our present study. A new theoretical derivation of the drag coefficient of drops settling in a strongly stratified liquid including Marangoni effects is beyond the scope of the present work. Thus numerical simulations covering the same parameter range are performed to provide independent verification. Details of the numerical simulations have been provided in § 2.2. Here, we discuss the numerical results and compare them with the experiments.

8.1. Numerical results

The numerical simulations are initialized by placing a 100 cSt drop ($250\ \mathrm {\mu } {\rm m} \leq R \leq 400\ \mathrm {\mu } {\rm m}$) above its density-matched position, i.e. $h = 0$, into a linearly stratified ethanol–water mixture ($20\ {\rm m}^{-1} \leq \mathrm {d}w_{e}/\mathrm {d}y \leq 140\ {\rm m}^{-1}$), with corresponding uniform top and bottom layers. As the drop descends towards this equilibrium position, at a certain moment in time, Marangoni effects start to become dominant and the bouncing is initiated, the same as in the experiments. Since the simulations demand a rather small time step due to the high Péclet number, and thus considerable CPU time, only the first complete bouncing cycle is considered when analysing the numerical results. This is in contrast with the experiments, where the third bouncing cycle is used.

Typical snapshots of the numerical simulation as the drop reaches its peak velocity during rising and sinking are shown in figures 7($a$,$b$), respectively. We show the mass density of the mixture and how mixing occurs in close proximity to the drop, as well as the velocity magnitude inside and outside the drop. In the latter, clear vortical structures arise that are induced by the baroclinic torque following the deflection of the isopycnals. As the drop rises, a strong Marangoni flow is visible at the interface, causing strong internal circulations in the drop; see figure 7($a$). The isopycnals are compressed below the drop as lighter fluid is advected downwards by the Marangoni flow. During sinking, the high-velocity region has shifted towards the apex of the drop, and the strong internal circulation in the drop has stopped; see figure 7($b$). This leads to the conclusion that during sinking, Marangoni effects are indeed rendered ineffective; see Appendix B. The deflection of the isopycnals has become more significant as the drop is dragging liquid down.

Figure 7. Snapshot of the numerical simulation showing the mass density of the mixture (left) and the velocity magnitude (right) around an ($a$) rising and ($b$) sinking 100 cSt oil drop as it reaches its peak velocity. ($c$–$e$) Temporal evolution of the drop's height $h(t)$ with $R = 250\ \mathrm {\mu } {\rm m}$ for three different stratifications.

In addition, the trajectories of three bouncing drops with the same size but in different stratifications are plotted in figures 7($c$–$e$). Qualitatively, the numerically simulated bouncing cycles show great resemblance to the experimental observations; see figure 2. A quick rise is followed by a slower descent, where the asymmetry in the bouncing cycle is less dramatic, as for high-viscosity oils. Also, the bouncing period seems to shorten with stronger stratifications, and the bouncing amplitude becomes smaller, the same as in the experiments. The second rise in figure 7($c$) for a small drop in a relatively weakly stratified ethanol–water mixture has also been observed experimentally – see figure 2($d$) – although less pronounced.

8.2. Experimental and numerical comparison

Now we are in a position that allows us to compare the experimental and numerical results one-on-one. For this specific case, a 100 cSt oil drop with $R = 280\ \mathrm {\mu } {\rm m}$ is placed inside a stratified ethanol–water mixture where $\mathrm {d}w_{e}/\mathrm {d}y = 35\ {\rm m}^{-1}$. Figures 8($a$,$b$) compare the drop's trajectories and velocity profiles, respectively. Qualitatively, again, the agreement is good. Quantitatively, there are some differences. For example, the overall bouncing period is shorter in the simulations ($T = 52.3$ s) compared to the experiments ($T = 88.8$ s), and the peak velocities reached in the numerics ($\dot {h}_{{p,rise,num}} = 0.24$ mm s$^{-1}$ and $\dot {h}_{{p,sink,num}} = -0.13$ mm s$^{-1}$) exceed those of the experiment ($\dot {h}_{{p,rise,exp}} = 0.08$ mm s$^{-1}$ and $\dot {h}_{{p,sink,exp}} = -0.06$ mm s$^{-1}$). In addition, the top and bottom positions of the bouncing curve from the simulation ($h_{{top,num}} = 2.3$ mm and $h_{{bot,num}} = 1.3$ mm) are shifted slightly above the experimentally determined ones ($h_{{top,exp}} = 1.9$ mm and $h_{{bot,exp}} = 1.1$ mm). Both discrepancies might be explained by the fact that the Marangoni force, caused by Marangoni advection at the interface, is overestimated in the simulations. This would give rise to a quicker ascent, causing the drop to reach greater heights and thus reaching larger velocities during its descent, as the drop is further away from its density-matched position. Comparing the background concentration profiles and the corresponding concentration gradients in figure 8($c$) shows that although the concentration profiles look very similar, the local concentration gradient varies in the experiment, whereas it remains constant in the numerical simulation. Consequently, the concentration gradient $\mathrm {d}w_{e}/\mathrm {d}y$ that the drop ‘feels’ at $h_{bot}$ is smaller in the experiments than in the numerics. Since $V_{M} \sim \mathrm {d}w_{e}/\mathrm {d}y$ (see (4.4)), Marangoni advection will be slightly larger in the numerical simulation.

Figure 8. Experimentally measured (black) and numerically determined (red) ($a$) trajectory $h(t)$ and ($b$) velocity $\dot {h}(t)$ of a 100 cSt oil drop with $R = 280\ \mathrm {\mu } {\rm m}$ and $\mathrm {d}w_{e}/\mathrm {d}y = 35\ {\rm m}^{-1}$. (See supplementary movie 1, available at https://doi.org/10.1017/jfm.2023.415, for more details.) ($c$) The background concentration profiles and corresponding gradients.

Finally, (6.8) is evaluated with the experimentally and numerically determined time intervals and drop velocities for 100 cSt oil drops during rising and sinking. As in § 7, the aim is to determine the exponent $q$ to obtain the scaling relation of the drag coefficient $C_D^S \sim Ri^q / Re$. In figures 9($a$,$b$), the results are shown for both rising and sinking, respectively. The empty symbols represent the experimental results, and filled circles represent the numerical results.

Figure 9. Non-dimensionalized ($a$) rising and ($b$) sinking times of 100 cSt silicone oil drops of different radii in linearly stratified ethanol–water mixtures with indicated stratification strengths as functions of $Ri_{p}$ on a double logarithmic scale. The black (red) line shows the best fit through the experimental (numerical) data.

Based on the best fit through the data, indicated by the solid lines, for the experiments it holds that

(8.1a,b)\begin{equation} C^{S,{exp}}_{D,{rise,100cSt}} \sim \frac{Ri^{0.67}_{p}}{Re_{p}}, \quad C^{S,{exp}}_{D,{sink,100cSt}} \sim \frac{Ri^{0.68}_{p}}{Re_{p}}, \end{equation}

and for the simulations that

(8.2a,b)\begin{equation} C^{S,{num}}_{D,{rise,100cSt}} \sim \frac{Ri^{0.65}_{p}}{Re_{p}},\quad C^{S,{num}}_{D,{sink,100cSt}} \sim \frac{Ri^{0.54}_{p}}{Re_{p}}. \end{equation}

Considering these results, our study suggests that the drag coefficient of the high-viscosity drops in strongly stratified liquids scales as

(8.3a,b)\begin{equation} C^S_{D,{rise,100cSt}} \sim \frac{Ri^{0.66 \pm 0.01}}{Re}, \quad C^S_{D,{sink,100cSt}} \sim \frac{Ri^{0.61 \pm 0.07}}{Re}. \end{equation}

It is found that the scaling relations differ from the analytical predictions by Zvirin & Chadwick (1975) and Candelier et al. (2014), as well as from the empirical relation obtained by Yick et al. (2009).

An argument for this discrepancy has already been addressed in § 6. There, three different stratification regimes are introduced, depending on the relative magnitudes of the viscosity length scale $l_{\nu }$, the diffusivity length scale $l_{\mathrm {\kappa }}$, and the stratification length scale $l_s$ (Zhang et al. 2019). As mentioned, the experimental and numerical results by Yick et al. (2009) mix two of these asymptotic regimes, and the obtained scaling relation is therefore an ad hoc approximation. Provided that $Re \ll 1$ and $Pe \gg 1$, where $Pe = Re\,Pr$ with $Pr = \nu /D$, the domain of existence for each regime, respectively, is (Zhang et al. 2019)

(8.4a–c)\begin{equation} {Fr} \ll Re^{1/3}\,Pr^{{-}1/6}, \quad Pr^{{-}1/2} \ll {Fr} \ll Re^{{-}1}, \quad {Fr} \gg Re^{{-}1}. \end{equation}

Evaluating these conditions, it becomes evident that in the present study of the 100 cSt drops, the conditions of Regimes 1 and 2 are met simultaneously, and that a mix of these regimes occurs, hence the deviation from the theoretical predictions by Zvirin & Chadwick (1975).

9. Conclusion and outlook

To summarize, the bouncing dynamics of drops of different viscosities in stably stratified liquids with the presence of Marangoni flow is studied theoretically, experimentally and numerically. The main characteristics of the bouncing cycle have been discussed, and the importance of the findings by Young et al. (1959) regarding the prediction of the maximum and minimum bouncing positions has become evident. Based on our derivation of the scaling relation of the governing time scale in which the drop reaches its equilibrium position through either rising or sinking, it has become apparent that the scaling relation of the drag coefficient in the stratified liquid $C_D^S$ is of great importance. To this end, the experimentally determined quantities of the rising and sinking times, as well as the peak velocities reached during rising and sinking, are used to obtain the appropriate scaling relation of $C_D^S$ for drops in strongly stratified liquids. For low-viscosity oil drops, it was found that the drag coefficient follows the scaling relations obtained from literature. The significant difference between the relations obtained for rising and sinking explains the high asymmetry of the bouncing cycle. For the high-viscosity oil drops, the scaling relation of the drag coefficient is not available in the literature. Thus, to seek independent verification, numerical simulations are performed, mimicking the experiments for $1.7 \leq Ri_{p} \leq 27.6$. This also allowed for a one-to-one comparison between experiments and numerics, where it has been found that qualitatively, the agreement between the bouncing cycles is good. For both results, experimental and numerical, scaling relations are obtained for the high-viscosity oil drops during rising and sinking in strongly stratified liquids.

It is also found that when in the same parameter range, the scaling of the drag coefficient of a solid sphere could be applied to that of a drop, with or without Marangoni flow. This is supported by the results for 5 cSt drops. Thus the extensive knowledge on drag coefficients of solid spheres in stratified media can be of help to that of drops, which has rarely been explored. We also found that in the parameter space of the 100 cSt drops, the drag coefficient under stratified conditions has a scaling $C_D^S\sim Ri^{0.66\pm 0.01}/Re$.

As our work has shown, new insight has still to be discovered on the rising and sinking drops and bubbles in strongly stratified liquids. In particular, the dominant mechanism behind the drag enhancement, discovered only recently by Zhang et al. (2019), shows the complexity and richness of such hydrodynamical systems. Whereas the drag coefficient of spherical objects in homogeneous media has been investigated extensively for more than half a century (Clift, Grace & Weber 2005), it might be worthwhile to extend this research towards (strongly) stratified liquids.