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. Reference Li, Meijer and Lohse2022), which is a slightly modified version of the double-bucket method (Oster Reference Oster1965). 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. Reference Lin, Leger, Kunkel and McCarthy2013; Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019). 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. Reference Khattab, Bandarkar, Fakhree and Jouyban2012), 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. Reference Young, Goldstein and Block1959), 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. (Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019), 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 Reference Phillips1970; Wunsch Reference Wunsch1970), 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. Reference Li, Diddens, Prosperetti and Lohse2021). 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 Reference Diddens2017). 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. Reference Cairncross, Schunk, Baer, Rao and Sackinger2000), 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 Reference Brooks and Hughes1982). 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 Reference Heil and Hazel2006) 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. Reference Li, Meijer and Lohse2022), 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. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019, Reference Li, Diddens, Prosperetti and Lohse2021, Reference Li, Meijer and Lohse2022), 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. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019). Also, because of the enhanced drag caused by the stable stratification, the sinking velocity decreases exponentially (Zvirin & Chadwick Reference Zvirin and Chadwick1975; Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019). Later, this entrained layer breaks due to diffusion, and the Marangoni flow velocity increases exponentially (Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019), 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. Reference Yick, Torres, Peacock and Stocker2009); 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. Reference Li, Diddens, Prosperetti and Lohse2021, Reference Li, Meijer and Lohse2022), $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. (Reference Young, Goldstein and Block1959), 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. Reference Li, Diddens, Prosperetti and Lohse2021), 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. Reference Li, Meijer and Lohse2022). According to Young et al. (Reference Young, Goldstein and Block1959), 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 (Reference Zvirin and Chadwick1975), (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 Reference Zvirin and Chadwick1975; Li et al. Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019) 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 (Reference Zvirin and Chadwick1975). 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. Reference Young, Goldstein and Block1959) 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 Reference Hadamard1911; Rybczynski Reference Rybczynski1911). 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 (Reference Zvirin and Chadwick1975) 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. (Reference Candelier, Mehaddi and Vauquelin2014) 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. (Reference Mehaddi, Candelier and Mehlig2018) 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. Reference Srdić-Mitrović, Mohamed and Fernando1999; Torres et al. Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000; Yick et al. Reference Yick, Torres, Peacock and Stocker2009; Zhang et al. Reference Zhang, Mercier and Magnaudet2019), 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. (Reference Yick, Torres, Peacock and Stocker2009), 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. (Reference Zhang, Mercier and Magnaudet2019). 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. (Reference Candelier, Mehaddi and Vauquelin2014). Otherwise, if $l_{D} \ll l_{s} \ll l_{\nu }$ (Regime 2), then Zvirin & Chadwick (Reference Zvirin and Chadwick1975) 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. Reference Zhang, Mercier and Magnaudet2019). They argue that since the experimental and numerical results by Yick et al. (Reference Yick, Torres, Peacock and Stocker2009) 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. Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000; Yick et al. Reference Yick, Torres, Peacock and Stocker2009; Doostmohammadi et al. Reference Doostmohammadi, Dabiri and Ardekani2014; Zhang et al. Reference Zhang, Mercier and Magnaudet2019; Mandel et al. Reference Mandel, Waldrop, Theillard, Kleckner and Khatri2020). 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 (Reference Zvirin and Chadwick1975) for $Pe \rightarrow \infty$, and by Candelier et al. (Reference Candelier, Mehaddi and Vauquelin2014) 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.