# Self-assembly and Particle Aggregation in Stratified Fluids

#### Experimental Discovery and First-principle Mechanisms of Diffusion-induced Flows in Exterior Domains

The motion of particulates in the ocean as they fall under gravitation constitutes a critical component of the carbon cycle. Phytoplankton can photosynthesize dissolved carbon dioxide gas into dense, solid carbon debris called “marine snow,” which facilitates the transport of carbon from the surface waters to the ocean depths. Over the course of eons, the carbon eventually turns into oil. The prevailing view maintains that particulate aggregation occurs through random turbulent oceanic motions that generate collisions between small-scale sticky particulates, which then form aggregates [4].

However, we recently discovered that such particles—when suspended at similar depths in stratified water—can slowly form aggregates in the absence of both external random flows and adhesion [7]. Surprisingly, these particles experience an effective horizontal attractive force that originates from the no-flux boundary condition. This attraction leads to self-assembly, during which particles form molecule-like assemblages that subsequently appear to dynamically solve jigsaw-like puzzles while en route to large-scale particulate discs.

Figure 1a depicts a time-series montage of the formation of a “particulate molecule” that eventually merges with a larger aggregate, and Figure 1b portrays a schematic of the experiment. Figures 1c and 1d present an initial image of the particulate distribution and a visual from 250 minutes later. Our findings offer insights into new mechanisms for large-scale aggregate formation in a wide variety of physical systems, including the ubiquitous stratified ocean waters, potential microplastic organization mechanisms, magnetic dynamos on Mercury and Ganymede, and even origin-of-life questions (such as matter’s self-organization in deep ocean vents).

**Figure 1.**Experimental self-assembly.

**1a.**Time-series montage that illustrates the formation of a molecule. Five different output times—with 10 minutes between each instance—appear from left to right (artistic blurring depicts the motion). The spheres’ radii and densities are respectively \(0.025\)-\(0.05\) cm and \(1.05\) g cc

^{−1}; the top fluid is fresh water of density \(0.997\) g cc

^{−1}and the bottom is a sodium chloride water solution of density \(1.1\) g cc

^{−1}, with a transition thickness of approximately \(2\) cm.

**1b.**Schematic of the experimental setup.

**1c.**Initial cluster.

**1d.**Final self-assembled cluster. Figure adapted from [7].

The discovery itself arose through a fortunate series of laboratory mistakes that ultimately led to a “eureka” moment. During a routine outreach demonstration, a setup error revealed this self-assembly phenomenon. Our demo—which seeks to exemplify stratification’s ability to trap particulates within the transition layer between top and bottom densities for a short time period—drops hundreds of thousands of glass spheres (with radii of 50 microns) into a stratified tank. The particles form a trapped cloud within the transition layer, and each particle is coated in the buoyant upper fluid; the collective cloud suspends for a transient period until the particles diffusively shed their lightweight upper fluid coatings. The system then begins to “rain” through a Rayleigh-Taylor-like multiphase instability, and all suspended particulates eventually reach the tank bottom. This phenomenon is closely related to our previous discovery that a dense sphere can actually bounce off of an internal density transition layer due to related mechanisms [1], which led to experimental, computational, and analytical studies of single body behavior [5, 6].

During one demonstration of this raining phenomenon, we made the bottom layer too salty and the particulate was buoyant with respect to the lower fluid. As such, no rain-out phase was evident after 10 minutes of observation. We checked the densities with our densitometers and confirmed the error. But rather than immediately clean up the mess, we left the system overnight; our procrastination was incredibly fortunate because the extremely slow dynamics had time to play out. The next morning, the diffuse particulate system had self-assembled into a large-scale aggregate disc. We knew that we had stumbled onto an intriguing observation and consequently launched an intense experimental, theoretical, and computational campaign over the next several years.

We suspected that the underlying mechanism pertained to a so-called *diffusion-induced flow*. Independent of one another, Owen Phillips and Carl Wunsch analytically discovered such flows for the simple case of a tilted infinite plane in a linearly stratified fluid [8, 9]. Their works derived mathematically exact, steady solutions of the incompressible Navier-Stokes equations coupled to a solute field that is governed by an advection-diffusion equation. The initially horizontal isolines of salt concentration do not orthogonally intersect the tilted plane, as required by the no-flux boundary condition. This situation causes a turning of the isolines, which then creates a buoyancy mismatch that drives a flow up the inclined plane along a narrow boundary layer near the surface.

**Figure 2.**Comparison of the experiment and simulation for the flow field in a vertical plane that slices the spheroid through its north and south poles. Here, \(Pe=52\), \(\mu=0.016\) poise, \(\kappa=1.5\times10^{-5}\) cm

^{2}s

^{−1}, the spheroid vertical radius \(= 0.5\) cm, the width \(= 2\) cm, and the density gradient \(\sigma=0.002\) g cm

^{−4}.

**2a.**Flow speed from particle imaging velocimetry.

**2b.**Flow speed from COMSOL.

**2c.**Horizontal flow velocity from particle imaging velocimetry.

**2d.**Horizontal flow velocity from COMSOL. Figure courtesy of [7].

The Phillips/Wunsch solution reveals some interesting scales for diffusion-induced flows. The characteristic velocity and length scales for the Phillips solution are \(U = \kappa \left( \frac{g\sigma}{\kappa \mu} \right)^{1/4}\) and \(L = \left( \frac{\kappa\mu}{g \sigma} \right)^{1/4}\), where \(g\) is the gravitational acceleration, \(\mu\) is the dynamic fluid viscosity, \(\kappa\) is the salt diffusivity, and \(\sigma\) is the slope of the background density field (all assumed to be constant). The length scale \(L\) defines the boundary layer’s thickness above the tilted flat plane, and the velocity scale \(U\) defines its strength. We use this velocity scale to nondimensionalize the equations of motion for the velocity field \(\boldsymbol{u}\), pressure \(P\), density \(\rho\) (which varies only during the evolution of the salt concentration), position \(\boldsymbol{x}\), and time \(t\) (for a single sphere of radius \(a\)) via \(\tilde{\boldsymbol {x}} =\frac{\boldsymbol{x}}{a}\), \(\tilde{\boldsymbol {u}} = \frac{\boldsymbol{u}}{U}\), \(\tilde{\rho} =\frac{\rho}{\sigma a}\), \(\tilde{P} = \frac{Pa}{\mu U}\), and \(\tilde{t} = \frac{t \kappa}{a^2}\). The resulting nondimensional partial differential equation (PDE) system (dropping tildes and primes) for the incompressible fluid velocity and concentration field is respectively

\[Re \ \rho \left[ \frac{1}{Pe} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \nabla \boldsymbol{u} \right] = - \nabla P - Pe^3 \rho \hat{z} + \Delta \boldsymbol{u}\tag1\]

\[\frac{\partial \rho}{\partial t} + Pe \ \boldsymbol{u} \cdot \nabla \rho = \Delta \rho.\tag2\]

Here, the Reynolds number is \(Re = \frac{\sigma a^2 U}{\mu}\) and the Péclet number is \(Pe = \frac{a}{L}\). For all of our experiments, \(Re<0.001\), so we employed the so-called Stokes approximation to set the left side of \((1)\) to zero while retaining nonlinear effects through nonzero \(Pe\) in the advection-diffusion of the salt concentration. As such, the system’s only remaining parameter is the Péclet number. We assume that the boundary conditions are no-slip for the velocity and no-flux for the tracer on any physical boundary.

Given our experimental observations of diffusion-induced flows that self-assemble into large-scale compact discs, we first examined the case of a stationary single oblate spheroid in a linearly stratified fluid. We used particle tracking velocimetry to measure the exterior flows to the compact domain, and the finite element package COMSOL generated solutions of the coupled PDEs. Figure 2 provides a strong quantified comparison of the experimental and computational flows. The plots reveal an attractive flow into the body’s equator, which is spun up and down and ejected at the north pole as a stronger horizontal flow that fans radially outward. This equatorial flow drives the self-assembly in Figure 1, which causes nearby particles to coalesce into a larger disc.

Next, we considered the behavior of the two-body problem. We first studied some three-dimensional COMSOL simulations to compute the diffusion-induced flow that is induced by two fixed spheres in steady state [7]. Figure 3a illustrates the velocity distribution in the symmetry plane that cuts vertically through the north and south poles of both spheres. We computed the viscous stress integral around one of the spheres to determine the horizontal force that the diffusion-induced flow exerts on the sphere for a range of fixed separation distances (see Figure 3b). We fit the data to a power law force \(\frac{F}{F_0}=\left(\frac{d}{a} + 0.3249\right)^{-1.466}\) with a prefactor \(F_0=2.0865 \times 10^{-6}\) dynes, where \(d\) is the separation distance and \(a\) is the particle radii. The size of prefactor \(F_0\) is reasonable; balancing Stokes drag with this force yields particle velocities on the order of 2.5 microns per second, which is roughly the average speed that we observed experimentally.

**Figure 3.**COMSOL force calculations for three-dimensional simulations of two identical spheres of radius \(a\) that are separated by an amount \(d\) in linear stratification at \(Pe=5.35\).

**3a.**Steady-state density and flow field for two spheres at \( d=0.045\) cm in a vertical plane that slices the north and south poles.

**3b.**Steady attractive force on one sphere, which we compute by numerically integrating the projected stress tensor. The physical parameters are \(a = 0.045\) cm, \(\kappa=1.5\times10^{-5}\) cm

^{2}/s, \(\sigma = 0.03\) g/cm

^{4}, \(g = 981\) cm/s

^{2}, and \(\mu=0.0113\) poise.

**3c.**Comparison of Stokesian dynamics and two-sphere experiments for two linear background density fields: \(\sigma = 0.03\) g/cm

^{4}and \(\sigma=0.025\) g/cm

^{4}. Figure courtesy of [7].

We then modified a well-tested Stokesian dynamics numerical algorithm [7]—based on the theory for multibody low Reynolds number dynamics [3]—by adding this computationally measured force. Figure 3c depicts the simulation’s output for an experiment with two attracting bodies that collapses in finite time for two slightly different stratifications. Note the reasonable quantitative agreement. For the \(N\) body problem, use of the same modified Stokesian dynamics algorithm finds that self-assembly and particle aggregation strongly resemble the experimental observations in Figure 1 [7]. The low Péclet limit has mathematically exact flows and density perturbations for spheres with both insulating and absorbing boundary conditions. These solutions confirm that a complete flow reversal occurs for sufficiently diffusive spheres, revealing the possibility of much more complex self-assemblies.

Despite our discovery, we have barely scratched the surface of this rich class of problems. In the immediate future, we hope to assess the role of particle geometry, add non-insulating boundary conditions with different diffusivities between the solid and fluid, improve first-principles derivation of the effective force field, and study the details of the finite-time collapse that creates the self-assemblies. We are actively exploring the ramifications of many of these problems and anticipate that exciting new research in this field will keep us busy for many years to come.

**Acknowledgments:** We thank Daniel Harris, Robert Hunt, and Zeliha Kilic, our co-authors on the original work that appeared in* Nature Communications *[7]. Additionally, we received funding through the National Science Foundation (grant DMS-1910824) and the Office of Naval Research (grants ONR N00014-18-1-2490 and ONR N00014-23-1-2478).

**References**

[1] Abaid, N., Adalsteinsson, D., Agyapong, A., & McLaughlin, R.M. (2004). An internal splash: Levitation of falling spheres in stratified fluids. *Phys. Fluids*, *16*(5), 1567-1580.

[2] Allshouse, M.R., Barad, M.F., & Peacock, T. (2010). Propulsion generated by diffusion-driven flow. *Nat. Phys.*, *6*, 516-519.

[3] Brady, J.F., & Bossis, G. (1988). Stokesian dynamics. *Ann. Rev. Fluid. Mech.*, *20*(1), 111-157.

[4] Burd, A.B., & Jackson, G.A. (2009). Particle aggregation. *Ann. Rev. Mar. Sci.*, *1*(1), 65-90.

[5] Camassa, R., Falcon, C., Lin, J., McLaughlin, R.M., & Mykins, N. (2010). A first-principle predictive theory for a sphere falling through sharply stratified fluid at low Reynolds number. *J. Fluid Mech.*, *664*, 436-465.

[6] Camassa, R., Falcon, C., Lin, J., McLaughlin, R.M., & Parker, R. (2009). Prolonged residence times for particles settling through stratified miscible fluids in the Stokes regime.* Phys. Fluids*, *21*(3), 031702.

[7] Camassa, R., Harris, D.M., Hunt, R., Kilic, Z., & McLaughlin, R.M. (2019). A first-principle mechanism for particulate aggregation and self-assembly in stratified fluids. *Nat. Commun.*, *10*, 5804.

[8] Phillips, O.M. (1970). On flows induced by diffusion in a stably stratified fluid. *Deep Sea Res. Oceanogr. Abstr.*, *17*(3), 435-443.

[9] Wunsch, C. (1970). On oceanic boundary mixing. *Deep Sea Res. Oceanogr. Abstr.*, *17*(2), 293-301.

Roberto Camassa is the Kenan Distinguished Professor of Mathematics at the University of North Carolina (UNC) at Chapel Hill. He received the Kruskal Prize Lecture from the SIAM Activity Group on Nonlinear Waves and Coherent Structures in 2022. | |

Richard M. McLaughlin is a professor of mathematics at UNC Chapel Hill. He has chaired the Department of Mathematics for more than 10 years. Together, McLaughlin and Camassa run the Joint Fluids Laboratory at UNC, which hosts a 120-foot modular wave tank with state-of-the-art-scientific instrumentation and wave-making capabilities, as well as smaller devices that probe intriguing fluid phenomena. |