Retention of rising droplets in density stratification

TL;DR

This study experimentally investigates how single oil droplets rise in a two-layer density stratification, identifying key timescales and parameters influencing their retention, and introduces a theoretical model for drag enhancement that aligns with recent computational work.

Contribution

First experimental study on liquid droplet retention in stratification, providing quantitative timescales and a theoretical formulation for drag enhancement.

Findings

Retention timescales depend on Froude and Reynolds numbers.

Significant trapping occurs when Fr ≲ 1.

Theoretical model for drag enhancement matches recent computational results.

Abstract

In this study, we present results from experiments on the retention of single oil droplets rising through a two-layer density stratification. These experiments confirm the significant slowdown observed in past literature of settling and rising particles and droplets in stratification, and are the first experiments to study single liquid droplets as opposed to solid particles. By tracking the motion of the droplets as they rise through a stratified fluid, we identify two timescales which quantitatively describe this slowdown: an entrainment timescale, and a retention timescale. The entrainment timescale is a measure of the time that a droplet spends below its upper-layer terminal velocity and relates to the length of time 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. Both timescales are found to depend on the Froude and Reynolds numbers of the system, Fr $=U_u/(Nd)$ and Re $=\rho_u U_u d/\nu$. We find that both timescales are only significantly large for Fr $\lesssim1$, 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 drag enhancement $\Gamma$, the ratio between the maximum stratification 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 et al. [J. Fluid Mech. 875, 622-656 (2019)] on the drag enhancement on a solid sphere settling in a stratified fluid, and provides the first experimental data supporting their approach.

Tables

Table 1. Table 1 : Definition and range of parameters covered in laboratory experiments. Variables with the subscript f 𝑓 f represent that quantity in either the upper layer (e.g., ρ u subscript 𝜌 𝑢 \rho_{u} and U u subscript 𝑈 𝑢 U_{u} ) or lower layer (e.g., ρ l subscript 𝜌 𝑙 \rho_{l} and U l subscript 𝑈 𝑙 U_{l} ).

Parameter	Symbol	Definition	Range of values
Drop density	${\rho }_d$	-	0.9375 – 0.9927 g cm-3
Ambient fluid density	${\rho }_f$	-	0.9972 – 1.117 g cm-3
Terminal drop speed	$U_f$	-	0.33 – 13.3 cm s-1
Drop diameter	$d$	-	0.15 – 0.78 cm
Transition thickness	$h$	-	3.0 – 9.0 cm
Buoyancy frequency	$N$	$\sqrt{-\left(g/{\rho }_u\right)\left(\partial \rho /\partial z\right)}$	3.6 – 7.5 s-1
Dynamic viscosity of water	${\mu }_u$	-	0.01 g cm-1 s-1
Dynamic viscosity of oil	${\mu }_d$	-	0.093 – 0.10 g cm-1 s-1
Drop Reynolds number	Red	${\rho }_d{U}_{u}d/{\mu }_u$	5.4 – 540
Reynolds number	Ref	${\rho }_f{U}_{f}d/{\mu }_u$	5.4 – 1060
Archimedes number	Arf	$g\left({\rho }_f-{\rho }_d\right){\rho }_f{d}^3/{\mu }_u^2$	140 – 700,000
Froude number	Fr	$U_u/\left(Nd\right)$	0.38 – 4.2
Prandtl number	Pr	$\nu /\kappa$	$\sim 600$
Morton number	Mo	$g{\mu }_u^4\left({\rho }_f-{\rho }_d\right)/\left({\rho }_f^2{\gamma }^3\right)$	${10}^{-11}-{10}^{-12}$
Relative density	$\Delta {\rho }_u$	$\left({\rho }_u-{\rho }_d\right)/{\rho }_u$	0.0045 – 0.061
Relative layer thickness	-	$h/d$	4.5 – 52

Equations

$$\rho (z)={\rho }_u+\frac{2{\rho }_u}{g}\frac{1}{\beta L_{ruler}^2}{\int }_{z_u}^z\Delta zd{z}^{\prime }$$

$${\text{Ar}}_f=g\frac{({\rho }_f-{\rho }_d){\rho }_f{d}^3}{{\mu }_u^2}$$

$${\text{Re}}_f=\frac{{\rho }_f{U}_{f}d}{{\mu }_u}$$

$${\text{Re}}_f=0.60{\text{Ar}}_f^{0.57}$$

$${\log }_{10}\text{Re}=-1.81391+1.34671W-0.12427W^2+0.006344W^3,$$

$$\text{MAPE}=\frac{100\%}{M}\sum _{i=1}^M\left|\frac{y_{i,data}-y_{i,fit}}{y_{i,data}}\right|$$

$$C_D=\frac{4}{3}\frac{g({\rho }_f-{\rho }_d)d}{{\rho }_f{U}_f^2}$$

$$t_r=\frac{d_r}{U_u}$$

$$d_r\sim t_e\Delta U,$$

$$t_r=t_e\left(0.46\frac{\Delta U}{U_u}\right).$$

$$\frac{d_{\text{wake}}}{d}\propto \exp \left(kNt\right),$$

$$F_{tot}$$

$$F_{bh}({\rho }_f)=\frac{\pi }{6}\left({\rho }_f-{\rho }_d\right)g{d}^3=F_0{\text{Ar}}_f,$$

$${\text{Re}}_f\approx \alpha {\text{Ar}}_f^{\beta },$$

$${\text{Ar}}_f(U_f)={\left(\frac{U_f}{U_0}\right)}^{1/\beta }$$

$$F_{dh}\left(U_f\right)=-F_0{\left(\frac{U_f}{U_0}\right)}^{1/\beta }.$$

$$F_{s,max}=-F_b-F_d$$

$$F_b=F_{bh}({\rho }_u)=F_0{\text{Ar}}_u=F_0{\left(\frac{U_u}{U_0}\right)}^{1/\beta }.$$

$$F_d=F_{dh}(U_{min})=-F_0{\left(\frac{U_{min}}{U_0}\right)}^{1/\beta }.$$

$$\Gamma =\frac{F_{s,max}}{F_d}=\frac{-F_b-F_d}{F_d}=\frac{-F_b}{F_d}-1.$$

$$\Gamma ={\left(\frac{U_u}{U_{min}}\right)}^{1/\beta }-1.$$

$$F_{\rho \omega }\sim {\text{Fr}}^{-1}{\text{Re}}^{-1/2}.$$

Full text

I Introduction

There are many examples of droplets, bubbles, and particles interacting with stratified fluids, including atmospheric and marine pollution (Turco et al., 1990), oil spills (Camilli et al., 2010; Kessler et al., 2011; Socolofsky and Adams, 2002, 2005; Gros et al., 2017), oil seeps Clark et al. (2000); Leifer et al. (2009), falling leaves (Lam et al., 2019), and marine snow MacIntyre et al. (1995); Kindler et al. (2010); Prairie et al. (2013). As a result, understanding when and how stratification affects rising and settling is of significant interest within a variety of fields.

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 Horowitz and Williamson (2010); Auguste and Magnaudet (2018) to basic statistics such as terminal rise velocity of droplets Wegener et al. (2010); Bäumler et al. (2011); Bertakis et al. (2010). However, density stratification adds an additional level of complexity. When a droplet rises it entrains denser fluid with it, altering the effective buoyancy of the droplet and reducing its upward speed (figure 1). 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 Zhang et al. (2019), indicating the complex nature of this basic flow phenomenon.

Many prior experiments 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. Bayareh et al. (2013) showed that the drag coefficient of a settling spherical drop was enhanced in linearly stratified fluids with drop Froude numbers in the range $4\leq\text{Fr}\leq 16$ (where Fr is the ratio of the buoyancy timescale to the inertial timescale). Shaik and Ardekani found that at low Reynolds number, stratification and inertial forces both can increase the drag on a drop, and that this drag enhancement depends on the dynamic viscosity ratio between the drop and the surrounding fluid Shaik and Ardekani (2020). Other numerical work has studied two-droplet (Bayareh et al., 2016) and swarm-scale (Dabiri et al., 2015) 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 Blanchette and Shapiro (2012) 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. Srdić-Mitrović et al. (1999) 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. Abaid et al. Abaid et al. (2004) 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 Hanazaki et al. (2009); Okino et al. (2017); Akiyama et al. (2019) have studied the wake of spheres moving vertically 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.

Experimental and theoretical studies of settling porous spheres have also predicted increased drag or prolonged retention times at sharp density transitions (Camassa et al., 2009, 2010, 2013; Prairie et al., 2013; Panah et al., 2017), 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 $\ll 1$). Similar studies on spheres in linear stratification (Yick et al., 2009; Mehaddi et al., 2018) 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 (Socolofsky et al., 2011; Weber et al., 2012), indicating the necessity of studies to be performed beyond the low Reynolds number regime.

Finally, recent theoretical and computational work by Zhang et al. Zhang et al. (2019) addressed the “stratification drag” force that is commonly used as a catch-all for the contribution of stratification to changes in a sphere’s motion, including increased residence time or significant slowdown Srdić-Mitrović et al. (1999); Verso et al. (2019). The authors decomposed contributions of stratification to enhanced drag into two components: (1) entrainment 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 (nonalignment of density and pressure gradients), 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. Using laboratory experiments, we examine flow and retention for a range of drop sizes, drop densities, and ambient stratification profiles. In section II, we discuss the nondimensional parameters relevant to this problem. In section III, we will describe the experimental setup and measurements taken. We will discuss our results in section IV, beginning with analysis of the drop’s position and velocity, and then introduce two timescales relating to fluid entrainment and droplet retention. We find that these timescales are related, and that they are dependent upon the Froude and Reynolds numbers of the system, with significant fluid entrainment or drop retention only occurring for Fr $\lesssim 1$. We then develop a theoretical formulation for the drag enhancement, $\Gamma$, induced by stratification. The scaling of $\Gamma$ compares very favorably to the work of Zhang et al. Zhang et al. (2019) and provides the first experimental evidence supporting their approach. Finally, we will close with a discussion of the implications of this work and future directions in section V.

II Nondimensional parameters

We covered a range of the parameter space relevant to this problem, particularly in the intermediate Reynolds number regime. Table 1 lists the parameter definitions and ranges covered in this study, which spanned 179 different droplet experiments. Non-dimensional parameters with the subscript $f$, such as Ref, encompass two separate nondimensional numbers for a single drop’s behavior in the upper and lower layers of ambient fluid, Reu and Rel. 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, $\rho_{f}$ represents fluid density, $U_{f}$ represents the terminal drop velocity in a given homogeneous-density region, $d$ is the drop diameter, $\mu$ is the dynamic 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. The buoyancy frequency, $N$, is defined as $N=\sqrt{-g/\rho_{u}(\partial\rho/\partial z)}$, where $\partial\rho/\partial 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 1, the Reynolds number (Ref or Red) represents the ratio of inertial to viscous forces. The Archimedes number (Arf) 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$). The Prandtl number, the ratio of momentum diffusivity to salt diffusivity, remained fixed at approximately 600 for all experiments. The Morton number is used to characterize the shape of a drop or bubble in a surrounding fluid, and is a function of the interfacial tension between the droplet and surrounding fluid, $\gamma$ (on the order of 30 mN/m for our experiments).

III Experimental approach

III.1 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 figure 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 ($\rho_{l}=1.106$ to $1.117$ g/cm3) and then with fresh water ($\rho_{u}=0.9972$ to $0.9981$ g/cm3), 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 hours, 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 hours. This two-layer filling method yielded error-function shaped density profiles, as seen in figure 3(a,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 programmed to yield a range of transition layer thicknesses. Examples of such density profiles are shown in figure 3(c,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) in order to study a range of drop densities $\rho_{d}$, from 0.9375 to 0.9927 g/cm3. 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 minutes between drop releases was chosen to ensure the tank was quiescent for each experiment.

III.2 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 high-speed 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 figure 4. In some cases, the drops exhibited slight out-of-plane motion, which we estimate to contribute 1 percent or less error in measured vertical position, based on a camera distance of $\sim 1$ meter 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 (Allan et al., 2016), which uses center-of-mass detection to determine droplet position. Example tracked paths and velocities for five drops are shown in figures 5, 6 and 7. Instantaneous velocities were obtained following the methods of Srdic-Mitrovic et al. Srdić-Mitrović et al. (1999), in which a least-squares line was fit to a window of 7 points of vertical position ($\sim$0.06 sec 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 (Sutherland et al., 1999; Dalziel et al., 2000). 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 $\Delta 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 figure 3) was obtained by integrating the gradient from the known upper layer density using the equation below,

[TABLE]

where $\Delta z$ is the measured apparent displacement field (the difference between the curve shown in figure 4(b) and vertical height), and $\beta\simeq 1.88$ s2 cm*-1* following Sutherland et al. (1999). 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 $\pm 0.8$ cm to yield a profile whose constant upper and lower layer densities matched those measured with the hand-held densitometer.

Finally, for 5 of the above experimental cases, shadowgraph experiments were performed to visualize the wakes of droplets. A schematic of this setup is shown in figure 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 in order to avoid total internal reflection within the transition region, which would have occluded $\sim$ 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 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 minutes 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 discuss the terminal velocity of the droplets in the upper and lower layers, as a nondimensional predictor of this behavior will prove useful when discussing drag enhancement in stratification in section IV.5. In the next two subsections, we analyze the timescales over which 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 connect these timescales to flow visualizations of the droplet’s wake. Finally, we present a theoretical formulation of the drag enhancement based on our experimental data and compare this with the results of Zhang et al. Zhang et al. (2019).

IV.1 Drop paths and velocities

Shown in figures 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 ($\rho_{d}=0.9927$ g/cm3) and rising through a 4.6 cm transition region. Relative shading indicates variation in the second derivative of density (Settles, 2012). 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 grey 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 section IV.4.

To demonstrate how varying different experimental parameters affect drop motion, figure 7 shows sample drop paths and velocities over time for four example experimental droplet cases. The first case, shown in figure 7(a,e), is the droplet shown in figure 6. We use this small, dense droplet passing through a 4.6 cm transition region as a reference case for comparison with: (b,f) a similarly-sized, lighter droplet in similar ambient stratification; (c,g) a larger, dense droplet in similar stratification; and (d,h) 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) 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 section IV.3.

IV.2 Terminal velocity

Droplets reach a constant, terminal velocity in both the upper and lower layers. This terminal velocity is governed by buoyancy, viscosity, and inertia. The drop’s Archimedes number in the upper and lower layer is plotted versus its Reynolds number in each layer in figures 8(a,b). As noted in table 1, the definitions of these two nondimensional numbers is as follows:

[TABLE]

[TABLE]

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 Arf and Ref. Measurements of Arl and Rel are represented as black circles, while values of Aru and Reu are shown as blue triangles.

A power-law fit was found to describe the relationship between these two parameters from experimental measurements:

[TABLE]

with $R^{2}=0.99$. This compares reasonably well with the relationship for solid spheres in the range $435<$ Ar $\leq 1.16\times 10^{7}$ and $12.2<$ Re $\leq 6.35\times 10^{3}$ Clift et al. (2005):

[TABLE]

where $W=\log_{10}(4/3$Ar), shown as the dot-dashed line in figures 8(a,b). It should be noted that this empirical equation for terminal velocity is for a solid sphere, not a deformable liquid drop.

In figure 8, as well as later figures, we assessed deviation between data and the empirical fits using the mean absolute percent error (MAPE), defined as

[TABLE]

where $y$ is the variable on the vertical axis and $M$ is the number of experimental cases ($M=179$).

Also shown in figures 8(a,b) are terminal velocities from experiments and numerical simulations of drops rising in homogeneous-density water from several other studies with varying Morton number. The Morton number (table 1) describes the importance of gravity and viscosity relative to interfacial tension and inertia, and is used to characterize the shape of drops and bubbles. Shown are velocities for toluene droplets (Mo = $1.95\times 10^{-11}$) Wegener et al. (2010), n-butanol acetate droplets (Mo = $2.87\times 10^{-10}$) Bäumler et al. (2011), and n-butanol droplets (Mo = $1.23\times 10^{-6}$) Bertakis et al. (2010). For some of these cases in the literature, the drop aspect ratio (the ratio of vertical drop diameter to horizontal drop diameter) was much less than 1 due to a higher Morton number, which explains some of the divergence from our experimental data, which we estimate to have a Morton number on the order of $10^{-11}$ to $10^{-12}$. Although there exist many more studies that compute a terminal velocity relation for droplets or nonspherical solids in homogeneous fluid (e.g. Wallis, 1974; Ganguly, 1990), the majority are for low Reynolds number flows and not applicable here.

Figure 8(c) also shows the approximate drag coefficient of drops in each homogeneous-density layer, computed as

[TABLE]

following Clift et al. (2005). The observed drag coefficient tends to be higher than that predicted by the standard drag curve for solid spheres (Clift et al., 2005, Table 5.2).

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 equation 4. This terminal velocity will be used later in section IV.5 for our theoretical formulation of drag enhancement in stratification.

IV.3 Entrainment and retention timescales

As shown in figures 5, 6, and 7, some droplets experienced a significant slowdown as they passed through the transition region. Previous work has shown that this slowdown is in part due to the entrainment of denser fluid Srdić-Mitrović et al. (1999). It is thus useful to look at metrics of drop retention that measure both the length of time that denser fluid is entrained, as well as the length of time that the drop is physically retained in the transition layer.

We will consider the first metric to be the entrainment time, $t_{e}$. This timescale is demonstrated in figure 9, using both vertical drop position and drop velocity. The entrainment time measures the amount of time the drop spends below its upper layer terminal velocity $U_{u}$, i.e., the length of time over which entrained denser fluid is appreciably slowing the droplet’s motion. It is computed as the time between when the drop first slows to $U_{u}$ (followed by significant further slowdown), and when the drop again speeds up to this upper-layer terminal velocity. The point where the drop has asymptotically reached $U_{u}$ again is chosen as the point at which the distance between the drop position (the black line in figure 9(a)) and a line representing the upper layer terminal speed (the grey dashed line in figure 9(a)) is less than 5 percent of the drop diameter. This timescale is very similar to the delayed settling time (DST) used by Prairie et al. Prairie et al. (2013).

The second metric is a retention time, $t_{r}$, the time that the droplet is retained in the transition layer. This time is shown in figure 10(a), and is computed as the time between when the drop actually passes an upper threshold ($z\approx 10$ cm, where $z=0$ is the center of the transition region), and when it would have passed the upper threshold if it had not slowed down once it first reached $U_{u}$. Figure 10(b) compares values of $t_{e}$ and $t_{r}$, which follow a power-law relationship of $t_{r}=0.079t_{e}^{1.7}$ with an $R^{2}$ value of 0.93. There is obviously a strong relationship between these two metrics, with droplets that have a long entrainment timescale also being retained the longest.

The retention time can also be expressed in terms of the retention distance, $d_{r}$, as shown in figure 10(a), and the upper layer terminal velocity:

[TABLE]

and $d_{r}$ can be scaled as

[TABLE]

where $\Delta U=U_{u}-U_{\text{min}}$ (figure 9(b)). We found a best-fit relation of $d_{r}=0.46t_{e}\Delta U$, with an $R^{2}$ value of 0.81, as shown in figure 11(a). This expression can then be substituted into equation 8, yielding

[TABLE]

The retention time is, therefore, the entrainment time multiplied by a constant, $c=0.46$, and by a factor representing the relative magnitude of the drop’s slowdown, $\Delta U/U_{u}$. This relationship is shown in figure 11(b), with an $R^{2}$ value of 0.92. The physical retention of a drop at the transition region is thus a function of the length of time that denser fluid is appreciably entrained, $t_{e}$, and the relative drop slowdown, $\Delta U/U_{u}$, which may be due to both the amount of dense fluid entrained as well as any contributions of baroclinic vorticity generation to the stratification force.

In order 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 12 and 13 show nondimensional entrainment time, $\tau_{e}=t_{e}N$, and nondimensional retention time, $\tau_{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, indicating the importance of stratification in drop retention. 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 figures 12(c) and 13(c), are plotted against Fr in figure 14. 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 $\approx 1$, and we can split droplet behavior into two regimes:

For Fr $\geq 1$, the drop is rising on a timescale less than the buoyancy timescale, and the drop behaves as if it were in a homogeneous fluid. Retention times are close to zero, as inertial forces dominate over buoyant forces. 2. 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 Froude number is the most important parameter to consider when predicting trapping or slowdown in a two-layer stratification.

We observed a significant slowdown for a wide range of Reynolds number, with $\tau_{e}>10$ observed for Reu up to 300. Significant delay of this magnitude was only observed for Fr $\lesssim 1$ and and $\Delta\rho_{u}\lesssim 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. Srdić-Mitrović et al. (1999) 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. Torres et al. (2000) 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. (Yick et al., 2009), found that enhanced drag scales with Ri0.51 in the very small Reynolds number regime, where the Richardson number Ri $=\frac{N}{\nu/d^{2}}\frac{1}{\text{Fr}}$. Recent work by Zhang et al. Zhang et al. (2019) contributed the most exhaustive study of the parameter space to date, and delineated the dominant mechanisms of drag enhancement within different regimes. Their modeling and numerical results have not been experimentally verified, and did not quantify delay in terms of time scales, an aspect of this increased stratification drag that is particularly relevant in environmental applications. We will discuss the applicability of our results to their work on drag enhancement in more detail in section IV.5.

IV.4 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 discuss visualizations of the wakes of the droplets first shown in figures 5 and 6, and show that there is correlation between the decay of the trailing denser fluid and the entrainment timescale.

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 ($\rho_{d}$ = 0.9927 g/cm3) 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 Horowitz and Williamson (2010), 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 Hanazaki et al. (2009); 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 sec; $\tau_{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.

Although the far-field wake of the drop seems to have no impact, the gradual bleeding away of the local tail of fluid carried by the drop may instead play a role in drop retention. In figures 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 figure 15(a). Because bright and dark regions of the shadowgraph indicate regions of strong concavity in the density field (i.e., large values of $\nabla^{2}\rho$), 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 figure 15(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

[TABLE]

yielding an average best-fit decay coefficient of $k=-0.60\pm 0.06$ across all five cases, with an average $R^{2}$ of 0.97. The time at which the drops have asymptoted to their upper-layer terminal velocity, denoted by the vertical colored lines in figure 15(b-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 grey lines in figure 15(b-d) and the vertical lines of different shading), and supports our interpretation of this timescale as an indicator of fluid entrainment.

To our knowledge, this is the first experimental work to quantitatively measure the wake of a spheroid in a density transition, and relate it to the spheroid’s slowdown and retention. 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, and hope that a similar approach could be used to physically understand the force due to entrained fluid studied previously Srdić-Mitrović et al. (1999); Torres et al. (2000); Yick et al. (2009); Verso et al. (2019); Zhang et al. (2019).

IV.5 Comparison to other work on drag enhancement in stratification

While we have framed our observations above in terms of timescales, the experimental data from this study also allow for computation of enhanced drag forces in the presence of stratification. Here, we will present a theoretical formulation for this drag enhancement and compare it to the parameter space and force balances presented by Zhang et al. Zhang et al. (2019) for numerical studies of rigid spheres.

The forces acting on the particles or droplets passing through a stratified fluid are typically decomposed into several different terms Clift et al. (2005); Srdić-Mitrović et al. (1999); Verso et al. (2019):

[TABLE]

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. 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 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 Verso et al. (2019).

First, we will examine the functional form of the buoyancy and drag forces in a homogeneous fluid, $F_{bh}$ and $F_{dh}$, far away from the transition region. The buoyancy force can be computed as a function of the local fluid density $\rho_{f}$, which varies with depth $z$ Srdić-Mitrović et al. (1999):

[TABLE]

where $F_{0}=\frac{\pi\mu^{2}}{6\rho_{f}}$. We can relate this buoyancy force to the local terminal velocity, $U_{f}(z)$, using the relationship measured in section IV.2,

[TABLE]

and the definition of Reynolds number, Re${}_{f}=\rho_{f}U_{f}d/\mu$, and rearrange to obtain an expression for the Archimedes number in terms of the local droplet velocity,

[TABLE]

where $U_{0}=\frac{\alpha\mu}{\rho_{f}d}$, and the empirical measurements discussed in section IV.2 gave $\alpha=0.60$ and $\beta=0.57$.

In a homogeneous fluid, a droplet will reach its terminal velocity when $F_{dh}=-F_{bh}$. Using this relation and equations 13 and 15 we can obtain a drag law for the drag force in a homogeneous fluid as

[TABLE]

We now use the relations given in equations 13, 15, and 16 in a simplified force balance in a stratified fluid. The force balance is considerably more tractable if we focus on the force enhancement at the point of minimal velocity ($t=t_{min}$, $u=U_{min}$), which is when the stratification force is at its maximum Verso et al. (2019). At this point, $\frac{\partial u}{\partial t}=0$, $F_{a}$ is a function of $\frac{\partial u}{\partial t}$ and so is also zero, we assume $F_{h}$ is negligible Verso et al. (2019), and the droplet is usually entirely in the upper fluid layer (see e.g. figures 5, 6, 7). The force balance from equation 17 can thus be assumed to simplify to

[TABLE]

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 equations 13 and 15:

[TABLE]

Using equation 16, the local drag force scales with the local velocity to the power of $(1/\beta)$,

[TABLE]

Thus, we can write a drag enhancement ratio, $\Gamma$, which compares the maximum stratification force experienced by the droplet to the drag force:

[TABLE]

Substituting the expressions from equations 18 and 19, we find that

[TABLE]

This drag ratio $\Gamma$ is shown in figure 16 as a function of Reu and Fr. A power-law fit was found, with $\Gamma=19$ Fr*-0.90* Re${}_{u}^{-0.55}$. These exponents agree remarkably well with the dependence of the drag force due to vorticity induced by baroclinic torque, $F_{\rho\omega}$, proposed by Zhang et al. Zhang et al. (2019). In their Regimes 2 and 3 (Pr${}^{-1/3}\ll$ Fr $\ll 1$ and Fr $\gg 1$, respectively, and in which all of our experimental measurements fall), they found that

[TABLE]

In figure 16, the largest errors in $\Gamma$ occur for the smallest Froude numbers, which fall towards the border of Regime 1 (for large Re and Pr, Regime 1 is delineated as Fr $\ll$ Pr*-1/6*, or Fr $\ll 0.34$ in our case). Our smallest Froude number is 0.38. The scaling of drag enhancement from our experimental data therefore agrees best where we expect it to: farther from the edges of the regimes delineated in Zhang et al. Zhang et al. (2019).

In summary, these experimental and analytical results regarding the drag enhancement due to stratification show that the baroclinic vorticity contribution is important in our experiments, and provide the first experimental data supporting the approach of Zhang et al Zhang et al. (2019). 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_{\rho\rho}$), and the force due to modification of the local vorticity field due to baroclinic torque ($F_{\rho\omega}$) – are predicted to scale with Fr*-1* in the parameter space covered by our experiments. Zhang et al. showed that $F_{\rho\omega}$ should dominate in the Reynolds number regime studied here, which results in our observed scaling of $\Gamma$ with both Fr and Reu. We do expect $F_{\rho\rho}$ 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_{\rho\omega}$. 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_{\rho\omega}$ is significant in the regime studied here.

V 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 $\Delta\rho_{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. Srdić-Mitrović et al. (1999), who found that experiments at higher Reu ($>15$) showed no significant change in drag as a sphere settled through a density gradient, we observed significant drop retention and slowdown for Reu up to 300, and for a wide range of Rel. As noted in previous work Abaid et al. (2004); Verso et al. (2019), 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.

Although our results regarding decay of a droplet’s wake are promising, a more systematic study is required to draw stronger conclusions about the relationship between the entrainment timescale and the diameter, length, and overall shape of the trailing fluid column, as well as about the vortex shedding that happens before and during passage through the transition region. 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. However, these forces are implicitly included in $\Gamma$, so any interfacial tension effects are captured in our measurement of the additional stratification-related forces. Our results regarding the scaling of drag 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 (Adams and Weatherly, 1981) and dispersal of pollutants in the atmosphere (Turco et al., 1990).

VI Conclusions

In this study, we characterized the dynamics governing retention of a single droplet at a transition in density between two homogeneous fluids. 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}\sim$ Ar${}_{f}^{\beta}$, where $\beta=0.57$.

We explored two metrics measuring the timescale of drop delay at a density transition. 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, $t_{r}=0.46t_{e}(\Delta U/U_{u})$, where $\Delta 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<\Delta\rho_{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 $\lesssim 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 the trailing column of dense fluid entrained by a droplet decays over time. The point at which the wake diameter becomes on the order of 1/4 to 1/2 of the drop diameter coincides with the drop reaching its upper-layer terminal velocity, supporting our interpretation of the entrainment timescale as the time required for denser trailing fluid to bleed away. This effect appears to be independent of the type of large-scale wake (zigzagging or vertical) the droplet has when first entering the transition region.

Finally, based on a force balance at the point of the droplet’s minimum velocity, we derived a drag enhancement ratio $\Gamma=\frac{F_{s,max}}{F_{d}}=(U_{u}/U_{min})^{1/\beta}-1$, and found a power-law fit to our experimental data of $\Gamma\sim$ Fr*-0.90*Re${}_{u}^{-0.55}$. This compares very favorably with the behavior of the baroclinic vorticity force proposed by Zhang et al., and suggests that drag forces generated by the baroclinic torque, likely in addition to the more traditionally-studied buoyancy-related entrainment forces, are important in this problem.

Acknowledgements.

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.