Retention of rising droplets in density stratiﬁcation

In this study, we present results from experiments on the retention of single oil droplets rising through a two-layer density stratiﬁcation, with the goal of quantifying and parametrizing the impact of stratiﬁcation on timescales that describe the delay in rising. These experiments conﬁrm the signiﬁcant slowdown observed in past literature of settling and rising particles and droplets in stratiﬁcation, and these are the ﬁrst experiments to study single liquid droplets as opposed to solid particles or bubbles. By tracking the motion of the droplets as they rise through a stratiﬁed ﬂuid, we identify two new timescales which quantitatively describe this slowdown: an entrainment timescale and a retention timescale. These timescales measure dynamics that were not captured in previous timescale discussions, which primarily focused on the timescale to the velocity minimum ( U min ). The entrainment timescale is a measure of the time that a droplet spends below its upper-layer terminal velocity and relates to the duration over which the droplet’s rise is affected by entrained dense ﬂuid. The retention time is a measure of the time that the droplet is delayed from reaching an upper threshold far from the density transition. These two timescales are interconnected by the magnitude of the slowdown ( U u − U min ) relative to the upper-layer terminal velocity ( U u ), as well as a constant that reﬂects the approximately universal form of the recovery of a droplet’s velocity from U min to U u . Both timescales are found to depend on the Froude and Reynolds numbers of the system, Fr = U u / ( Nd ) and Re = ρ u U u d /ν . We ﬁnd that both timescales are only signiﬁcantly large for Fr (cid:2) 1, indicating that trapping dynamics in a relatively sharp stratiﬁcation arise from a balance between drop inertia and buoyancy. Finally, we present a theoretical formulation for the force enhancement (cid:4) , the ratio between the maximum stratiﬁcation-induced force and the corresponding drag force on the droplet, based on a simple force balance at the point of the velocity minimum. Using our experimental data, we ﬁnd that our formulation compares well with recent theoretical and computational work by Zhang et al. [J. Fluid Mech. 875 , 622 (2019)] on the force enhancement on a solid sphere settling in a stratiﬁed ﬂuid, and provides the ﬁrst experimental data supporting their approach.

In this study, we present results from experiments on the retention of single oil droplets rising through a two-layer density stratification, with the goal of quantifying and parametrizing the impact of stratification on timescales that describe the delay in rising. These experiments confirm the significant slowdown observed in past literature of settling and rising particles and droplets in stratification, and these are the first experiments to study single liquid droplets as opposed to solid particles or bubbles. By tracking the motion of the droplets as they rise through a stratified fluid, we identify two new timescales which quantitatively describe this slowdown: an entrainment timescale and a retention timescale. These timescales measure dynamics that were not captured in previous timescale discussions, which primarily focused on the timescale to the velocity minimum (U min ). The entrainment timescale is a measure of the time that a droplet spends below its upper-layer terminal velocity and relates to the duration over which the droplet's rise is affected by entrained dense fluid. The retention time is a measure of the time that the droplet is delayed from reaching an upper threshold far from the density transition. These two timescales are interconnected by the magnitude of the slowdown (U u − U min ) relative to the upper-layer terminal velocity (U u ), as well as a constant that reflects the approximately universal form of the recovery of a droplet's velocity from U min to U u . Both timescales are found to depend on the Froude and Reynolds numbers of the system, Fr = U u /(Nd ) and Re = ρ u U u d/ν. We find that both timescales are only significantly large for Fr 1, indicating that trapping dynamics in a relatively sharp stratification arise from a balance between drop inertia and buoyancy. Finally, we present a theoretical formulation for the force enhancement , the ratio between the maximum stratification-induced force and the corresponding drag force on the droplet, based on a simple force balance at the point of the velocity minimum. Using our experimental data, we find that our formulation compares well with recent theoretical and computational work by Zhang

I. INTRODUCTION
There are many examples of droplets, bubbles, and particles interacting with stratified fluids, including atmospheric and marine pollution [1], oil spills [2][3][4][5][6], oil seeps [7,8], marine snow [9][10][11], and transport and motion of microplastics and marine organisms [12][13][14][15]. As a result, understanding when and how stratification affects rising and settling is of significant interest within a variety of fields, particularly in the environment. Further, many industrial processes, although often involving immiscible fluids, rely on the effects of a two-layer stratification to control the motion of drops and particles [e.g., Refs. 16,17], as reviewed by Magnaudet and Mercier [18].
There has been extensive prior work examining particles, drops, and bubbles rising and settling in homogeneous-density fluids, from the vortex shedding and wake dynamics of a sphere [19,20] to basic statistics such as terminal rise velocity of droplets [21][22][23]. However, density stratification adds an additional level of complexity, in the form of additional forces acting on a particle or droplet. When a droplet or particle moves through a fluid, it experiences a variety of forces, including buoyancy, drag, added mass, and history forces. When a droplet rises vertically through a stable stratification, it entrains denser fluid with it, altering the effective buoyancy of the droplet and reducing its upward speed ( Fig. 1) [24]. Furthermore, recent numerical work has suggested the presence of an additional aspect of the stratification-induced force, due to the specific structure of vorticity generated within a stratified fluid [25,26]. This study seeks to quantify and explore the impact of these stratification-induced phenomena on the timescales of delay of a droplet rising to the surface. What is the net effect of stratification on the rising and settling of droplets and particles in a stratified fluid? What are the timescales associated with the induced delay, and how can we parametrize them? And finally, how does a droplet's motion connect to the physics of the problem?
Many prior experimental and numerical studies have been conducted investigating both solid and porous particles settling in stratified fluids, and will be described in detail below. However, to our knowledge, all work on the rising of droplets in stratified fluids has been computational in nature. The current study will thus focus on experimental work, to provide an expansion and validation of previous numerical studies.
The numerical simulations of Bayareh et al. [27] showed that the drag coefficient of a settling spherical drop was enhanced in linearly stratified fluids with drop Froude numbers in the range 4 Fr 16 (where Fr is the ratio of the buoyancy timescale to the inertial timescale). A sharp two-layer stratification was also studied briefly in their Appendix, for both a rigid particle and a drop, and good agreement was found with prior experimental and numerical results. Shaik and Ardekani found that at low Reynolds number, stratification and inertial forces both can increase the drag on a drop, and that this enhancement depends on the dynamic viscosity ratio between the drop and the surrounding fluid [28]. Other numerical work has studied two-droplet [29] and swarm-scale [30] interactions among droplets in linear stratification. These works did not involve a sharp transition between two homogeneous-density fluids, so trapping dynamics were not studied. In a two-layer stratification, the simulations of Blanchette and Shapiro [31] found that the dynamics of oil droplets in stratified fluids may have additional complexities, namely Marangoni forces. These authors found that in a sharp transition in stratification, a drop may either suddenly accelerate through the transition region or be prevented from crossing into the next layer, depending on the relative interfacial tension between the drop and the two layers.
While experimental work on drops is limited, there is significant prior experimental work that has looked at the small-scale dynamics of rigid spheres settling in stratified fluids. The experiments of Srdic-Mitrovic et al. [24] studied the gravitational settling of solid particles through a sharp two-layer stratification and found that stratification drag-that is, an increase in the drag coefficient and an associated deceleration occurring with entry into a stratified layer-was significant in only a narrow range of Reynolds numbers, 1.5 < Re < 15. Otherwise, experiments at Reynolds numbers outside of this range showed no significant change in drag as particles passed through the interface, instead behaving similarly to a particle in a homogeneous fluid. Verso et al. [32] obtained similar experimental results, and further observed that the minimum velocity occurred during a particle's exit from the density transition. Abaid et al. [33] also observed a velocity minimum for a sphere passing through a sharp transition between two homogeneous-density fluids, and in some cases reversal of the sphere's motion. Other experiments [34][35][36] have studied the wake of spheres moving at a constant vertical velocity in stratification, and have found that varying trailing jet structures emerge and contribute to fluid entrainment and mixing. However, these studies were conducted in linear stratification, so the effects of a sharp density transition are unknown in these regimes, and the fluid structures that emerge in a continuous stratification may be suppressed in a two-layer stratification. Studies on rigid spheres in two-layer or linear stratification [37][38][39] have also been conducted at very low Reynolds number. Oceanic applications, such as rising bubbles or oil droplets from spills in the ocean may have Reynolds numbers ranging from intermediate to high, depending on estimates of drop diameter and drop velocity [40,41], indicating the necessity of studies to be performed beyond the low Reynolds number regime and over a range of density transitions.
Experimental and theoretical studies of settling porous spheres have also predicted increased drag or prolonged retention times at sharp density transitions [11,[42][43][44], due to either diffusion of lighter fluid into the settling porous particle, or to entrainment of lighter fluid from above. However, most of these studies were limited to the Stokes regime (Re 1), and the observed retention was primarily driven by diffusive processes.
Finally, recent theoretical and computational work by Doostmohammadi et al. [25] and Zhang et al. [26] addressed the "stratification drag" force that is commonly used as a catch-all for the contribution of stratification to changes in an object's motion, including increased residence time or significant slowdown [24,32]. Doostmohammadi et al. identified an additional mechanism, the generation of vorticity by baroclinic torque (nonalignment of density and pressure gradients), which leads to increased shear at a particle's surface. Zhang et al. built upon this work by methodically decomposing contributions of stratification to enhanced force into two components: (1) modified buoyancy forces, due to the relative buoyancy of entrained fluid dragged behind a sphere, and (2) modification of the local vorticity field due to baroclinic torque, inducing an increased shear stress at the surface of the sphere, which we will refer to here as the baroclinic vorticity force. They presented a rigorous scaling of these two forces in different Prandtl, Reynolds, and Froude number regimes. While previous literature has focused heavily on the contributions of this first force due to the buoyancy effects of entrained fluid, the baroclinic vorticity contribution to stratification forces had not previously been identified, and experimental measurements corroborating this approach are lacking. The present work aims to quantify and explain the retention of single oil droplets at the transition between two homogeneous-density fluids by a methodical study of two new timescales describing a droplet's delay, as well as the forces contributing to those delays. Using laboratory experiments, we examine motion and retention for a range of drop sizes, drop densities, and ambient stratification profiles. In Sec. II, we discuss the nondimensional parameters relevant to this problem. In Sec. III, we will describe the experimental setup and measurements taken. We will discuss our results in Sec. IV, beginning with analysis of the drop's position and velocity, the timescale associated with the drop's initial deceleration, and then introduce two timescales relating to fluid entrainment and droplet retention. We find that these entrainment and retention timescales are related, and that they are dependent upon the Froude and Reynolds numbers of the system, but largely uncorrelated with the timescale associated with deceleration. Further, significant fluid entrainment or drop retention only occurs for Fr 1. In Sec. V, we develop a theoretical formulation for the force enhancement, , induced by stratification. The scaling of compares very favorably to the work of Zhang et al. [26] and provides the first experimental evidence supporting their approach of decomposing the stratification forces into modified buoyancy and baroclinic torque-induced shear. Finally, we will close with a discussion of the implications of this work and future directions in Sec. VI.

II. NONDIMENSIONAL PARAMETERS
Given the importance of density stratification in both environmental and industrial processes, we covered a range of the parameter space relevant to this problem, particularly in the intermediate Reynolds number regime. Table I lists the parameter definitions and ranges covered in this study, which spanned 179 different droplet experiments. Nondimensional parameters with the subscript f , such as Re f , encompass two separate nondimensional numbers for a single drop's behavior in the upper and lower layers of ambient fluid, Re u and Re l . The subscript f in the given definition is thus replaced by the corresponding upper (u) or lower (l) layer quantity. The subscript d represents the corresponding nondimensional number or parameter for the drop fluid properties. In these definitions, ρ f represents fluid density, U f represents the terminal drop velocity in a given homogeneous-density region, d is the drop diameter, μ is the dynamic viscosity of the fluid, ν is the kinematic viscosity of the fluid, and h is the thickness of the transition region, computed as the height encompassing 95% of the density variation between the upper and lower layers. frequency, N, is defined as N = √ −g/ρ u (∂ρ/∂z), where ∂ρ/∂z is computed as the least-squares slope of a 0.6 cm-wide region in the transition region of the density profile at rest, centered at z = 0.
Following the definitions given in Table I, the Reynolds number (Re f or Re d ) represents the ratio of inertial to viscous forces. The Archimedes number (Ar f ) is the ratio of buoyant forces to viscous forces. The Froude number (Fr) can be thought of in several ways: (1) as the ratio of flow inertia to external gravitational forces; (2) as the ratio of the buoyancy timescale (1/N) to the drop motion timescale (d/U u ); or (3) as a ratio of the speeds at which various information about the flow is propagating, i.e., the ratio of droplet speed to an internal wave speed. In experiments, these parameters were varied by changing the drop diameter, the drop density, and the transition region thickness (which in turn changes N). In particular, by varying the transition region thickness, we are able to vary the Froude number independently of the Reynolds number. The Prandtl number, the ratio of momentum diffusivity (i.e., kinematic viscosity ν) to salt diffusivity (κ), remained fixed at approximately 600 for all experiments.

A. Experimental setup
Experiments were conducted in a 61 cm tall acrylic tank with a width and depth of 30.5 cm by 30.5 cm. A schematic of this tank is shown in Fig. 2(a). Sodium chloride (Morton Canning & Pickling Salt) was used as the stratifying medium. Fluids for the two layers were prepared in two 35-gallon tanks with recirculating pumps, which were filled with reverse osmosis water. Salt was added to one tank and dissolved. Both tanks were left to circulate at room temperature to eliminate convection in the filled experimental tank.
The experimental tank was filled using two methods: (1) a two-layer filling method that yielded error function-type density profiles, and (2) a computer-controlled method, yielding linear density profiles in the transition region. For the first method, the tank was filled first with salt water (ρ l = 1.106 to 1.117 g/cm 3 ) and then with fresh water (ρ u = 0.9972 to 0.9981 g/cm 3 ), with a sponge float acting as a diffuser to reduce mixing between the two layers. To obtain a thin transition region (3-4 cm), the tank was allowed to sit and diffuse for 2-3 h, until the optical distortion caused by the difference in refractive indices between the two layers had reduced. To obtain a thicker transition region (7-8 cm), the tank was allowed to diffuse another 18 h. This two-layer filling method yielded error-function shaped density profiles, as seen in Figs. 3(a) and 3(b). The second filling method allowed more precise control of layer thickness. For this method, two computer-controlled peristaltic pumps (New Era Pump Systems NE-9000) feeding from a fresh water bucket and a salt water bucket were linearly ramped up and down to generate a linear stratification between the upper and lower layers. These could be programed to yield a range of transition layer thicknesses. Examples of such density profiles are shown in Figs. 3(c) and 3(d). The gray region in these density profiles represents the layer thickness h, where 95% of the density variation between the upper and lower layers occurs.
Oil droplets were composed of a mixture of 10 cSt silicone oil (Clearco Products Co.) and halocarbon oil (Sigma Life Science Halocarbon oil 27) to study a range of drop densities ρ d , from 0.9375 to 0.9927 g/cm 3 . Oil fluorescent tracer (Risk Reaction DFSB-K175 UV Orange) was also added for contrast. Drops were released individually by dispensing a small amount of oil using a syringe pump or handheld syringe, which fed into a 19 gauge needle inserted through a flange in the base of the tank. A waiting time of at least 15 min between drop releases was chosen to ensure the tank was quiescent for each experiment.

B. Experimental measurements
Physical characteristics of the ambient fluid and oil droplets were measured prior to experiments. An Anton Paar Lovis 2000 ME microviscometer and DMA 4100 M densitometer were used to measure the viscosities and densities of the fresh water, salt water, and droplet fluid. The densitometer also provided direct measurements of the fluid's temperature.
During experiments, images of the injected drops were captured at 120 to 125 fps using a highspeed camera (Photron FASTCAM SA3 at 1 MP, Point Grey Grasshopper3 at 5 MP) aligned with the plane of drop motion. A panel of light emitting diodes (LEDs) was placed behind the tank, along with a diffusive screen of vellum paper between the tank and lights, to increase the contrast between drops and the background. Before each drop was released, an image was taken of the field of view, including a calibration ruler placed in line with the needle and plane of droplet motion. Because of the tank's density stratification, the refractive index encountered by a light ray changes as light passes through the tank. The images captured by the camera thus have refractive distortion. This distortion was corrected by calibrating the drop position relative to the refracted ruler image, as shown in Fig. 4. In some cases, the drops exhibited slight out-of-plane motion, which we estimate to contribute 1% or less error in measured vertical position, based on a camera distance of ∼1 m and out-of-plane motion on the order of 1 cm.
Following distortion correction and subtraction of a mean background image, drop position over time was then tracked using the Trackpy software package [45], which uses center-of-mass detection to determine droplet position. Example tracked paths and velocities for five drops are shown in Figs. 5-7. Instantaneous velocities were obtained following the methods of Srdic-Mitrovic et al. [24], in which a least-squares line was fit to a window of seven points of vertical position (∼0.06 s of data) and the best-fit slope was assigned as the velocity of the center point in the window. Upper and lower layer terminal velocities, U u and U l , were computed as the least-squares slope of the drop trajectory in regions with constant speed in the upper and lower layers, respectively.
Drop diameters were measured manually from an image in the lower layer, and calibrated from pixels to centimeters. For 145 cases, manually measured and calibrated diameters were verified against images of drops taken using a telecentric lens; telecentric images yielded diameters that varied on average by 0.16 mm from the other method (an average relative difference of less than 5%), giving an estimate of the error in manual measurement and calibration.
In addition to correcting for refractive distortion, the quantified distortion of the calibration ruler was also used to determine density profiles using synthetic Schlieren [46,47]. Optical focusing and spreading of ticks on the ruler within the transition region were compared with the even spacing in the upper and lower layers, and the apparent displacement z of these ticks was then converted to a density gradient following assumptions of linearity, two-dimensionality, and small incident ray angles, as described in the above references. The full density profile (such as those shown in Fig. 3) was obtained by integrating the gradient from the known upper layer density using the following equation: where z is the measured apparent displacement field (the difference between the curve shown in Fig. 4(b) and vertical height), and β 1.88 s 2 cm −1 following Ref. [46]. Because precise measurement of the exact distance between the ruler and the tank side wall, L ruler , was difficult, and precision was difficult to maintain from experiment to experiment over 179 diffferent runs, this length was adjusted manually by ±0.8 cm to yield a profile whose constant upper and lower layer densities matched those measured with the hand-held densitometer. Finally, for five of the above experimental cases, shadowgraph experiments were performed to visualize the wakes of droplets. A schematic of this setup is shown in Fig. 2(b). Polyester drafting film (West Design Polydraw) was placed on one side of the tank, and a camera was placed facing the drafting paper so that drops and their wakes could be visualized via the focusing and defocusing of incoming light rays. A collimated light source (Thorlabs M450LP1 450 nm LED, in conjunction with an Edmund Optics 200 mm diameter, 800 mm focal length PCX condenser lens) was placed on the opposite side of the tank. The LED light source and collimating optics were set to angle downwards at about 30 degrees from the horizontal to avoid total internal reflection within the transition region, which would have occluded ∼1 cm of the droplet's path and also caused oversaturation in images. The images presented here thus show the projected fluid structures viewed at this angle, rather than a perfectly perpendicular view of the x-y plane. A tracking camera was placed perpendicularly to the shadowgraph camera, opposite the LED backlighting panel, which was set to emit green light. Each camera was equipped with a bandpass filter (Thorlabs FELH0550 long-pass filter with a cut-on wavelength of 550 nm, and Thorlabs FES0500 short-pass filter with a cut-off wavelength of 500 nm) so that illumination from the backlighting panel and the 450 nm LED could be separated. Both tracking and shadowgraph images were synchronized using a function generator (Siglent SDG1025) that triggered the two cameras externally. Before a set of shadowgraph experiments, a calibration image was taken for each camera with a clear acrylic ruler in the field of view; the ruler was then moved to the edge of the tank and the tank allowed to settle for approximately 30 min before releasing droplets. Shadowgraph images were post-processed by correcting for optical distortion and then subtracting a background image.

IV. RESULTS
We begin by discussing basic properties of the droplet motion. We will first demonstrate how droplet position and velocity vary with experimental conditions. We then present parametrizations of the terminal velocity of the droplets in the upper and lower layers, which will be useful for computing stratification-related force enhancement in Sec. V. In the next subsection, we analyze the timescales over which deceleration, fluid entrainment, and significant droplet retention occur, discuss their dependence on the nondimensional parameters of the system, and delineate when drops are significantly retained at the transition region. We then briefly connect these timescales to flow visualizations of the droplet's wake.

A. Drop paths and velocities
Shown in Figs. 5 and 6 are a sequence of shadowgraph images for (A) a larger (d = 0.40 cm) and (B) a smaller (d = 0.28 cm) droplet, both composed of the densest oil mixture (ρ d = 0.9927 g/cm 3 ) and rising through a 4.6 cm transition region. Relative shading indicates variation in the second derivative of density [48]. As the drops exit the transition region, an internal wave field is generated. Also shown in each figure are the corresponding vertical drop position, z, and velocity, u, as a function of time, t, for each set of shadowgraph images. The droplets slow as they pass through the transition region (snapshots (a)-(c) in both sets of shadowgraph images), and reach a velocity minimum [snapshots (e)-(f)] just above the transition region. The drop then eventually regains speed [snapshots (g)-(i)], asymptoting to its upper-layer terminal velocity U u , indicated as the dashed gray line in the velocity plots, by about snapshot (j). A complex wake structure can also be observed in the shadowgraphs, which will be discussed in detail in Sec. IV D.
To demonstrate how varying different experimental parameters affects drop motion, Fig. 7 shows sample drop paths and velocities over time and drop velocity as a function of vertical position for four example experimental droplet cases. The first case, shown in Figs. 7(a), 7(e), and 7(i), is the droplet shown in Fig. 5. We use this small, dense droplet passing through a 4.6 cm transition region as a reference case for comparison with: (b), (f), (j) a similarly sized, lighter droplet in similar ambient stratification; (c), (g), (k) a larger, dense droplet in similar stratification; and (d), (h), (l) a small, dense droplet passing through a thicker transition region of 7.1 cm. It can be seen that lighter droplets have significantly higher terminal velocities in the upper layer, and that small, dense droplets take significantly longer to traverse the field of view of the camera. The denser drops (e), (g), (h) and (i), (k), (l) also remain at a speed lower than their upper layer terminal velocity U u for an extended period of time, indicating that entrained ambient fluid plays a role in delaying the drop's upward motion. Analysis of these delays will be presented in Sec. IV C.  [24]. Empirical equations from Ref. [49] are shown for the terminal velocity of solid spheres as the gray dot-dashed line, and for the terminal velocity of liquid drops in air as the gray dotted line.

B. Terminal velocity
Droplets reach a constant, terminal velocity in both the upper and lower layers, and a parametrization of this terminal behavior will be useful when considering the force balance in stratification in Sec. V. This terminal velocity is governed by buoyancy, viscosity, and inertia. The drop's Reynolds number in the upper and lower layer is plotted versus its Archimedes number in each layer in Figs. 8(a) and 8(b). As noted in Table I, the definitions of these two nondimensional numbers is as follows: where the Archimedes number is the ratio of buoyant to viscous forces and the Reynolds number is the ratio of inertial to viscous forces. As before, the subscript f represents either the upper (u) or lower (l) layer property, and U f represents the terminal drop speed in that corresponding layer. Each drop thus yields two data points of Ar f and Re f . Measurements of Ar l and Re l are represented as black circles, while values of Ar u and Re u are shown as blue squares.
In Fig. 8, as well as in later figures, we assessed deviation between data and empirical fits using coefficient of determination (R 2 ) and the mean absolute percent error (MAPE). A definition of this error metric, as well as statistics on all fits presented in this paper, are given in Appendix A. A power-law fit was found that describes the relationship between the Archimedes and Reynolds numbers from experimental measurements: where α = 0.61 and β = 0.57, with R 2 = 0.99. The value of this power-law coefficient β relates to the scaling of the drag force with Reynolds number. As we will see in Sec. V, the drag force scales with U 1/β . From Stokes' law, for Re 1, the drag force is proportional to U 1 , so β = 1 [the leftmost sloped line on Fig. 8(a)]. For fully turbulent flow, the drag force follows the quadratic drag law, and is proportional to U 2 . Thus, for high Re, β = 1/2, shown as the rightmost sloped line on Fig. 8(a). Our data fall in the intermediate Reynolds number regime, so a scaling result of β = 0.57, in between 0.5 and 1, is expected.
In Fig. 8, we have also plotted the Reynolds and Archimedes number data for solid particles (of different densities and diameters, classified by P1-P4) presented in Verso et al. [32] and Srdic-Mitrovic et al. [24], all of which are summarized in Tables 4-8 of Ref. [32]. These data agree well with our interpretation that particles or droplets at lower Reynolds number have a value of β closer to 1, shown as the slope in the leftmost portion of Fig. 8(a) which lies primarily near those data points.
This power-law fit compares reasonably well with the relationship for solid spheres in the range 435 < Ar 1.16 × 10 7 and 12.2 < Re 6.35 × 10 3 ([49], Table 5.3): where W = log 10 ((4/3)Ar), shown as the dot-dashed line in Figs. 8(a) and 8(b). A similar equation for liquid drops in air from [50] and [49] (Eq. 7-1) is also plotted, which in the regime studied here behaves almost identically to that for solid spheres. In summary, for drops in the parameter space studied here, if the viscosity and density of the ambient fluid, as well as density and diameter of the drop are known, then the drop's terminal speed can be predicted with reasonable accuracy using the relation given in Eq. (4). This relationship, whose exponent β is determined by the drop's Reynolds number and describes the drag force scaling with velocity, will be used later in Sec. V for our theoretical formulation of force enhancement in stratification.

C. Timescales characterizing droplet slowdown and delay
As shown in Figs. 5-7, some droplets experienced a significant slowdown as they passed through the transition region. First, we will look at a timescale that reflects the deceleration of the droplet as it enters the transition region. In their studies of solid particles, Srdic-Mitrovic et al. [24] defined a timescale t min as the time between a particle's entrance into the stratified layer, and when it achieves its minimum velocity, as shown in Fig. 9. This deceleration happens quite rapidly, and in general t min is less than 2 s (true for > 86% of experimental cases). Following Srdic-Mitrovic et al., we 124803-13  Fig. 10.
In Fig. 10(b), we observe the same Re −1.7 l scaling as has been observed in solid particles [24,32], in which the timescale of deceleration was primarily determined by the properties of the homogeneous layer the particle was exiting (in our case, the lower layer). Interestingly, Verso et al. [32] observed two power laws in their data: −1.7 and −3.4, and attributed this to other dimensionless quantities. It is unclear what physical process led to the −3.4 scaling observed in Verso et al., and we do not observe that power law in the data studied here.
The deceleration timescale t min collapses when scaled with a viscous diffusion timescale, d 2 /ν, indicating the importance of viscous processes in this deceleration. As shown in Figs. 10(a) and 10(c), this timescale is uncorrelated with the Froude number, indicating that the characteristics of the ambient stratification, i.e., N, do not control this deceleration process. Further, the power law scaling of Re l stays more or less constant (≈ −1.7) when fitting alone or with Fr l [ Fig. 10(b) vs. Fig. 10(c)]. Thus, t min provides an incomplete picture of overall retention effects.
To get a complete understanding of droplet trapping due to stratification, it is therefore useful to look at metrics of drop retention that measure both the duration of time that the drop's motion is affected by stratification, as well as the duration of time that the drop is physically retained by the transition layer. These net retention effects are particularly important for the environmental, oceanic, and industrial processes discussed earlier.
We will consider the first metric to be the entrainment time, t e . This timescale is demonstrated in Figs Fig. 11(a)] and a line representing the upper layer terminal speed [the gray dashed line in Fig. 11(a)] is less than 5% of the drop diameter. This timescale is very similar to the delayed settling time (DST) used by Prairie et al. [11].
The second metric is a retention time, t r , the time duration that the droplet is retained and slowed from rising to the surface. This time is shown in Fig. 11  To understand the physical relationship between t e and t r , a physical and geometrical argument can be constructed using the scaling of the retention distance d r [shown in Fig. 11(c)], the extra distance the droplet would have traveled had it not been retained. By definition, the area of the blue shaded region in Fig. 11(b) between U u and u(t ) is equal to the retention distance d r : where t 0 is the time marked by the first blue triangle in Fig. 11(b), and U = U u − U min . Dividing both sides by t e U , we obtain By definition in Fig. 11(c), This factor, d r /(t e U ), is a key parameter with a physical interpretation as the ratio of the actual retention distance to the retention distance had the droplet traveled at a speed U min for the duration of t e . In our data, this parameter appears to be approximately constant, with a mean value of 0.47 [ Fig. 12(a)]. We can examine the area of the blue shaded region in Fig. 11(b) that was represented by d r , and imagine drawing a rectangle of length t e and height U [ Fig. 12(b)]. The fact that this parameter is equal to approximately 1/2 across all experiments, and that U min is reached quickly (i.e., t min is small compared to t e ), implies there is an approximately universal behavior in the recovery from to U min to U u .
We can now return to Eq. (8), and plug in our mean value for d r /(t e U ), c = 0.47: The retention time is, therefore, the entrainment time multiplied by the constant, c = 0.47, which represents the universal form of the recovery of the drop's velocity, and by a factor representing the relative magnitude of the drop's slowdown, U/U u . Note that this does not imply a linear relationship between t e and t r , as U/U u varies with the experimental parameters as well. This relationship is shown in Fig. 13. The physical retention of a drop at the transition region is thus 124803-15 determined by the duration of time that denser fluid is appreciably entrained, t e , and the relative drop slowdown, U/U u , which may be due to the amount of dense fluid entrained as well as any other contributions due to the distortion of isopycnals.
To understand how environmental conditions affect these timescales, the retention and entrainment times were then compared against nondimensional parameters governing the drop's rise. Figures 14 and 15 show the nondimensional entrainment time, τ e = t e N, and nondimensional retention time, τ r = t r N, as a function of the Froude number and Reynolds number in the upper layer. These timescales are more strongly correlated with the Froude number than the Reynolds number [see panels (a) and (b) in Figs. 14 and 15], indicating the importance of stratification in drop retention. Since the timescales show significant correlation to both Re u and Fr, we conclude that a scaling that includes both, shown in panel (c) of Figs. 14 and 15, is most appropriate. Other nondimensional numbers, including Re l , h/d, and a lower-layer Froude number, were also compared but did not yield significant collapse of the data and so are not presented here.
Nondimensional entrainment and retention times, with their associated Reynolds scaling from Figs. 14(c) and 15(c), are plotted versus Fr in Fig. 16. Both of these temporal metrics have a power law relationship with Froude number, with larger nondimensional retention and entrainment times occurring for small Fr. A delineation can be drawn at Fr ≈1 (shown as the dotted lines in Fig. 16), and using the definition of the Froude number as the ratio of inertial to buoyancy forces, we can split droplet behavior into two regimes: (1) For Fr 1, the drop is rising on a timescale (d/U u ) less than the buoyancy timescale (1/N), and the drop's motion is not significantly impacted by interactions with the ambient stratification. Retention times are close to zero, as inertial forces dominate over buoyant forces.
(2) For Fr < 1, significant drop retention is observed. As the buoyancy timescale (1/N) becomes less than the inertial timescale (d/U u ), the drop is more significantly retained in the transition region. This corresponds to buoyancy forces being strong enough to counter the drop's inertial forces. Conceptually, one may also consider that for Fr < 1, the drop's velocity is less than a characteristic internal wave velocity, Nd, and the drop does not have the kinetic energy required to "punch through" the transition region.
These results demonstrate that in the regime covered by our experiments, the dynamics of droplet trapping are primarily governed by a balance between buoyancy and inertia, and that the upper-layer Froude number is an important parameter to consider when predicting trapping or slowdown in a two-layer stratification.
Finally, we have compared the timescale of deceleration used in previous literature, τ min , with the two timescales developed in this work, τ e and τ r . Figures 17(a) and 17(b) show little correlation between τ min and our approach, while as demonstrated in the preceding discussion, τ e and τ r are strongly related [ Fig. 17(c)]. In effect, the processes of fluid entrainment and drop retention are not a function of how long it takes for a droplet to slow down. τ e and τ r are functions of the Reynolds and Froude numbers of the upper layer (i.e., the layer the droplet is entering); τ min is a function of Reynolds number of the lower layer (the layer the droplet is exiting), and is independent of the Froude number (and therefore properties of the ambient stratification). Thus, these two approaches to studying the timescales of drop motion are capturing significantly different physics. The deceleration time τ min captures viscous effects, including how long it takes for the droplet to respond to entering the transition region, and generally happens quite rapidly. The entrainment and retention times capture the other half of droplet behavior that had largely remained unquantified: the overall effect of entrained fluid and interaction with the ambient stratification on a droplet's rise, occurring over a longer timescale. Our overarching goal in this work was to quantify net trapping effects. In environmental or engineering applications where the overall rise time is of interest, the new timescales and their parametrizations developed here capture the physics of stratification-related retention well.

D. Fluid entrainment
Above, we have given a conceptual model of how stratification ultimately affects the residence time of droplets at a density transition. Here, we briefly discuss visualizations of the wakes of the droplets first shown in Figs. 5 and 6. Figures 5 and 6 showed a sequence of shadowgraph images for (A) a larger and (B) a smaller droplet, both composed of the densest oil mixture (ρ d = 0.9927 g/cm 3 ) and rising through a 4.6 cm transition region. The most obvious difference between these two cases is the asymmetry of the droplet's wake. The larger droplet shown in case A, with Re l = 370, appears to be shedding vortices in a zigzag pattern, while the smaller droplet with Re l = 196 has a highly symmetric wake structure. This aligns with the delineation of the effect of Reynolds number on the wakes of rising and falling spheres in a homogeneous fluid presented by Horowitz and Williamson [19], who found that wake structures transition from vertical to oblique at Re = 210, and from oblique to zigzag at Re = 260. The results are also qualitatively very similar to those for spheres in a linear stratification [34]; in our experiments, the finite transition between two homogeneous-density fluids appears to constrain the wake structure to the transition region.
While the vortex structures shed by the drop as it leaves the homogeneous lower layer vary significantly between these two experimental cases, the two cases have very similar Froude numbers (0.55 and 0.53), and indeed very similar entrainment and retention times (for example, t e = 4.95 and 4.18 s; τ e = 29.5 and 24.9, respectively). This variation in wake structure when entering and initially exiting the transition region thus appears to have little effect on fluid entrainment and drop retention. In Appendix B, we present some preliminary measurements of the local wake diameter from these shadowgraph images, which suggest that there is correlation between the decay of the trailing denser fluid and the entrainment timescale.

V. FORCE ENHANCEMENT BY STRATIFICATION
In this section we will discuss the forces acting on the droplets and how they scale with the parameters of the problem. In particular, we wish to estimate the stratification force on the droplet, F s , which is the excess force experienced by the droplet in a stratified fluid, relative to the force on a droplet in a homogeneous fluid. The timescales developed above measure the duration of delay, and this stratification force is what causes this delay. Here, we will present a theoretical formulation for this force enhancement due to stratification and compare it to the parameter space and force balances presented by Zhang et al. [26] for numerical studies of rigid spheres.
The forces acting on particles or droplets passing through a stratified fluid are typically decomposed into several different terms [24,32,49]: where F b is the buoyancy force, F d is the drag force, F a is the added mass force, F h is the Basset history force, and F s is any additional force attributable to interactions with the ambient stratification, including the modified buoyancy force, the baroclinic vorticity force, and potentially Marangoni forces (the latter discussed in more detail in Sec. VI and Appendix C). Here, we assume F b and F d correspond to the buoyancy and drag forces, respectively, which would be experienced by a droplet moving through the undisturbed ambient fluid. (Therefore, the density corresponds to the value for the undisturbed background fluid at the current location of the droplet.) Other experimental work in a similar parameter regime has shown the added mass and history forces are typically much weaker than other forces in the system, and can be neglected [32]. This behavior was further observed in numerical studies in linear stratification, where the history force was negligible for all time, and the added mass force was generally small [25]. In Appendix C, we provide a thorough estimate of the order of magnitude of each of these forces to support our assumption that they are negligible throughout most of the droplet's motion. First, we will examine the functional form of the buoyancy and drag forces in a homogeneous fluid, F b and F d , far away from the transition region. The buoyancy force can be computed as a function of the local fluid density ρ f , which varies with depth z [24]: where F 0 = πμ 2 6ρ f . We can relate this buoyancy force to the local terminal velocity, U f (z), using the relationship measured in Sec. IV B, 124803-19 and the definition of Reynolds number, Re f = ρ f U f d/μ, and rearrange to obtain an expression for the Archimedes number in terms of the local droplet velocity, where U 0 = αμ ρ f d is the characteristic velocity if Ar f = 1, and the empirical measurements discussed in Sec. IV B gave α = 0.61 and β = 0.57.
In a homogeneous fluid, a droplet will reach its terminal velocity when F d = −F b . Using this relation and Eqs. (11) and (13) we can obtain a drag law for the drag force in a homogeneous fluid as As mentioned in Sec. IV B, we expect this exponent, 1/β, to range between 1 at low Reynolds number (where drag scales as U ) and 2 for high Reynolds number, turbulent flow, where drag scales as U 2 . As expected, in the intermediate Reynolds number regime studied here, this exponent falls between these two limiting values, with 1/β = 1.75.
When considering the full force balance, including the stratification force F s , we will focus on the point of minimal velocity (t = t min , u = U min ), which previous studies have shown is when the stratification force is at its maximum [25,32]. This approach is similar to that of Doostmohammadi et al. [25], who examined the force balance on a particle at its peak velocity in a linear stratification. In addition to being more tractable at this single point in time, the instantaneous acceleration cannot be accurately measured from the experimental data due to error propagation in the second derivative, so the time dependence of the inertial, added mass, and history forces is difficult to compute. Regardless, at this point of velocity minimum, Du Dt = 0, the added mass is defined using the acceleration (F a = − 1 2 C A ρ f V p Du Dt ) and so is also zero, we assume F h is negligible [25,32, Appendix C], and the droplet is usually entirely in the upper fluid layer (see, e.g., Figs. 5-7). The force balance from Eq. (10) can thus be assumed to simplify to where F b and F d are the local buoyancy and drag forces when the droplet has just entered the upper layer, and is moving at its minimum velocity U min . This local buoyancy force is computed for a droplet that has fully entered the upper layer using Eqs. (11) and (13), Using Eq. (14), the local drag force scales with the local velocity to the power of (1/β ), Thus, we can write a force enhancement ratio, , which compares the maximum stratification force experienced by the droplet to the drag force, and allows for direct comparison to previous results on solid particles, This quantity can be understood in the following manner: When the droplet speed is at a minimum, i.e., a velocity deficit relative to the upper layer terminal velocity, the absolute value of the drag force F d is less than that of the buoyancy force F b . However, at the velocity minimum, the droplet is not accelerating, so an additional force must be provided by the the stratification force, given in Eq. (18).
For our data, a power-law best fit was found as a function of Fr and Re u , with = 19 Fr −0.89 Re −0.56 u (Fig. 18). These exponents agree remarkably well with the dependence of the drag force due to vorticity induced by baroclinic torque, F ρω , proposed by Zhang et al. [26]. In their Regimes 2 and 3 (Pr 1/3 Fr, in which all of our experimental measurements fall), they found that Some of our data fall toward the border of Zhang et al.'s Regime 1; for large Re and Pr, Regime 1 is delineated as Fr Pr −1/6 , or Fr 0.34 in our case. Our smallest Froude number is 0.38 (Fr is indicated by the colorbar in Fig. 18). It appears that even at the lower limit of this regime, there is fairly good agreement.
In summary, we see that the governing velocity scales relating to the droplet slowdown (U u and U min ), which govern the behavior of t e and t r , are also the primary variables in a measurement of the force enhancement due to stratification, (U u /U min ) 1/β − 1. These experimental and analytical results regarding the force enhancement show that the baroclinic vorticity contribution is likely important in our work, and provide the first experimental data supporting the approach of Zhang et al. [26]. Both of the relevant forces that contribute to "stratification drag" discussed by Zhang et al.-the force due to relative buoyancy of entrained fluid (F ρρ ), and the force due to modification of the local vorticity field due to baroclinic torque (F ρω )-are predicted to scale with Fr −1 in the parameter space covered by our experiments; moreover, F ρρ does not vary with Reynolds number, but F ρω does. Zhang et al. showed that F ρω should dominate in the Reynolds number regime studied here, which results in our observed scaling of with both Fr and Re u . We do expect F ρρ to play a role in this regime based on our timescale results and flow visualization of the droplet wake; however, its contribution is likely less than that of F ρω . It should be noted that the simulations of Zhang et al. were conducted in a linear stratification, so their definitions of Reynolds and Froude numbers vary from the definitions used in this study. Regardless, these experimental results demonstrate that the baroclinic vorticity force F ρω is very likely significant in the regime studied here, and that experimental measurements  of the velocity and vorticity fields surrounding drops and particles in stratification are warranted to confirm this.

VI. DISCUSSION
We have studied the retention and entrainment dynamics of droplets in the regime 0.38 < Fr < 4.2, 59 < Re l < 1060, 5.4 < Re u < 580, with relative drop densities ρ u ranging from 0.0045 to 0.061, and h/d ranging from 4.5 to 52. Counter to the results of Srdic-Mitrovic et al. [24], who found that experiments at higher Re u (>15) showed no significant change in drag as a sphere settled through a density gradient, we observed significant drop retention and slowdown, with τ e > 10 observed for Re u up to 300, and over a wide range of Re l . As noted in previous work [32,33], the parameters of the layer that a spheroid is entering appear to be the most critical for observing "levitation" or significant delay. This holds true in our study, in which Froude number and Reynolds numbers based on the upper layer terminal velocity are the governing parameters.
Significant retention was primarily observed for Fr 1 and and ρ u 0.035. Within existing literature that has examined deceleration and retention of particles and solid spheres in linear or sharp stratification, the parameter regime in which a significant delay occurs has not been well quantified. Srdic-Mitrovic et al. [24] only observed significant slowdown for Reynolds numbers (based on the velocity on entry into the stratified layer) between 1.5 and 15. The numerical simulations of Torres et al. [51] at intermediate Reynolds number found drag to strongly increase with Fr −1 for Fr < 20; however, it was found that this increase in drag was due to a rear buoyant jet that persists in a continuous stratification, and may not be applicable to the relatively sharp two-layer stratification studied here. Other work, including Yick et al. [37], found that enhanced drag scales with Ri 0.51 in the very small Reynolds number regime (where the Richardson number Ri = Re/Fr 2 , implying the same Fr −1 scaling). We also explored only a certain Froude and Reynolds number regime in this study; it remains to be seen whether this scaling applies to drops that are in the Stokes regime, or for extremely high Reynolds number spheroids, such as rising bubbles.
Finally, we have not systematically varied the Marangoni forces in this study. A brief orderof-magnitude estimate of this force is given in Appendix C. More importantly, because we are treating the stratification-related force F s as the remainder of our force balance, any Marangoni force effects due to the ambient stratification are captured in our measurement of the stratification drag enhancement . Our results regarding the scaling of force enhancement with Reynolds and Froude number are comparable to the literature from solid particles, suggesting that this is a valid approach.
The dynamics explored in this study are applicable to a range of environmental scenarios, including oil spills and natural oil seeps. Although our results are for liquid droplets in an ambient stratification, we expect some of these findings to hold for solid particles as well, and may have applications in sediment suspension in benthic boundary layers [52] and dispersal of pollutants in the atmosphere [1].

VII. CONCLUSIONS
In this study, we characterized the dynamics governing retention of a single droplet at a transition in density between two homogeneous fluids, focusing on timescales that describe the net retention effects, and the excess force causing this retention. We examined fluid flow and droplet retention for a range of drop sizes, drop densities, and ambient stratification profiles, allowing us to characterize drop behavior for a range of Reynolds and Froude numbers. We found that far from the density transition, within the homogeneous fluid layers, the droplets followed a balance between buoyancy, inertial, and viscous forces, and that Re f ∼ Ar β f , where β = 0.57. We first studied the timescale to velocity minimum, t min , which has been the focus of previous work [24,32] but only captures the initial rapid deceleration process. We developed two metrics  measuring the timescale of drop delay at a density transition which focus on the velocity deficit and overall recovery to terminal velocity. The first metric, the entrainment time t e , measures the amount of time that denser fluid is appreciably entrained, reducing the drop's speed. The second, the retention time t r , measures the degree to which the droplet's rise is delayed. The retention time is related to the entrainment time by a simple linear relation involving the magnitude of the slowdown the droplet experiences and a constant (approximately 1/2) that describes the universal form of the recovery of droplet velocity, t r = 0.47t e ( U/U u ), where U = U u − U min . The timescales and simple scaling arguments from our data lend themselves particularly well to experimental quantification of the trapping dynamics, as they are readily measurable from kinematic data and do not require taking derivatives or quantifying various forces at all points in time.
In the regime covered by our experiments (0.38 < Fr < 4.2, 59 < Re l < 1060, 5.4 < Re u < 580, 0.0045 < ρ u < 0.61), nondimensional entrainment and retention times were found to depend on the Froude number and Reynolds number. Significant retention with either timescale only occurred for Fr 1, suggesting that retention is primarily a function of the ratio of the buoyancy timescale (1/N) to the inertial timescale (d/U u ), and that trapping dynamics are dominated by the effects of stratification. We also found that stratification-related retention appears to be independent of the type of large-scale wake (zigzagging or vertical) the droplet has when first entering the transition region.
Finally, we examined the forces that drive the entrainment and retention time metrics developed here. Based on a force balance at the point of the droplet's minimum velocity, we derived a force enhancement ratio = F s,max F d = (U u /U min ) 1/β − 1, and found a power-law fit to our experimental data that compares very favorably with the behavior of the baroclinic vorticity force proposed by Zhang et al This suggests that forces generated by the baroclinic torque, likely in addition to the more traditionally studied entrainment-related modified buoyancy forces, are important in this problem and in the regimes studied here.

ACKNOWLEDGMENTS
We are grateful to the reviewers for their thoughtful and constructive comments on this manuscript. S.K. acknowledges support from the Hellman Faculty Fellows Fund. The authors thank Thomas Bellotti and Joshua Roe for contributing to numerous discussions on this project. The Photron cameras used in this study were generously lent by Jian-Qiao Sun. The authors also acknowledge the Joint Applied Math and Marine Sciences Fluids Lab at the University of North Carolina at Chapel Hill, where S.K. was a postdoctoral researcher, including Roberto Camassa, Richard McLaughlin, Holly Arrowood, Lauren Colberg, Elaine Monbureau, Sarah Spivey, and Arthur Wood for contributing to preliminary stages of this project.

APPENDIX A: STATISTICS
In this section we briefly summarize the statistical fits in the main text of the paper (Table II). In addition to the coefficient of determination (R 2 ), mean absolute percent error (MAPE) is used to assess goodness of fit and is defined as where y is the variable on the vertical axis and M is the number of experimental cases (M = 179).

APPENDIX B: DECAY OF ENTRAINED FLUID
In Sec. IV D, we showed that the far-field wake of the drop seems to have minimal impact on retention and droplet delay. Instead, the gradual bleeding away of the local tail of fluid carried by 124803-23  Fig. 13(a) Retention, entrainment times t r = t e (0.47 U/U u ) 0.91 22% Fig. 13(b), Eq. (9) Entrainment timescale τ e = 38Fr −0.81 Re −0.21 u 0.77 34% Fig. 14(c), Fig. 16(a) Retention timescale τ r = 13Fr −1.5 Re −0.38 u 0.82 47% Fig. 15(c), Fig. 16 Fig. 18 the drop may play a dominant role in drop retention. In Figs. 5 and 6, the width of the tail of fluid dragged by the drop (denoted by changes in ambient illumination in the shadowgraph) slowly decreases over time. Although shadowgraphs are generally a qualitative tool for understanding density variations, we were able to estimate the approximate diameter of the trailing fluid carried by the drop as it rises through and past the transition region. These values were measured manually for the two cases shown here-as well as three others for which tracking and shadowgraph data were available-approximately one diameter below the bottom of the drop as it rises, as shown in Fig. 19(a). Because bright and dark regions of the shadowgraph indicate regions of strong concavity in the density field (i.e., large values of ∇ 2 ρ), actual perturbations in the density field persist slightly farther than can be observed in the shadowgraph. The estimated diameter is thus some fixed fraction of the actual wake diameter; however, we believe this is an adequate analog to examine trends in the wake over time. The ratio of wake diameter to drop diameter is plotted in Fig. 19(d) versus time nondimensionalized by the buoyancy frequency N. Nondimensionalized position and velocity for each case are also included in (b) and (c) for easy comparison. Time series of wake diameter are shorter than those of tracked position, as the shadowgraphs images were zoomed in closer to see fine details and the drop thus remained in the field of view of the camera for a shorter period of time.
An exponential decay can be fit to all five sets of data. Each time series was fit to a function of the form yielding an average best-fit decay coefficient of k = −0.60 ± 0.06 across all five cases. The time at which the drops have asymptoted to their upper-layer terminal velocity, denoted by the vertical colored lines in Figs. 19(b)-19(d), coincides with when the wake diameter is equal to between 1/4 and 1/2 of the drop diameter. Once the tail has become about two to four times smaller than the droplet, the drop reaches its homogeneous upper-layer behavior. The decay of this tail of denser fluid correlates well with the measured entrainment timescale for each case (i.e., the time between the vertical dotted gray lines in Figs. 19(b)-19(d) and the vertical lines of different shading), and supports our interpretation of this timescale as an indicator of fluid entrainment, albeit for a narrow range of Fr. While these visualization results are preliminary and limited to only a single drop density, we find that this is a promising path to explore further. We hope that a similar approach, over a wider range of Fr and Re u , could be used to physically understand the forces due to entrained fluid studied previously (including modified buoyancy force, as well as baroclinic vorticity force) [24][25][26]32,37,51], the majority of which have focused on the deceleration timescale t min and have not quantified the diameter of the trailing entrained fluid over time.

APPENDIX C: ORDER-OF-MAGNITUDE ESTIMATES OF ADDED MASS, HISTORY, AND MARANGONI FORCES
In Sec. V, we presented a simplified force balance at the point of velocity minimum, and used the scaling of the drag and buoyancy forces with velocity to estimate the excess force acting on the droplet due to interactions with the ambient stratification. Here, we develop order-of-magnitude estimates of the forces we neglected (the added mass and history forces) or incorporated into the stratification force (the Marangoni force).
The added mass and history forces result from the unsteady acceleration of the fluid surrounding the droplet as it changes velocity. These two forces depend on the acceleration of the droplet; the instantaneous acceleration cannot be accurately measured from our experimental data, and so the time dependence of these quantities is difficult to compute. While the added mass force should by its definition be zero at the point of velocity minimum (our focus in Sec. V), it is difficult to separate 124803-25 the added mass and history forces, and so here we give an order-of-magnitude estimate of their sum in two phases: in the deceleration preceding U min , and the more gradual acceleration following the velocity minimum. We focus on behavior surrounding the velocity minimum as this is where previous studies have shown the stratification force is at its maximum [25,32].
To estimate the orders of magnitude for various forces at the velocity minimum, we will use the experimental run shown in Fig. 5 as a reference case. For this run, we have h = 4.6 cm, d = 0.40 cm, ρ d = 0.9927 g/cc, ρ u = 0.9972 g/cc, ρ l = 1.054 g/cc. In general, we will assume ρ f = ρ u where needed, as the relevant dynamics take place when the drop is in the upper layer. We can estimate the magnitude of these forces following previous work [24,32,49], where the dimensionless acceleration parameter M a used in Refs. [24,32,49] is 0.1, so we have assumed the parameter values C a ≈ 2.1 and C h ≈ 0.48. In Eq. (C2), we have assumed the particle decelerates at a constant rate for a finite time, t, for a total velocity change of U . Significant acceleration occurs in two phases: (a) Phase 1: As the particle passes through the transition region, it decelerates rapidly (t = 2.5-3.5 s in Fig. 5, indicated by the shaded region). For our reference drop in this phase, we can estimate that U ∼ −8 cm/s and t ≈ 1 s, resulting in F a + F h ≈ 0.6 dyn. Although this is larger than the buoyancy force in the upper layer (0.14 dyn), the buoyancy force in the lower layer is much larger, 1.95 dyn.
(b) Phase 2: Just after the particle passes through the transition region, it reaches its minimum velocity and slowly returns to the upper layer terminal velocity (t = 4-8 s in Fig. 5); this was the primary delay process that the entrainment and retention times quantified (Sec. IV C). Here we estimate that U ≈ 1 cm/s and t ≈ 4 s, resulting in F a + F h ∼ −0.03 dyn.
Thus, the added mass and history forces are comparable to the buoyancy force only during the initial deceleration regime, as the droplet first interacts with the transition region. After the droplet passes through the transition region it slowly returns to the upper layer terminal velocity; here the acceleration is 1-2 orders of magnitude lower, and so the added mass and history forces are much smaller than other forces in the problem. Since the entrainment and retention times are primarly related to the gradual velocity recovery (Sec. IV C), we do not expect the added mass or history forces to have a significant effect on the overall retention or delay studied here.
We note three limitations when using these formulas: (1) these forces have previously been computed for solid particles, not droplets, (2) the forces are computed assuming homogeneous background fluids, and (3) the above equations are for somewhat lower Reynolds number, 0 < Re < 62. However, as we are not aware of any work on liquid droplets in this moderate Reynolds number regime, we believe this is the best estimate currently available. Due to the fact that the droplets are nearly spherical (e.g., Figs. 5 and 6), we believe the order of magnitude should at least be correct. Additionally, this analysis is consistent with previous work on solid particles passing through stratification layers, which have found the added mass and history forces are negligible, especially after the droplet has entered the final fluid layer [32], as well as previous work on solid particles in a linear stratification [25]. Finally, we find good agreement between our force enhancement ratio and the scaling of Zhang et al.'s baroclinic vorticity force F ρω , suggesting that we are indeed capturing primarily stratification-related forces in our estimate of .
The above estimates are for forces that occur for both solid particles and droplets. For liquid drops, the Marangoni force-due to the variation of interfacial tension between the droplet and surrounding fluid with salt concentration-can also play a role. In Sec. V we are treating F s as the remainder of the force balance; therefore, Marangoni forces are implicitly included in our measurements of F s , and in turn, . Nonetheless, we give an order-of-magnitude estimate of this force, with the caveat that this will vary depending on the salt used or the presence of surfactants.
In general, we expect different interfacial tension between the drop and the ambient fluid in the upper and lower layers, but these small changes are difficult to measure directly, and have not previously been characterized for silicone oil/salt water solutions. For hydrocarbon oils, the interfacial tension γ typically increases with NaCl concentration, m, at a rate of dγ /dm ∼ 1.5 dyne/cm/mol/L (although caution is warranted; the exact rate depends on chemical species, and in rare cases the interfacial tension can actually decrease) [53]. Assuming linear stratification, the Marangoni force is F M = (−dγ /dz)(π d 2 ), and so the maximum force can be estimated as where for our salt water solutions, the NaCl concentration in the lower layer is m l ≈ 1 mol/L (and m u ≈ 0 mol/L), and the above force is estimated for the reference case. This indicates that when the drop is in the transition region, the Marangoni forces can be comparable to or even slightly larger than the buoyancy. Note that (assuming the interfacial tension increases with concentration, which is true for hydrocarbon oils in NaCl solutions) this force would act in the same direction as buoyancy. Thus, it should decrease the retention time of the droplet. Regardless, the net effect of Marangoni forces, if they are significant, are included in the force enhancement . Unfortunately, the experimental setup did not allow for systematic variation of the surface tension gradient. Thus, we were not able to isolate the effect of Marangoni forces on the dynamics of the droplet, and exploration of this effect will be left for future work.