## 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.* Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019, Reference Li, Diddens, Prosperetti and Lohse2021). The transition from levitating to bouncing is caused by an oscillatory Marangoni instability (Li *et al.* Reference Li, Diddens, Prosperetti and Lohse2021), which has two different mechanisms depending on the drop viscosity (Li, Meijer & Lohse Reference Li, Meijer and Lohse2022). 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.* Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019), 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.* Reference Prairie, Ziervogel, Arnosti, Camassa, Falcon, Khatri, McLaughlin, White and Yu2013, Reference Prairie, Ziervogel, Camassa, McLaughlin, White, Dewald and Arnosti2015). 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 Reference Jacobson1999; Huneeus, Chevallier & Boucher Reference Huneeus, Chevallier and Boucher2012). 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 Reference Burns and Meiburg2012, Reference Burns and Meiburg2015; Sutherland, Barrett & Gingras Reference Sutherland, Barrett and Gingras2015), 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 Reference Zvirin and Chadwick1975; Candelier, Mehaddi & Vauquelin Reference Candelier, Mehaddi and Vauquelin2014; Mehaddi, Candelier & Mehlig Reference Mehaddi, Candelier and Mehlig2018), numerically (Torres *et al.* Reference Torres, Hanazaki, Ochoa, Castillo and van Woert2000; Doostmohammadi, Dabiri & Ardekani Reference Doostmohammadi, Dabiri and Ardekani2014; Lee, Fouxon & Lee Reference Lee, Fouxon and Lee2019; Zhang, Mercier & Magnaudet Reference Zhang, Mercier and Magnaudet2019), experimentally (Srdić-Mitrović, Mohamed & Fernando Reference Srdić-Mitrović, Mohamed and Fernando1999; Abaid *et al.* Reference Abaid, Adalsteinsson, Agyapong and McLaughlin2004; Hanazaki, Kashimoto & Okamura Reference Hanazaki, Kashimoto and Okamura2009; Camassa *et al.* Reference Camassa, Ding, McLaughlin, Overman, Parker and Vaidya2022) or by a combination of the above (Yick *et al.* Reference Yick, Torres, Peacock and Stocker2009; Camassa *et al.* Reference Camassa, Falcon, Lin, McLaughlin and Mykins2010).

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 (Reference Blanchette and Shapiro2012) 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.* (Reference Mandel, Waldrop, Theillard, Kleckner and Khatri2020) 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 Reference Young, Goldstein and Block1959; Leven & Newman Reference Leven and Newman1976; Chen & Stebe Reference Chen and Stebe1996) or did not consider the Marangoni effect (Bayareh *et al.* Reference Bayareh, Doostmohammadi, Dabiri and Ardekani2013; Shaik & Ardekani Reference Shaik and Ardekani2020). It is also worth mentioning that the effects of surfactants on drop motion have also been studied (Levich Reference Levich1962; Leven & Newman Reference Leven and Newman1976; Chen & Stebe Reference Chen and Stebe1996; Martin & Blanchette Reference Martin and Blanchette2017). For a more extensive summary regarding the motion in stratified liquids, the reader is referred to the reviews by Magnaudet & Mercier (Reference Magnaudet and Mercier2020) and More & Ardekani (Reference More and Ardekani2022). For a review on further physicochemical hydrodynamical phenomena of drops, we refer to Lohse & Zhang (Reference Lohse and Zhang2020).

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.* (Reference Yick, Torres, Peacock and Stocker2009) and Candelier *et al.* (Reference Candelier, Mehaddi and Vauquelin2014). 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.* 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$.

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.

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

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

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

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$).

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

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

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:

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

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

*a*,

*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.5*a*,*b*) and the volume of the drop into (4.3), we obtain

*a*,

*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.6*a*,*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

where

*a*–

*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

where

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.6*a*,*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.

In both cases, the density differences follow a linear trend initially, as predicted by (4.6*a*,*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.6*a*,*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

the Péclet number, defined as

and the Richardson number, defined as

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

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

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

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.

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

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

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

*a*,

*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.* (Reference Yick, Torres, Peacock and Stocker2009), 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.

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.* (Reference Candelier, Mehaddi and Vauquelin2014), $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.* Reference Li, Diddens, Prosperetti, Chong, Zhang and Lohse2019). 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.* Reference Lee, Fouxon and Lee2019; Dandekar, Shaik & Ardekani Reference Dandekar, Shaik and Ardekani2020; Shaik & Ardekani Reference Shaik and Ardekani2020) – 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.

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.

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.

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

*a*,

*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

*a*,

*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

*a*,

*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 (Reference Zvirin and Chadwick1975) and Candelier *et al.* (Reference Candelier, Mehaddi and Vauquelin2014), as well as from the empirical relation obtained by Yick *et al.* (Reference Yick, Torres, Peacock and Stocker2009).

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.* Reference Zhang, Mercier and Magnaudet2019). As mentioned, the experimental and numerical results by Yick *et al.* (Reference Yick, Torres, Peacock and Stocker2009) 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.* Reference Zhang, Mercier and Magnaudet2019)

*a*–

*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 (Reference Zvirin and Chadwick1975).

## 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.* (Reference Young, Goldstein and Block1959) 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.* (Reference Zhang, Mercier and Magnaudet2019), 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 Reference Clift, Grace and Weber2005), it might be worthwhile to extend this research towards (strongly) stratified liquids.

## Supplementary movies

A supplementary movie is available at https://doi.org/10.1017/jfm.2023.415.

## Funding

We acknowledge support from the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), a NWO Gravitation programme funded by the Ministry of Education, Culture and Science of the government of the Netherlands, an ERC-Advanced Grant under project no. 740479, the Balzan Foundation. Y.L. acknowledges financial support from the Fundamental Research Funds for the Central Universities and the Natural Science Foundation of China under grant no. 12272376. C.D. kindly acknowledges financial support by the Industrial Partnership Programme (IPP) of the Netherlands Organization for Scientific Research (NWO). This research programme is co-financed by Canon Production Printing Holding B.V., University of Twente and Eindhoven University of Technology.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Physical properties of the ethanol–water mixture

Physical properties of the ethanol–water mixtures for different ethanol weight fractions are taken from Khattab *et al.* (Reference Khattab, Bandarkar, Fakhree and Jouyban2012) and Par *et al.* (Reference Par, Guevara-Carrion, Hasse and Vrabec2013); see figures 10($a$–$c$). The interfacial surface tension $\sigma (w_{e})$ between the silicone oil and the ethanol–water mixture is measured using the pendant-drop method on a goniometer and shown in figure 10($d$). The markers indicate the average value of six measurements, where the standard deviation is of the order of the size of the markers.

## Appendix B. Contours of the velocity magnitude

Figure 11 shows the contours of the velocity magnitude inside and outside the drop as the drop reaches its peak velocity during rising (left) and sinking (right). The Marangoni flow is strong as the drop rises, causing the contours to be closely packed at the interface of the drop. Additionally, a strong internal circulation inside the drop is visible. On the other hand, during sinking, the Marangoni flow is very weak, and the contour lines are not as closely packed. It can be seen that the internal circulation has vanished.