Gravity driven current during the coalescence of two sessile drops

Coalescence of liquid drops is critical in many phenomena such as emulsion stability, inkjet printing, and coating applications. For sessile drops on a solid surface, the coalescence process is more complicated than the coalescence of drops suspended in a fluid medium as a result of the coupling of the contact line motions to the fluid flow. In this paper, we use video microscopy to track the evolution of the interfaces and contact lines as well as the internal fluid motion within a merged sessile droplet. In this study, the fluids in the coalescing drops are miscible and have similar surface tensions and drop volumes but different viscosities and densities. Coalescence occurs in three stages. During the first stage, rapid healing of the bridge between the drops occurs just after they touch. In the second stage, slower rearrangement of the liquids occurs. We show that these intermediate rearrangements are driven by gravity even for density differences of the two fluids as small as 1%. For the systems exam...


INTRODUCTION
Coalescence of drops is critically important to numerous technological and natural processes.2][3] While the coalescence of drops suspended in a fluid medium has been examined extensively, 4,5 coalescence of two sessile drops on a solid surface adds the feature that the contact lines of the drops must move during the merging process. 6hus, contact line dynamics add richness and complexity to the coalescence of sessile drops.
Coalescence of sessile drops of identical fluids shows a rapid bridge healing process and a slow contact line relaxation.0][11] In the slow process, contact lines move and the drop relaxes to a circular shape over longer time scales, slowed by the dissipation at the moving contact lines and fluid flow. 9,10,12][20] Different mechanisms have been used to cause the drops to touch and coalesce.Drops are allowed to spread spontaneously into one another, 7,19 one drop is deposited onto another, 13,18,20 or one drop is blown into another 17 on a surface until they meet and merge.11]16 In still other cases, drops are forced to collide due to wettability gradients on the surface, [21][22][23] where wettability gradients cause the merging drops to exhibit different contact angles just prior to coalescence even when the drops contain identical liquids.
Analyses of the initial bridge healing have included scaling arguments as well as numerical simulation based on a lubrication approximation. 15,24More detailed numerical simulations have used the volume-of-fluid method 25 and lattice Boltzmann method. 18These studies demonstrate two main features of the initial stage of coalescence.First, the flow during bridge healing can be divided into two regimes: a visco-capillary and an inertio-capillary regime 11,13,19,26,27 depending on the relative importance of viscosity, inertia, and capillarity.Second, the height and width of the bridge grow with power law behavior in time, where the power law exponents depend on whether viscous stresses or inertia are comparable to capillarity. Beysens and coworkers found that the relaxation of the drop to a spherical shape is usually five or six orders of magnitude slower than coalescence of drops in a bulk fluid because of the contact line motion. 6,9,10This slow relaxation process is attributed to a very small Arrhenius factor resulting from a liquid-vapor phase change in the vicinity of the contact line.
1][32] In this case, coalescence is strongly delayed by Marangoni forces. 33,34The delayed coalescence behavior is also observed in the merging of ethanol and water drops, 35 where ethanol evaporation leads to self-propelling drops.Delayed coalescence behaviors are unique to situations involving significant Marangoni forces and do not occur when the asymmetry is due to different deposition times for the same fluid. 36n microfluidic devices, advective mixing during coalescence is used to control chemical reactions and physical dispersions. 37As suggested by Ottino and coworkers, 38,39 the design of chaotic micromixers can enable effective mixing even at low Reynolds number.In the confined geometries of microfluidic devices, advective mixing after coalescence has been achieved using herringbone structures in the channels 40 as well as using electro-wetting 41 or electro-kinetic based mixing. 42n contrast, the efficiency of mixing in open geometries such as sessile drops is less well investigated.Experiments examining coalescence of sessile drops, supplemented by simulation, show little mixing without applying wetting conditions or external forces. 432][23] Electro-wetting induced drop oscillation by application of an external AC field can enhance mixing. 44n this paper, we investigate coalescence of two sessile drops containing fluids with different densities and viscosities, but nearly the same contact angle and volume, which results in the drops having the same shape just before they coalesce.Adapting a laser induced fluorescence method first used to visualize a drop coalescing with a bulk fluid, 45 we track the mobile interface between the two merging fluids as well as the composite drop/air interface during the coalescence process.With simultaneous side and top views, the contact line shape is also tracked.Similar to previous studies, a rapid bridge healing process is observed, which we term stage 1.This stage of the coalescence process is governed by minimization of the capillary energy for drops with either the same densities or different densities.During the slower contact line relaxation process, we see that there is internal fluid motion that occurs despite imperceptible changes in the shape of the external interface.The internal interface formed between the two initial liquid volumes demonstrates an advective motion for initial drops with density difference as small as 1%.7][48][49] While the rapid bridge healing process has been studied in detail before, internal advective motion over longer time scales has not.The advective motion observed in stage 2 of the coalescence process is the primary focus of this investigation.Using dimensional analysis as well as a lubrication formulation, we identify the characteristic time scale and the dimensionless groups controlling the intermediate, gravity-driven fluid motion of stage 2. We also note that there is a final stage of coalescence, stage 3, in which diffusive motion leads to complete mixing of the fluid within the merged drop.For the systems we examine, the time scales for the three stages of coalescence are well separated.Further, during stages 1 and 3, there is little advective motion for the systems we examine, while the bulk of the advective fluid motion occurs during stage 2.

MATERIALS AND METHODS
The substrates are fabricated from polydimethylsiloxane (PDMS) using a standard mixture ratio of 10:1 for the silicone elastomer base and the elastomer curing agent (Dow Sylgard 184).The PDMS is first mixed and defoamed using a centrifugal mixer (Thinky AR100) before pouring the mixture into Petri dishes to cure and form a flat surface.Once poured, the mixture is degassed in the Petri dish for about 10 min and cured in a 60 • C oven overnight.The PDMS sample is removed from the Petri dish and the flat surface is used as the solid substrate surface for the coalescence experiment.Before experiments are performed, the solid PDMS substrate is rinsed with ethanol, acetone, and distilled and deionized (DI) water.The measured advancing and receding water contact angles on the PDMS surfaces are θ a = 105 • ± 9 • and θ r = 62 • ± 10 • , respectively.These contact angle values are in agreement with previously reported measurements on similar surfaces and suggest that a ridge forms in the soft PDMS below the contact line. 50If force is applied to the contact line on such soft surfaces, the contact line moves, pulling the ridge along but with enhanced dissipation. 51he fluids tested are composed of mixtures of purified water (Milli-Q, 18 MΩ cm, organic content < 10 ppb) with varying concentrations of glycerol (Sigma-Aldrich ≥ 99.5%), allowing the adjustment of the mixture viscosity with minimal change in surface tension.NaCl (Sigma-Aldrich ≥ 99.5%) and CsCl (Sigma-Aldrich ≥ 99%) are used to adjust the density of the glycerol/water mixtures with minimal viscosity and surface tension changes.Fluorescein dye (0.05 mM, Acros Organics) is used as a fluorescent marker in one of the fluid drops to be merged in the coalescence experiment.The solutions carrying the fluorescein are buffered to pH = 9.00 ± 0.02 (pHydrion Buffers, Micro Essentials) to maintain a high quantum yield of the dye.To ensure that the buffer salt does not alter the fluid properties with and without dye, all the aqueous fluids are prepared using the same pH 9 buffer solution.Even though pure glycerol is hygroscopic, the further absorption of water in the prepared solutions is observed to be slow.We store the solutions in capped beakers and periodically measure the densities, which do not change for weeks.
The density, viscosity, and surface tension of each solution used are listed in Table I.Densities are determined by measuring the mass (Denver Instrumental Company Balance, Model XE-100A, accuracy 10 −4 g) of a known volume contained in a volumetric flask (Pyrex No. 5640, 25 ± 0.03 ml).The viscosities of the solutions are determined using a Cannon-Fenske viscometer (Model K655-450) following the standard procedure. 52urface tension is measured using a Wilhelmy pin connected to a balance on a MicroTroughX Langmuir trough. 53The Pt pin is washed following a standard cleaning protocol and burned with a Bunsen burner to remove residual organics, after which the Pt pin is placed back on the Langmuir trough balance and the balance is zeroed in the air environment.The apparatus is calibrated using the known PDMS surface tension (19.8 mN/m) to obtain the calibration constant and again checked for consistency by measuring the surface tension of DI water.The uncertainties in the surface  tension values reported in Table I arise from statistical fluctuation in the measurement of both the fluid surface tension and the calibration factor.
To probe the influence of asymmetric properties on the coalescence behavior, we select fluid pairs with various density differences as well as different viscosities.The fluid pairs used in the experiments are listed in Table II.The characteristic parameters are density difference ∆ρ = | ρ b − ρ t | and viscosity ratio λ = µ t /µ b , where the subscript t denotes the fluid property of the smaller density fluid and the subscript b denotes the fluid property of the larger density fluid.The Bond number, Bo = ρ t gL 2 /γ (L is the length of the composite drop), describing the relative importance of gravity versus surface tension is also listed in Table II, with the range of values corresponding to those observed for the composite drop after coalescence.The surface tensions of the two fluids are similar.To be consistent, we use the surface tension of the lighter fluid to calculate the Bond number.As seen in Table II, density differences range from 1% to almost 40%, and the viscosity ratio varies from 0.2 to 10 such that the larger density fluid can be either more or less viscous than the smaller density fluid.Bond numbers range from 1 to 10, suggesting that gravity will flatten the external interface of the composite drop.
The fluid feed system is shown schematically in Figure 1(a).Flat-tipped syringe needles are pushed through the PDMS substrate, with the spacing between two feed points controlling the sizes of the drops just prior to coalescence.Fluid drops are pumped from underneath the substrate into the syringe needles which are connected to two syringes by polyvinyl chloride tubes.The syringes are attached to the same micrometer stage, ensuring that the two initial drops are of nearly equal volume as they touch and merge.The rate of drop formation is controlled manually by the micrometer at approximately 5-20 µl/min, producing contact line speeds of approximately 0.02-0.05mm/s.These contact line speeds are maintained until just before the drops touch, and then the final approach is achieved by natural relaxation of the contact lines after pumping has

EXPERIMENTAL RESULTS
Observations of coalescence of the fluid drop pairs listed in Table II are captured using the optical setup and video capture system described in Materials and Methods section.The resulting image sequences are used to demonstrate that there are three primary stages of coalescence, with each stage governed by different forces and phenomena.In Figure 3, frames (a.1) through (a.5) show coalescence of two drops with the same density.One of the initial drops contains fluorescent dye.In Figure 3, frames (b.1) through (b.5) show coalescence of two drops with a density difference of 0.074 g/ml.The lighter fluid contains fluorescent dye.In both sequences, the first frame shown is just prior to contact between the drops.The second frame shown corresponds to t = 0.03 s after the drops touch.By this time, the rapid bridge healing process is nearly complete in all experiments.In some experiments, we observe capillary waves resulting from the rapid bridge healing process.Frames (a.3) and (b.3) capture the end of the stage 1 bridge healing process, after which any capillary waves observed have decayed and the exterior shape of the merged drop no longer changes.The internal interface between the two merged fluids at the end of stage 1 is sharp and nearly vertical.The density differences do not change the bridge healing process during stage 1.
Based on the fluorescence imaging, we see no evidence of twisting of streamlines and therefore conclude that no advective mixing of the two initial fluids has occurred.Following the end of stage 1, frames (a.3) through (a.5) show that there is no internal flow within the merged drops of fluid pair BB within the first 3 s after the drops contact one another.Frames (b.3) through (b.5) show that the internal interface between fluid pair CB evolves over the first 3 s after contact until the interface is nearly horizontal, with the lighter, fluorescently dyed fluid stratified in the top layer.These two experiments taken together suggest that two fluids of different densities stratify due to a gravity current during stage 2 of the coalescence process.In stage 2, the gravity driven flow proceeds in time scales of the order of seconds until full stratification occurs, as shown in frame (b.5).At this late time in stage 2, the interface between the two initial fluid volumes is still sharp.As we will discuss later, at longer times, of the order of several minutes, stage 3 of the coalescence process occurs, in which diffusion blurs the interface and the two fluids mix.Throughout both stages 2 and 3 of coalescence, the contact line pinning, due to residual contact angle hysteresis, inhibits the composite drop from attaining the circular contact line footprint that surface tension alone would produce.

Stage 1-Bridge healing
For drop coalescence, the flow during bridge healing can be divided into two regimes: a visco-capillary and an inertio-capillary regime. 11,13,19,26The visco-capillary time scale is Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions.Downloaded to IP: 128.2.20.8On: Mon, 29 Feb 2016 18:15:58 τ ν ∼ µ 3 /ργ 2 ∼ 70 ns for a typical composite drop (µ ∼ 7 mPa s, ρ ∼ 1.2 g/ml, γ ∼ 70 mN/m) and ranges from 10 ns to 12 µs for all the experiments we performed.The inertio-capillary time scale is τ i =  ρR 3 0 /γ ∼ 50 ms for a typical composite drop (ρ ∼ 1.2 g/ml, R 0 ∼ 5 mm, γ ∼ 70 mN/m) and ranges from 10 ms to 250 ms for the fluid pairs we considered.In contrast to previous experiments on coalescence of drops with identical densities, 11 the fluids we examine are 20 times less viscous, resulting in a visco-capillary time scale 8000 times shorter than the previous experiments.Thus, the initial bridge-healing regime is dominated by capillarity and viscosity up to about 70 ns, which is significantly shorter than the camera frame rate and thus is not captured in our experiments.Beyond about 70 ns, the flow enters the inertio-capillary regime for time scales up to about 50 ms for drops a few millimeters in radius.Consistent with this estimate of the characteristic time scale, for most of our experimental conditions, the bridge between the coalescing drops heals and any capillary waves that appear on the drop surface decay within the first two video frames (≈60 ms).As pointed out in a recent experiment by Eddi, 27 capillary waves are observed for low viscosity drops which implies that the observations take place in the inertio-capillary regime.For very large drops and for the case with the largest density difference, observation of the side view of the merged droplet shows that drop shape oscillation can persist as long as 0.5 s.
During stage 1, the contact line at the bridge moves outward from the contact point while the contact lines along the long axis of the drop move slightly inward.At the end of stage 1, the internal interface between the fluids is sharp and tilts at most about 2 • away from vertical (Figure 3(b.3)).The conditions at the end of stage 1 set the initial condition for the gravity driven flow that occurs during stage 2 of the coalescence process.

Stage 2-Gravity current
In stage 2 of the coalescence process, the lighter fluid layer flows over the top of the denser fluid layer, even for fluid pairs AI and AH where the surface tension differences between these two fluids would have led to an initial Marangoni force opposing the gravity-driven flow and for fluid pair BJ where the density difference is only 1%.Experimentally, for all but the fluid pairs with the largest density differences, the shape of the external interface of the composite drop does not change appreciably after the end of stage 1 (see Figure 3(b.3) and later frames).
Figure 4(a) overlays edge profiles extracted from top view images near the beginning and end of stage 2. The nearly complete overlap of these edge profiles shows that the contact lines are pinned throughout stage 2, exhibiting only very minor local adjustments of the order of 0.05 mm-0.15mm for the duration of this stage of the process.Although refraction of the internally emitted fluorescent light as it leaves the drop may distort the perceived internal interface shape, we obtain a qualitative impression of the full three dimensional fluid movement during stage 2 by varying the lighting conditions.Figure 6(a), frames (a.1) through (a.5), shows the coalescence of drops when the spherical convex lens in Figure 1(b) is replaced with a concave lens, leading to a large laser light spot exciting fluorescein throughout the entire drop.In frame (a.6), the drop contact line is shown by illuminating the backlight lamp 1.The stratification is complete after 11 s, and a comparison of the upper images in frame (a.5), showing the extent of dyed, lighter fluid, and frame (a.6), showing the extent of the undyed, heavier fluid, demonstrates that the lighter fluid forms a layer occupying the top portion of the composite drop after stratification.It is evident that the dyed fluid flows in the y direction during the stratification process and that the coalescence induces three-dimensional flow.Under the modified lighting conditions, in frames (a.1) to (a.5), the observed intensity represents an integration of the light passing through the entire drop from the side and top.The bright line in the top view image, indicated by the arrows shown in frames (a.3) and (a.4), aligns with the position of the bottom contact line of the internal interface.This suggests that the heavier fluid has flowed under the lighter fluid but has not proceeded further to the left than the contact line in any cross section of the composite drop.
Figure 6(b) shows an image sequence of the coalescence process illuminated by the laser light sheet.In these images, the laser light sheet is scanned rapidly in the y direction, illuminating different three-dimensional cross sections of the drop during the internal flow of stage 2. The scanning occurs at ≈8 mm/s, a speed that is much faster than the speed of the fluid motion.Thus, the frames shown represent cross-sectional views of the fluid motion at essentially the same time point in the evolution of the stratification.The image sequences shown in Figure 6(a) (frames (a.3) and (a.4)) and Figure 6(b) taken together show that fluid in the central cross section moves at a faster velocity than the fluid occupying the cross sections closer to the drop edge.

Stage 3-Diffusive mixing
After the internal contents of the composite drop have become fully stratified at the end of stage 2, diffusion drives the final mixing of the initial fluid volumes.Figure 7 demonstrates the diffusive mixing of the fluids.The interface between the two fluids becomes visually blurred some time after stratification is complete.Line scans of intensity taken along a line perpendicular to the interface show that the composition gradient broadens diffusively.In a typical case, the width of the interface broadens about 450 µm over 50 s, where the width is defined as the distance between the locations along the line scan where the intensity is 10% and 90% of the maximum intensity.Using this broadening width and time scale, we estimate a diffusion coefficient of D ≈ 4 × 10 −6 cm 2 /s using the relationship D = l 2 broaden /2t D , where l broadening is the broadened width and t D is the time over which broadening occurs.Using this value of diffusion coefficient along with the measured velocity during the initial stages of motion of the gravity current, and the drop size as the characteristic length scale, results in an estimated Peclet number during stage 2 of Pe ≈ 2500.The large value of Peclet number suggests that diffusion has negligible impact on the transport during stage 2 before stratification is complete, consistent with our observations.

MODEL AND ANALYSIS
The experimental observations described above strongly suggest that fluid motion within the composite drop during stage 2 proceeds via a gravity current, promoting stratification of the initial fluid volumes into two layers.To verify the hypothesis that the fluid motion is a gravity current, we first conduct dimensional analysis to determine the scaling of the internal interface velocity with fluid properties and flow conditions.We assume that the hydrostatic pressure variation driving flow within the composite drop scales as ∆ρgH, and that the viscous stress resisting motion scales as µ b L/T H, where T is the characteristic time (use of the viscosity of the heavier fluid is suggested by Huppert 46 ).Equating the driving pressure gradient and the viscous stress, we obtain a characteristic time T ≈ µ b L/∆ρgH 2 .Noting that inertia has decayed by the end of stage 1, and therefore, the Reynolds number is small during stage 2, and also noting that the composite drops always exhibit aspect ratios H/L less than unity, we anticipate a time scaling typical of those used in lubrication analyses of thin film geometries.Thus, we modify the characteristic time scale by a factor of the aspect ratio, such that To verify the time scale indicated by dimensional analysis, we examine experimental data from various fluid pairs, isolating changes in density difference, drop geometry, and viscosities.Figure 8 compares the displacement of the top point on the internal interface, x f , for two fluid pairs with similar viscosities and sizes but different density differences ∆ρ.displacement as a function of time.The larger density difference drives faster displacement of the interface, leading to a more rapid approach to the final stratified interface position.Figure 8 shows the displacement of the internal interface, scaled by the length of the composite drop, as a function of the time scaled by the characteristic time T c given in Eq. ( 1).On normalized axes, the experimental curves collapse, suggesting that the data scale with density difference in a manner consistent with that expected for a gravity current.Figure 9 compares several displacement curves for a single fluid pair, varying the initial drop volumes and therefore varying the composite drop dimensions H and L. Figure 9(a) shows the displacement of the top point x f on the internal interface as a function of time for different drop volumes.While the approach to the final displacement occurs at approximately the same initial rate for all of the cases shown, the final displacement is larger and the time at which final displacement is reached is later for larger volumes.The final displacement of the internal interface depends on drop volume because the less dense fluid occupies a different layer thickness within the composite drop at the end of stratification.The final position of the interface depends on the three-dimensional composite drop shape, and thus the Bond number, as well as the drop volume.These factors render it difficult to determine a well-defined analytical expression for the final interface displacement.Figure 9(b) shows the normalized interface displacement scaled by the composite drop length L, as a function of time normalized by the characteristic time scale given in Eq. ( 1).The dimensionless curves of Figure 9(b) collapse to nearly a single curve with some spread in the shoulder region, indicating that the geometric dependence of the time scale given above works well, including the extra factor of the aspect ratio H/L.However, the value in the plateau region varies monotonically with drop volume or Bond number.The normalization with L accounts for the major effect of drop volume but does not account for the flattening of the drop with increasing Bond number.This residual effect of Bond number causes the observed spread in plateau values.
Figure 10 considers a wide variety of fluid pairs with varying density difference and fluid viscosities, keeping volume and thus the Bond number nearly fixed.Figure 10(a) shows that the displacement of the top point on the internal interface, x f , approaches a similar final value, as a result of the similar composite drop volumes.However, the rates of approach to the final displacement values vary considerably for the eight fluid pairs considered.While the data sets are not clearly ranked in terms of any individual fluid property, it is notable that the data set corresponding to the slowest approach to equilibrium corresponds to fluid pair BJ, which has the smallest density difference of the group.The data set with the fastest approach to equilibrium corresponds to fluid pair AH, which has the largest density difference.These observations are consistent with expectations for the gravity current scaling given by Eq. ( 1).For the data sets exhibiting intermediate time scales, it is clear that both fluid viscosities also play a role in the rate of approach to equilibrium.
Figure 10(b) shows the result of normalizing the displacement with the maximum displacement of the stratified interface, plotted as a function of time normalized by the characteristic time scale given by Eq. ( 1).The data sets collapse considerably to occupy nearly a single curve, indicating that the motion of the internal interface is well described by the expected scaling for a simple gravity current.However, the collapse of the data is not as good as in Figures 8 and 9, which considered only density difference and composite drop geometry.Specifically, examining the transition region from a rapidly advancing gravity current to complete stratification of the two layers, we see that the data sets deviate monotonically with viscosity ratio.The leftmost data set showing the fastest approach to equilibrium corresponds to fluid pair AB with viscosity ratio λ = 0.16, while the rightmost data set showing the slowest approach to equilibrium corresponds to fluid pair AH with viscosity ratio λ = 10.The monotonic variation with viscosity ratio is notable, since the viscosity of the heavier fluid in most of the fluid pairs considered is similar, falling between µ b ≈ 4 and 8 mPa s for all cases except fluid pair BA, in which µ b = 41 mPa s.Unlike in previous analyses of gravity currents, in which the characteristic time scale depends only on the heavier fluid viscosity, 46 it appears that the stratification of two fluid layers in a merging sessile drop is also influenced by the viscosity of the lighter fluid layer, which can also play a role in retarding the flow.
Motivated in part by the lack of complete collapse of the displacement curves with the scaling corresponding to a simple gravity current shown in Figure 10(b), we seek a more detailed description of the internal fluid flow within the composite sessile drop during stage 2 using the lubrication approximation to simplify the governing equations.The lubrication approximation has been previously applied to model similar phenomena including the internal interface motion for two layer systems 54,55 and a gravity current in a tube geometry. 48,49,56The two main assumptions of the lubrication approximation are (1) the reduced Reynolds number is small, Re(H/L) << 1, and (2) the square of geometric aspect ratio is small, (H/L) 2 << 1. 57 In our experiments, the aspect ratio of the composite drop H/L ranges from 0.2 to 0.4.The Reynolds number is maximum at the beginning of stage 2 and ranges from 0.1 to 1.5 for different experimental fluids considered here.Thus, the two conditions are satisfied: (1) the reduced Reynolds number varies from 0.03 < Re(H/L) < 0.6 and (2) the square of the aspect ratio varies from 0.04 < (H/L) 2 < 0.16.Experimentally, we observe that the external and internal interfaces meet the solid substrate at a high angle, potentially violating the lubrication approximation.The external interface contact angle in our experiments is about 90 • , suggesting that the wedge flow near the contact line could impact the validity of the lubrication approximation when the internal interface approaches the corner region.To mitigate these effects, we examine the interface evolution only at early times, when the internal interface is far away from the corner and the impact of the wedge flow is negligible.Even though the internal interface itself meets the solid substrate at a high angle, the lubrication approximation is still expected to yield good qualitative results.In cases similar to ours, [58][59][60][61] it has been shown that a violation of the lubrication approximation near the contact line does not contaminate the interface shape and propagation speed of the overall gravity current.In the Simulation section, we use the numerical results to explicitly assess the validity of the lubrication approximation at the internal interface in our case.
Based on experimental observations, we fix the contact line position and assume that the shape of the external interface of the composite drop does not vary with time.Because the fluids are miscible aqueous solutions and the measured air/liquid surface tensions are very similar, we assume the interfacial tension of the internal interface is zero.The differential equation describing the time evolution of the internal interface profile is given in Eq. ( 2) in dimensionless form.The derivation of Eq. ( 2) can be found in Appendix A. ( The function b(x, y,t) describes the internal interface profile, h(x, y,t) describes the external interface profile (see Figure 5), ρ t is the density of the lighter fluid, ∆ρ is the density difference between the two fluids, γ is the surface tension of the lighter fluid, and µ b is the viscosity of the heavier fluid.The x coordinate is scaled by the composite drop length L while the y coordinate is scaled by the drop width L 2 , and the z coordinate is scaled by the height of the composite drop H.The characteristic time scale that naturally arises from this scaling is given by T c = 3µ b L 2 ∆ρg H 3 .In addition to the evolution equation given by Eq. ( 2), mass conservation of the fluids in the top and bottom components of the composite drop must be invoked (see Figure 5(a) for the definition of the relevant domain), and where m t and m b are the masses of the lighter and heavier fluids, respectively, and dxd y and V 2 =  Equation (2) shows that the flow behavior in the x and y direction is similar if we permute the role of x and y in the right hand side of the equation.For the flow in the x direction, the first term on the right hand side ∂ ∂x (b 3 ∂b ∂x ) corresponds to the gravity current and the second term represents the coupling of the evolution of the external interface, h, to the internal interface, b.If the external interface is in static equilibrium, then the second term is identically zero, as shown by taking the first order derivative of the Young-Laplace equation for the equilibrium drop shape: In cases where the composite drop shape deviates from equilibrium, then depending on the extent of the deviation from the equilibrium shape and the value of ρ t /∆ρ, the coupling term can have a significant impact on the evolution of the internal interface.Our observation that the external interface does not change during stage 2 beyond our detection limits for almost all the fluid pairs suggests that the second term is negligible in these experiments.However, as the density difference becomes as large as 20%, we observe more significant changes in the external interface, which may influence the internal fluid motion.In later comparisons between experimental data and analysis, we neglect this small change.
Equation (2) shows that there are several relevant dimensionless groups, including a fractional density difference ρ t ∆ρ , a Bond number γ ρ t g L 2 , and an aspect ratio L L 2 (the ratio of length by width as shown in Figure 2).The time scale for the evolution of the internal interface, b(x, y,t), does not explicitly depend on the viscosity of the upper, lighter fluid in contrast to the experiments shown in Figure 10(b).However, the evolution of the external interface shape h(x, y,t) depends on both the heavier fluid viscosity and the lighter fluid viscosity, thus causing the evolution of the internal interface to implicitly depend on both viscosities through the coupling given in Eq. ( 2).Thus, the simplification that the external interface is in static equilibrium effectively removes any dependence on the lighter fluid viscosity.This suggests that the lack of full collapse shown in Figure 10(b) is a direct result of motion of the external interface, even though this motion may be smaller than our ability to detect it.We note, in particular, that the fluid pairs AH and AI shown in Figure 10 have the largest density differences.These cases exhibited measurable changes in the external interface shape, leading to a more prominent coupling between the internal interface b(x, y,t) and the external interface h(x, y,t) and the most prominent deviation from a simple gravity current.

SIMULATION
Since the drops we examine are neither axisymmetric nor planar, modeling the complete situation would require solving Eq. (2) for a fully three dimensional flow.Simulations of free-surface flows with fully three-dimensional flow are computationally intensive, 24 and so, we chose to solve the two-dimensional version of Eq. ( 2) for ease of simulation.A two dimensional simulation still provides qualitative information that can be compared with the experimental observations of the evolution of the displacement of the internal interface.The two-dimensional assumption simplifies the problem by reducing the aspect ratio for the characteristic lengths in the x and y directions to zero, as well as removing the complication of coupled flow in the x, y, and z directions due to mass conservation.The resulting governing equation is We assume a constant external composite drop shape h(x,t) = heqm (x), consistent with the experimental observations of Figure 4.This simplification decouples the mutual dependence of the internal interface b(x,t) and external interface h(x,t).We assume an equilibrium shape for the external interface given by the two-dimensional Young-Laplace equation subject to the lubrication approximation, γ heqm ′′ = ρg heqm .After scaling with the height and length of the composite drop such that heqm ∼ H h eqm , x ∼ Lx, the dimensionless equilibrium drop shape is where Bo = ρ t gL 2 /γ defines the Bond number for the composite drop.In the simulation, choosing the equilibrium shape for the external interface means the second term on right hand side of Eq. (2 ′ ) is zero and the internal interface equation then reduces to (Figure 3), we choose a steeply sloped linear profile given by b = 19.964x+ 0.499, which is the line connecting the two points (−0.025, 0) and (0.025, h eqm (0.025)).We use a slightly tilted line since a vertical line presents difficulties in taking derivatives of the interface profile.The initial interface for the simulation is shown by the dashed line in Fig. 11.Since Eq. (2 ′′ ) is a second order nonlinear diffusion equation, in addition to the initial condition b(x, 0) as specified above, we set the boundary conditions as b(x f ,t) = h eqm (x f ) and b(x b ,t) = 0 for all time.The coordinates x f (t) and x b (t) are defined as the contact positions of the internal and external interfaces, respectively, as shown in Figure 5(a).The specification of initial and boundary conditions along with the mass conservation condition, Eq. ( 3), guarantees the existence and uniqueness of the solution of Eq. (2 ′′ ). 61,62The positions x b (t) and x f (t) are determined by this unique solution to the evolution equation.The position x b (t) is taken to be the mesh point at which the interface intersects the solid substrate.The position x f (t) is taken to be the mesh point at which the interface intersects the external interface, adjusted to conserve the volume  x f (t) x b (t) b(x, y,t)dxdy +  1/2 x f (t) h(x,t)dxd y = V b of the heavier drop.For the simulation, the conservation of mass conditions that depend on h eqm (see Eq. ( 4)) become requirements for area conservation for each fluid component in two dimensions.As the differential equation will be discretized by a co-volume method as discussed later, this area conservation constraint is satisfied automatically.The interface height at the left boundary point b(x b ,t) is set to zero as the interface propagates, whereas the right boundary point b(x f ,t) is forced to be on the external boundary h eqm at discretized mesh points.This results in an error in area conservation of ∆x 2 /2, where ∆x is the mesh size for the discretization.
The differential equation Eq. (2 ′′ ) is a type of degenerate parabolic problem which can be solved numerically using a co-volume method suggested by Baughman and Walkington. 63We use the explicit finite difference scheme to discretize the interface profile b(x,t) as b n j , where j denotes the spatial index and n denotes the temporal index.The temporal derivative of b is given by in the explicit scheme, where ∆t is the time step.The second order derivative for b 4 can be written symmetrically using the scheme analyzed by Baughman and Walkington. 63 For this type of second order nonlinear diffusion equation, the stability criterion for the numerical scheme is given by an expression of the form ∆t ∼ (∆x) 2 . 63,64In addition to this stability criterion, b n j must be non-negative to be in the physical subspace.Since the value of b n j is non-negative from the initial condition in our simulation, this requirement is satisfied.Thus, Eq. ( 5) is solved by iteration of time steps to obtain the time evolution of the interface profile b using the above explicit scheme.In our simulation, we calculate the evolution of the interface from time t = 0 to t = 8.The time steps are determined using the stability criterion, once the spatial interval ∆x is selected using the relationship ∆x = 1 n x , where n x is the number of grid points.Figure 11 shows the evolution from t = 0 to t = 8 using a uniform grid ∆x = 1 1280 = 7.81 × 10 −4 and time step ∆t = 1.22 × 10 −7 with equilibrium external interface shape of h eqm .For times beyond t = 8, the left hand side of the internal interface collides with the external interface, where the recirculation in the wedge must be considered and Eq. ( 5) is rendered invalid, thus we terminate the simulation thereafter.Convergence of the numerical solutions is examined in detail in Appendix B.
As shown in Figure 11, the internal interface has a high slope at the contact point x b which brings the lubrication approximation into question.For the lubrication approximation to be valid, the dimensional terms Re(u ∂u ∂x ), Re(w ∂u ∂z ), and ∂ 2 u ∂x 2 must be small compared with ∂ 2 u ∂z 2 , where u and w are the horizontal and vertical components of the velocity, respectively.From the numerical solution of the evolution equation, Eq. (2 ′′ ), we calculate the velocities and relevant derivatives at an early time in the evolution of the gravity current.These comparisons show that the ratio of Re(u ∂u ∂x ) to ∂ 2 u ∂z 2 is smaller than 0.001 over the entire domain.The ratio of Re(w ∂u ∂z ) to ∂ 2 u ∂z 2 is smaller than 0.6 over the entire domain and declines quickly to smaller than 0.1 within 0.0006 spatial units (5 spatial grid points out of 6400) from x b .Unlike the other two ratios, the ratio of ∂ 2 u ∂x 2 to ∂ 2 u ∂z 2 is large at the contact point x b due to the large slope; however, the ratio of these two terms declines quickly to smaller than 0.1 within 0.002 spatial units (15 spatial grid points out of 6400) from x b .Therefore, the lubrication assumption fails very near the corner and is satisfied away from the contact point x b .The governing equation of motion Eq. (2 ′′ ) is of the same form as that derived in the gravity current of Huppert, with different boundary conditions. 61The solution obtained by Huppert also has a high slope at the front of the current, and the lubrication approximation is also violated near the moving contact line.Nevertheless, the lubrication solution describes the experimentally observed shape of the gravity current well and it describes the propagation speed of the gravity front very well.Gratton and Minotti 65 and Betelú et al. 66 address in more detail the agreement of the gravity current profiles with experiments even though the lubrication approximation fails at the corner.The prior results and our examination of our numerical results justify the use of the lubrication approximation to qualitatively describe the experimental observations of the evolution of the shape and propagation speed of the internal interface.12, the simulation qualitatively captures the displacement of the internal interface both at the top intersection with the external interface and at the bottom intersection with the solid surface.In both the experiment and simulation, the displacement of the point intersecting the external interface advances faster than that of the point intersecting the rigid substrate as would be expected from the contact line motion at a liquid or solid surface.In both the simulation and experiment, the intersection with the substrate moves farther than the intersection with the top interface before the final equilibrium of full stratification.
Since the external interface shape sets the initial and also boundary condition for the internal interface evolution, the gravity flattening of the external interface shape due to the Bond number variation (see Eq. ( 4)) will also affect the internal interface motion.Figure 13 shows the simulation examining the Bond number dependence of the top point of the internal interface.As seen in experiments in Figure 9(b), the simulation shows that the interface point evolves with almost the same initial rate while the displacement within the plateau region increases monotonically with respect to the Bond number.
Although the simulation results qualitatively agree with the experiments, the dimensional time scale in the simulation is faster by about an order of magnitude.This discrepancy may arise from the quantitative comparison of the lubrication solution with experiments for a highly sloped interface 58,59 but more likely arises because the simulation is two dimensional, while the experiments show three dimensional flow as seen in Figure 5. Similar to our results, in other studies of gravity driven currents, the observed characteristic velocity in three dimensional experiments is always smaller than that predicted by a two dimensional analysis because of the outer geometry and boundary conditions. 56,67,68Despite the discrepancy in the magnitudes of the time scales, the scaling of the characteristic time scales for the gravity current with fluid properties and drop geometry agrees well with the observed behavior, suggesting that the lubrication model captures most of the essential physics in the problem.

CONCLUSIONS
We have examined the coalescence of two sessile drops of miscible fluids with similar surface tensions and equal volumes.The fluid pairs can have different viscosities and densities, and the overall drop volumes vary.The fluid systems examined exhibit surface tensions of the order of 60 to 70 mN/m, viscosities from about 4 to 41 mPa s, density differences from 0.01 to 0.43 g/cm 3 , Bond numbers from 1 to 8, and contact angles of about 90 • .The substrates consist of cross-linked PDMS, which is known to form a ridge beneath the contact line.The ridges do not completely pin contact lines but do inhibit their motion.
The coalescence of the two initial sessile drops occurs in three stages of motion: (1) bridge healing, (2) a gravity current, and (3) long-term diffusive mixing.Stage 1 is governed by minimization of capillary energy with coincident movement of the contact lines.The bridge healing of stage 1 ends with a smooth external interface of the composite drop and no apparent twisting of streamlines that would cause advective mixing.The external interface achieves a nearly static shape within the resolution of our measurement for the remainder of the coalescence process, for density differences smaller than 0.2 g/ml.The time scale for stage 1 is well separated from that of stage 2, and thus, the end of this stage forms an initial condition for the subsequent stages.Stage 2 is controlled by gravity and viscosity even for density differences as small as 1%.The fluids move toward a simple stratified configuration with no advective mixing and little to no diffusive mixing, consistent with the estimated Peclet number Pe ≈ 2500.A three-dimensional lubrication analysis indicates that the predominant time scale for stage 2 is T c = 3µ b L 2 ∆ρg H 3 as long as the external interface is very near equilibrium.Our analysis suggests that even small deviations of the external interface from equilibrium may cause the flow to have significant contributions from capillarity as well.The flow in stage 2 is also governed by the additional dimensionless groups ρ ∆ρ and γ ρ t g L 2 .The experiments are in good agreement with the time scaling predicted by the analysis.The presence of small motions of the external interface leads to additional coupling with the flow of the upper, lighter fluid layer, and resistance due to the viscosity of this layer contributes to the motion of the internal interface.Thus, a systematic deviation of the scaled displacement curves is observed to monotonically vary with the viscosity ratio.These conclusions are also supported by numerical solutions of a two-dimensional version of the equations developed in the analysis.However, the two-dimensional numerical solutions overestimate the characteristic speeds of the flow.This is likely due to the fact that the flow is in reality three-dimensional as shown by experiments.It is also possible that the high contact angle results in an aspect ratio of the composite drop that is too large to be accurately treated by the lubrication approximation.
Stage 2 concludes with the complete stratification of the two fluids after time scales of the order of tens of seconds.Finally in stage 3, diffusive mixing occurs in the composite drop over time scales of the order of several minutes.
Experiments and simulations both show that once the initial bridge healing is complete, internal motion proceeds at low Reynolds number over very long time scales, even when the external interface appears to have long reached a static equilibrium.The fluid motion follows a well-defined gravity current even for very small density differences.These results show that imaging of only the external interface shape of merging sessile drops is not sufficient to ensure that mixing of two fluid volumes is complete, and that such mixing often occurs only over much longer time scales than expected.
The temporal evolution of the internal interface can be obtained by combining Eqs.(A5b) and (A6b), which is a kinematic boundary condition, to give Using the velocity profiles obtained from Eqs. (A5a)-(A6b), we obtain the differential equation Nondimensionalizing Eq. (A10a) using b = H b, h = H h, x = L x, y = L 2 ỹ, z = H z, where H, L, L 2 are shown in Figure 2 yields the following equation after dropping the tilde for the dimensionless variables: The necessity of a second length scale L 2 is discussed by Szeri 57 in deriving three dimensional lubrication theory.This equation suggests a time scaling of T c = 3µ b L 2 ∆ρg H 3 for the internal interface motion, which is identical to Eq. (1) except for the factor of three.The resulting dimensionless governing equation for the interface motion is given by Eq. ( 2), (2)

APPENDIX B: EXAMINATION OF CONVERGENCE FOR THE SIMULATION
To examine convergence, we first maintain a spatial mesh size of ∆x = 1 1280 = 7.81 × 10 −4 and vary the temporal step from ∆t = 9.77 × 10 −7 to ∆t = 6.10 × 10 −8 to check the error dependence on time steps. 74As this nonlinear diffusion equation does not have an analytical solution, we define the global error comparing the interface profile at t = 1 with the finest case where ∆x = 7.81 × 10 −4 and ∆t = 3.05 × 10 −8 using the square norm  thus the errors at the corners are the leading error in the numerical simulation. 63Baughman and Walkington 63 proved that for equations of this form, the estimated error approaches zero at least as rapidly as C(∆x 1/2 + ∆t 1/2 ), where C is a constant depending on the structure of the equation and boundary conditions.Figure 14(a) shows that the global error for our simulation decreases with first order behavior in the time step, which agrees with the convergence criteria.
Similarly, to examine the convergence with grid size, we calculate the global error for cases in which ∆t = 6.10 × 10 −8 is held fixed, but the grid varies from ∆x = 6.25 ) 2 ] 1/2 , (B2) where k = 5120/n x allows us to restrict the difference calculation at coarser grids to occur at the same x coordinate for different grid sizes.Figure 14(b) shows that the error of the simulation also decreases with nearly first order dependence on grid size.Convergence with spatial discretization is indicated for the numerical scheme with the same reasoning as for temporal discretization.
FIG. 1. Experimental setup.(a) Schematic diagram of droplet feed system.(b) Schematic diagram of optical train.

022101- 6 Zhang
FIG. 3. Time sequence of images of the coalescence process for two experiments with different fluid pairs.(a) Fluid pair BB for initial drop volumes of 20.2 µl.Both the fluorescent dyed fluid and the undyed fluid have a density of 1.145 g/ml.(b) Fluid pair CB for initial drop volumes of 22.7 µl.The fluorescent dyed fluid has a density of 1.071 g/ml while the undyed fluid has a greater density of 1.145 g/ml.
FIG. 4. (a) Overlay of edge profiles from top view images showing the position of the contact lines at two different times: just after the end of stage 1 (t = 0.13 s) and toward the end of stage 2 (t = 3 s) for fluid pair CB at a volume of 22.7 µl.The difference in contact line locations for these selected times is small and obscured by image resolution.(b) Overlay of two side view images at 0.13 s and 3 s.

FIG. 5 .
FIG. 5. (a) Schematic diagram of the internal interface within the composite drop.The position x = 0 is where the vertical internal interface is located at the end of stage 1, x b is the internal interface point in contact with the solid surface at time t, x f is the external interface point where the two fluids meet at the top of the drop at time t.(b) Measured x f as a function of time.The vertical line indicated on the plot at right corresponds to t = 0.1 s and represents the beginning of the gravity-driven flow of stage 2 of the coalescence process (fluid pair CB; volume 22.7 µl).For clarity, the error bars are not shown since the uncertainties for values of x f at different times are identically 0.045 mm due to the camera resolution, which falls within the symbols.

FIG. 7 .
FIG. 7. (a) Schematic diagram depicting the location of the line scan used to estimate the broadening of the interface.(b) Line scans of intensity across the internal interface for fluid pair CB with volume 8 µl, taken at t = 7 s and 57 s after the onset of coalescence.

Figure 8
FIG. 8. Displacement of the internal fluid interface as a function of time for varying density differences and similar composite drop size and viscosity.The inset shows unscaled displacement versus time for fluid pairs BJ and BK.The main figure shows the displacement versus time scaled with the characteristic time scale given in Eq. (1).

FIG. 10 .
FIG. 10.(a) Displacement versus time for several fluid pairs with different viscosities and density differences, keeping drop volume nearly constant.(b) Displacement normalized by the maximum extent of the internal interface versus time normalized by the characteristic time given in Eq. (1).Inset shows an expanded view of the transition region from a rapid gravity current to fully stratified layers.

S 3 h
(x, y,t)dxd y indicate the volumes occupied by the heavier fluid under the internal interface in the central portion of the composite drop and under the external interface in the end portion of the drop where the fluids do not overlap, respectively.
FIG. 11.Evolution of the dimensionless internal interface profile for a two-dimensional composite drop assuming an equilibrium composite drop shape.Internal interface profiles are shown for dimensionless time steps from 0 to 8 time units with 1 unit intervals.

2 .
The final discretized expression for Eq.(2 ′′ ) is then given by

FIG. 12
FIG. 12. (a) Measured absolute displacement of the two interface points in the composite drop.(b) Simulation tracking the intersections of the internal interface with the rigid substrate and the external interface.Experiments and simulations both correspond to fluid pair DG with λ = 1, Bo = 8.The time scale for these two graphs differs by a factor of 10.

Figure 12 (
Figure 12(b) shows a typical result of a converged simulation for values of fluid properties taken from the experiments.The displacement of the interface is qualitatively similar to that observed in experiments as shown in Figures 3 and 4(b).The simulation only tracks the gravity current until the bottom point on the interface reaches the intersection of the external interface and the solid surface (the contact line).As shown in Figure12, the simulation qualitatively captures the displacement of the internal interface both at the top intersection with the external interface and at the bottom intersection with the solid surface.In both the experiment and simulation, the displacement of the point intersecting the external interface advances faster than that of the point intersecting the rigid substrate as would be expected from the contact line motion at a liquid or solid surface.In both the simulation and experiment, the intersection with the substrate moves farther than the intersection with the top interface before the final equilibrium of full stratification.Since the external interface shape sets the initial and also boundary condition for the internal interface evolution, the gravity flattening of the external interface shape due to the Bond number variation (see Eq. (4)) will also affect the internal interface motion.Figure13shows the simulation examining the Bond number dependence of the top point of the internal interface.As seen in experiments in Figure9(b), the simulation shows that the interface point evolves with almost the same initial rate while the displacement within the plateau region increases monotonically with respect to the Bond number.Although the simulation results qualitatively agree with the experiments, the dimensional time scale in the simulation is faster by about an order of magnitude.This discrepancy may arise from the quantitative comparison of the lubrication solution with experiments for a highly sloped interface58,59 but more likely arises because the simulation is two dimensional, while the experiments show three dimensional flow as seen in Figure5.Similar to our results, in other studies of gravity driven currents, the observed characteristic velocity in three dimensional experiments is always smaller than that predicted by a two dimensional analysis because of the outer geometry and boundary conditions.56,67,68Despite the discrepancy in the magnitudes of the time scales, the scaling of the characteristic time scales for the gravity current with fluid properties and drop geometry agrees well with the observed behavior, suggesting that the lubrication model captures most of the essential physics in the problem.
temporal evolution of the external interface obtained from application of the material derivative is given by j | ∆t − b t=1 j | ∆t=3.05×10−8) 2 ] 1/2 , (B1)where b t=1 j | ∆t is the calculated interface height with the time step ∆t, and b t=1 j | ∆t=3.05×10−8 represents the interface height with the finest time step.Due to the finite propagation speed of the nonlinear diffusion equation, functions at corners are not as smooth compared to the linear diffusion equation,

FIG. 14
FIG. 14.(a) Global error as a function of time step (solid triangles represent the error defined by Eq. (B1) at different time steps).(b) Global error as a function of spatial discretization (solid triangles represent the error defined by Eq. (B2) at various grid sizes).

TABLE I .
Fluid properties of glycerol solutions.