# Improved Prediction of Settling Behaviour of Solid Particles through Machine
Learning Analysis of Experimental Retention Time Data

Liron Simon Keren<EMAIL_ADDRESS>Turbulence Structure Laboratory, Tel
Aviv University, Israel  Teddy Lazebnik  Department of Cancer Biology, Cancer
Institute, University College London, UK  Alex Liberzon  Turbulence Structure
Laboratory, Tel Aviv University, Israel

###### Abstract

The motion of particles through density-stratified interfaces is a common
phenomenon in environmental and engineering applications. However, the
mechanics of particle-stratification interactions in various combinations of
particle and fluid properties are not well understood. This study presents a
novel machine-learning (ML) approach to experimental data of inertial
particles crossing a density-stratified interface. A simplified particle
settling experiment was conducted to obtain a large number of particles and
expand the parameter range, resulting in an unprecedented data set that has
been shared as open data. Using ML, the study explores new correlations that
collapse the data from this, and previous work Verso  _et al._ [1]. The “delay
time,” which is the time between the particle exiting the interfacial layer
and reaching a steady-state velocity, is found to strongly depend on six
dimensionless parameters formulated by ML feature selection. The data shows a
correlation between the Reynolds and Froude numbers within the range of the
experiments, and the best symbolic regression is based on the Froude number
only. This experiment provides valuable insights into the behavior of inertial
particles in stratified layers and highlights opportunities for future
improvement in predicting their motion.

Inertial particles, Lagrangian trajectories, Density interface, Stratification
force

††preprint: PRF/123-QED

##  I Introduction

Inertial spherical particles (radius $a$, density $\rho_{p}$) crossing
stratified density interfacial layers appear in various engineering and
environmental fluid mechanics problems [2, 3, 4]. Modeling accurately the time
it takes for a particle to move through a stratified layer would be
beneficial, for example, in problems of marine snow aggregations [5, 6],
airborne or waterborne pollutant dispersion [7, 8, 5], and oxygen levels
regulation in the ocean [9, 10].

In some parameter regimes of particle velocity, size, and density versus the
density and viscosity of the fluids, there is a significant amount of fluid
(sometimes called a caudal wake) that follows the particle into the
interfacial layer of different densities. The coupled dynamics of the particle
and of the caudal wake leads to additional resistance on the particle. We will
denote this additional resistance by $F_{s}$. There is an active discussion in
the community about the origin, nature, and magnitude of the additional
resistance force or forces [2, 1, 11, 12, 13, 14].

In this study, we do not focus on the force. Instead, we address this
phenomenon from a different perspective. Whatever the details of the
additional resistance are, if it is non-negligible, it delays particle
crossing. Let us consider the problem at hand schematically in Fig. 1. It
shows a density profile of the two homogeneous fluid layers at densities
$\rho_{1}$, $\rho_{2}$ with $\rho_{2}>\rho_{1}$, and a continuous smooth
transition between the two in the stratified interfacial layer of thickness
$h$ (green curve). We plot the typical particle velocity that is expected to
change from $V_{1}$ to $V_{2}$ as a smooth monotonic function (red curve) if
the additional resistance is negligible. According to the known equation of
motion, the particle moves according to the balance of gravity, buoyancy, and
viscous drag forces [2]. If the additional resistance is non-negligible, it
could lead to a non-monotonic change of velocity with a local minimum near the
edge of the interfacial layer [13, 1, 11, 14] (blue curve). In other words,
the observed retention time, defined as the time interval ($t_{V_{1}-V_{2}}$)
during which particle terminal velocity is lower than the steady-state values
($V_{1},V_{2}$), will be longer than the expected crossing time (i.e., without
the additional resistance, marked here as $\hat{t}_{V_{1}-V_{2}}$). In
addition, it is also possible that the particle will levitate, reversing its
motion [12, 1]. With or without levitation, there is a non-negligible delay
time, defined as $\tau$ in Fig. 1. In this study, we focus on predicting
$\tau$ based on particle and fluid properties.

Figure 1: Schematic description of the problem of a particle crossing a
stratified interfacial layer. The density of the fluid $\rho(z)$ changes
vertically throughout the medium (green), from the first layer $\rho_{1}$ to
the second layer $\rho_{2}$. The particle velocity $V(t)$ is illustrated for
the case of negligible (red) or significant (blue) stratification force. The
settling velocity in each layer is named $V_{1}$ and $V_{2}$, and the
interfacial layer thickness is $h$. We mark three relevant time intervals: the
expected crossing time $\hat{t}_{V_{1}-V_{2}}$, the retention time
$t_{V_{1}-V_{2}}$, and the delay time $\tau$.

There are only a few studies that address the problem of inertial particles of
diameter $a$ crossing a relatively sharp interfacial stratified layer
($h/a\sim O(10)$) between two miscible fluids of different densities,
$\rho_{1}$ and $\rho_{2}$ at so-called “entrance” Reynolds
$Re_{1}=V_{1}a/\nu>10$ and Froude $Fr_{1}=V_{1}/Na<100$ numbers, where $N$ is
calculated as Eq. (1). In Fig. 2, we summarize the existing map of the results
in terms of $Re_{1},Fr_{1}$.

Figure 2: Summary of the parameter range of the experimental results of Verso
_et al._ [1], Abaid  _et al._ [12], Srdić-Mitrović  _et al._ [13] and this
study, in terms of entrance Froude versus entrance Reynolds numbers. A-
present study, B - Verso  _et al._ [1], C - [12], and D - [13].

The first experimental study of inertial particles settling through an
interfacial layer of finite thickness is by Srdić-Mitrović  _et al._ [13]. The
authors used a water-alcohol-brine system and particles in the ranges
$1.5<Re_{1}<15$, $3<Fr_{1}<10$. The authors attributed the particle slowdown
to the additional drag force due to caudal fluid. In addition, they mentioned
the plausible contribution of internal waves or modification of flow structure
around the particle due to the density gradient [13]. However, the authors
studied only the time it took the spheres to reach minimum velocity and did
not investigate the phenomena in the denser bottom layer after the particles
crossed the interface.

Abaid  _et al._ [12] performed similar experiments with sharper interfaces in
the ranges of $1\leq h/a\leq 5$, $20<Re_{1}<450$, and $5<Fr_{1}<20$. Note that
the range does not mean that authors changed parameters systematically, it
only marks the minimal and maximal values. In Fig. 2, the pink rectangle marks
the parameter range of the reported measurements of [12]. Despite the similar
range of Reynolds numbers, the authors observed different caudal wake effects
and reported a “temporal levitation” phenomenon, where spheres momentarily
reverse their motion upwards, against gravity. Abaid  _et al._ [12] formulated
a closure model as an additional virtual mass with its own degree of freedom,
mimicking a caudal wake that moves upwards when the particle moves downwards.
The authors did report prolonged time periods until the sphere regained steady
state motion in the bottom layer, but did not study any time scales.

Verso  _et al._ [1] estimated the retention time based on the results of an
experiment similar to Srdić-Mitrović  _et al._ [13]. The authors used a water-
alcohol-brine system in ranges of $2<Re_{1}<14$, $0.6<Fr_{1}<4$, and $h/a\sim
10$. Verso  _et al._ [1] did not observe the levitation phenomenon in their
parameter range. However, they proposed a model for the additional resistance
due to stratification and caudal fluid that helped to estimate the parameter
range of $Re,Fr$, and $h/a$ in which the levitation is possible. In addition,
Verso  _et al._ [1] observed very prolonged $t_{V_{1}-V_{2}}$ timescales and
developed their parametric model for the stratification force $F_{s}$. The
authors also show that the crossing time $t_{V_{1}-V_{2}}$ is inversely
proportional to the particle Reynolds number in the bottom layer,
$Re_{2}=V_{2}a/\nu_{2}$. Their parametric model also predicted the data from
the literature [13].

Recently, Wang  _et al._ [15] reported in their preprint a detailed experiment
on the bouncing effect of particles, similar to the levitation reported by
Abaid  _et al._ [12]. The authors measured and numerically simulated velocity
fields and suggested a model of the resistance force. Unfortunately, their
data is not yet available to compare to our funding.

These studies are based on detailed measurements of a relatively low number of
particles, $5-50$, mainly due to the technical difficulties involved in such
complex experiments. The major challenge is controlling the thickness and
location of the interface without mixing the fluid layers, handling small
spheres ($a\sim 10\div 10^{3}\;\mu$m), tracking the particles with high
spatial resolution, and in some cases enforcing refractive index matching.
Although these experiments provide insight into the fluid mechanics and
dynamics of the particles inside the interfacial layers, they do not create a
statistically sufficient dataset to model the retention time, and limited in
terms of parameter regimes.

We also suspect that there are more useful forms of Reynolds and Froude
numbers using different combinations of velocity and length scales. It appears
that the present typical distinction of parameter regimes used in the
literature in terms of $Re_{1}$ and $Fr_{1}$ [2, 16, 17, 18, 19] is
incomplete. Examination of the map presented in Fig. 2 raises questions about
its specificity regarding the effect of stratification resistance. We could
expect a map on which it is more clear which particles experience different
physics (feeling a significant caudal wake resistance) and which do not (i.e.,
crossing with quasi-steady-state velocity values).

In this study, we propose another approach: we deliberately simplify the
experimental design, avoiding the difficulties associated with estimating the
force component, $F_{s}$, mainly the necessity for index refraction matching.
As a result, we can significantly extend the parameter space and include
previously unexplored regimes. In Fig. 3, we plot our data together with the
only data and existing correlation for the delay and retention times by Verso
_et al._ [1]. It is clear that we arrived at a different parameter regime for
which previous correlations do not match.

Figure 3: Normalized recovery time (a) and normalized delay time (b) versus
$Re_{2}$. The graphs include measurements collected during our experiments (A)
and those of particle type 1 from Verso  _et al._ [1] (B).

Furthermore, the simplified experimental design allowed us to obtain an
unprecedented number of particle trajectories. We use this sufficiently big
data set in the unexplored parameter regime with the custom-designed ML-based
symbolic regression tool, SciMED [20]. We developed this tool to find symbolic
regression correlations of $\tau$, using hidden non-dimensional parameters,
that might have the potential to better explain the underlying physical
mechanisms. This approach leads to an opportunity to find a new parametric
predictive model for $\tau$, using a “data-driven” methodology [20].

## Materials and Methods

In this section, we describe in depth the materials, equipment, measures, and
analysis techniques that were used, following the order of the scheme
presented in Fig.4. The unprecedented number of particle trajectories
collected enables the use of machine learning methods. For that reason, we
developed and applied a specific method that we abbreviated SciMED , described
in detail in the recent publication [20]. Our main interest at the end of the
process is to find a new correlation that fits the available data of $\tau$ as
a function of particle and fluid parameters.

Figure 4: Scheme of the methodology applied in this research.

### Experimental Setup

We conducted experiments in a glass tank with a cross-section of $200\times
200$ mm2 and a depth of $300$ mm, creating a top lighter layer of water
(marked as layer 1) above a bottom heavier layer of water and Epsom salt
(MgSO4) solution (marked as layer 2). We first fill the lighter fluid and pump
the heavy fluid from a valve at the bottom of the tank, resembling the method
used in Verso  _et al._ [21, 1]. This results in an interfacial stratified
layer between two fluid layers, created through molecular diffusion. The
interfacial layer is growing very slowly with time (less than 10 mm per day),
and in the present experiments, it is in the range $h=10a\div 100a$. We varied
the ratio of densities of the fluid layers $\rho_{2}/\rho_{1}$ in the range of
1.040 to 1.200. For these low concentrations of Epsom salts, it is customary
to ignore the small variation in fluid viscosity. Alike previous studies [13],
we assume equal kinematic viscosity of both layers,
$\nu_{1}=\nu_{2}=\nu=10^{-6}$ m2 s-1. The values of fluid density
($\rho_{1},\rho_{2}$) and Brunt-Vaisala frequency ($N$) are shown in Table 1,
where in this study, $N$ is calculated as follows:

$N=\sqrt{\frac{2g}{\rho_{1}+\rho_{2}}\frac{\rho_{2}-\rho_{1}}{h}}.$ (1)

Table 1: Ranges the properties of the fluid layers varied in, across all experiments. $\rho_{1}$ |  $\rho_{2}$ |  T |  $\hat{h}$ |  $N$ |  $\nu$
---|---|---|---|---|---
(g cm-3) |  (g cm-3) |  (${}^{\circ}C$) |  (cm) |  (s-1) |  (m2 s-1)
0.997 - 1.002 |  1.038 - 1.200 |  16 - 23.5 |  2.94 - 5.11 |  3.05 - 6.51 |  $1\times 10^{-6}$

Both the upper and lower fluid layers of the resulting medium are sufficiently
deep for particles to reach terminal velocity. We used one type of high-
quality commercially available spherical particles (Cospheric Inc., Santa
Barbara, CA). These spherical polystyrene particles have a density of 2.5 g
cm-3 and a diameter range of 1000–1180 $\mu$m. We verified particle parameters
using microscopy imaging and custom image processing code. We have found that
the particles are close to perfectly spherical, with a diameter that is
normally distributed within the prescribed range.

We released particles at the center of the tank, below the free surface, using
a syringe pump and a long glass pipette filled with water and particles. The
diameter of the pipette exit is slightly larger than the particle diameter,
ensuring that only single particles can exit. We controlled the rate of the
syringe pump and ensured that there was a sufficient time interval between the
releases. Thus the majority of the settling particles crossed an undisturbed
steady-state interface. Sometimes, however, we were unable to prevent
particles from concentrating at the exit of the pipette and exiting with
insufficient time intervals between them. In these cases, particles settle one
after another with vertical distances of about 10 - 100 diameters. The
insufficient time between two consecutive particles means that some particles
enter the interface at the same horizontal location, and the interface itself
is distorted by the previous particle.

We repeated the process 17 times, with each experiment taking between 3 and
5.5 hours (from the moment of filling the tank until the moment of the last
recorded measurement). We filmed the motion of the settling particles using a
digital high-speed camera (Photron Nova, 1024 $\times$ 1024 pixels) with a
frame rate of 250 frames per second.

### Trajectory tracking and data collection

We tracked the location of particles using the well-known particle tracking
algorithm of OpenPTV [22, 1]. As our experiments were conducted with a
relatively sharp interfacial layer, they have a significant drawback of
substantial variation of refractive index in the interfacial layer. As a
result, images of small particles from this region are significantly
distorted, as demonstrated in zone (b) in Fig. 5. Similar to other studies
[23], it is impossible to obtain particle center in this region with an
uncertainty smaller than a particle radius. Therefore, we could not measure
particle velocity within the interfacial layer and focused only on measuring
the information relevant to the delay time. Thus, we measured the time instant
at which each particle starts slowing down, slightly before it enters the
interfacial layer and the time instant at which it reaches the terminal
velocity at the bottom layer, far below the layer of the strong refractive
index gradient.

Figure 5: Digital snapshots of particle positions taken at different moments
throughout its motion, moving from left to right sequentially. Spherical
particle shapes in the left and right images indicate that the particle is
currently settling through a homogeneous fluid, while the smeared elliptical
shape in the middle images implies that the particle is in the virtual
interface layer $\hat{h}$.

In addition, we use the particles to estimate the position of the interface.
We mark the entrance and exit points for each trajectory based on particle
velocity (outside the region of strong refractive index gradients). The entry
point is marked at the location where $V(z)=0.99V_{1}$, and the exit point
where either the particle reaches its minimum velocity, or if it doesn’t
experience a minimum, the terminal velocity $V_{2}$, similar to Refs. [12, 13,
21, 1]. We average these position values for groups of particles that were
released in the same run and use this estimate as the approximate interface
thickness $\hat{h}$. Fig. 6 demonstrates the individual entrance and exit
points (dashed lines) of three archetypes of particle trajectories that we
measured, alongside the interface thickness $\hat{h}$ calculated on the basis
of them (blue rectangle). To the best of our knowledge, this is the first time
this method for estimating the position and thickness of the interface has
been used. Note that it is as accurate as other definitions of interface
thickness, such as sampling of fluids at fixed depths, intrusive measurements
with a conductance probe, or imaging methods using dye or Schlieren optical
methods [13, 15]. Our method avoids mixing and significantly prolongs the
experimental run time, measuring thousands of trajectories.

Figure 6: Three typical plots that can be found in the dataset, showing the
dimensional velocity of a particle versus the distance traveled from the top
of the viewing. These particular trajectories are taken from the same group of
particles that were released successively in the presence of a density ratio
of $\rho_{2}/\rho_{1}=1.040$. The dashed lines indicate the stratification
region as determined by each particle’s motion, while the colored area (zone
b) indicates the interfacial layer determined by the average of the located
entrance and exit points. Zones marked as a-c indicate the upper and lower
fluid layers.

As seen in Fig. 6, the three archetypes of particle trajectories are particles
with a clear local minimum velocity ($I$), particles without a minimum ($II$),
and particles that increase their velocity before entering the interfacial
layer ($III$). Due to our inability to define particle position in zone (b),
we “tracked” the particles in the sense that they are still linked to the same
object, but there is no useful quantitative information in this region.
Therefore, we masked the values tracked within the individually measured
interface and used trajectories only to obtain $V_{1}$, $V_{2}$, the position
of the minimum velocity, and the actual retention time $t_{V_{1}-V_{2}}$.
Specifically, to determine $V_{1},V_{2}$, we averaged the first and last
20-100 measurement points (at 250 fps), with the number of measurements to be
averaged depending on the position of the interfacial layer entrance in the
frame. For trajectories of archetype $III$, we adjusted the number of averaged
measurement points to avoid the time with increased velocities. For
trajectories of archetype $I$, where the trajectory fails to stabilize within
the limits of the frame, we determined $V_{2}$ as the average of the last five
measurement values of each trajectory. Moreover, the retention time is
determined by the time interval between the moment the particle returns to a
terminal velocity by reaching $0.99\leq V(t)/V_{2}\leq 1.01$ and the time
instant of departing from $V_{1}$ before entering the interfacial layer.

Next, similar to previous work [1], we estimate each particle diameter and
density using the measurements of its terminal velocities in the top and
bottom layers, and applying a force balance between the immersed weight,
$F_{WB}=(\rho_{p}-\rho_{f})\text{\makebox[0.0pt][l]{\ooalign{\hfil\hphantom{$V$}\hfil\cr\kern
0.80002pt--\hfil\cr}}}Vg$, and the drag $F_{D}=C_{D}\frac{1}{2}\rho_{f}|V|VA$.
Here $\rho_{f}$ is the density of the surrounding fluid (either $\rho_{1}$ or
$\rho_{2}$ in the homogeneous layers),
$\text{\makebox[0.0pt][l]{\ooalign{\hfil\hphantom{$V$}\hfil\cr\kern
0.80002pt--\hfil\cr}}}V=\pi a^{3}/6$ is the particle’s volume, $A=\pi a^{2}/4$
is the projected surface area, and $C_{D}$ is the drag coefficient. The drag
coefficient is estimated using the correlation [24, 1]:

$C_{D}=0.4+\frac{24}{Re}+\frac{6}{1+\sqrt{Re}}$ (2)

We estimate the diameter and particle density using constrained optimization
within the range provided by the manufacturer using the force balances in the
upper and lower layers simultaneously, similar to Ref. [1].

In analogy to previous studies looking for correlations [13, 12, 1, 25], we
also used two average properties: the average density of the fluid medium
$\overline{\rho}$, and the average of terminal velocities
$\overline{V}=0.5(V_{1}+V_{2})$. From these properties, we are able to
calculate the expected “crossing time”
$\hat{t}_{V_{1}-V_{2}}=\hat{h}/\overline{V}$ that represents the time it takes
the particle to cross $\hat{h}$ assuming that $F_{s}$ is negligible. Finally,
we calculated the delay times $\tau$ as the delta between the actual retention
time $t_{V_{1}-V_{2}}$ and the expected crossing time $\hat{t}_{V_{1}-V_{2}}$
($0\leq\tau<\infty$).

### Data cleaning and preparation of non-dimensional dataset

Footage of 3,264 settling particles is freely available as open access data,
divided into 16 experiment dates (one of the days had 2 experiments, thus 17
experiments in total). We processed and extracted the dimensional properties
and velocity trajectory of 2,039 particles from this database. The remaining
1,225 trajectories could not be tracked due to the technical limitations of
the image processing and tracking algorithms. For example, due to the sharp
refractive index gradient in zone (b) in Figs. 5-6, some trajectories were
misidentified as two separate trajectories; one starting from above the
interface and ending within it, and the other starting somewhere within the
interface and continuing until the end of the frame. We did not attempt to
match fractions of trajectories together and deleted them from the dataset.
Other examples of deleted trajectories are of particles that resulted in a
calculated expected time ($\hat{t}_{V_{1}-V2}$) that is greater than the
measured retention time ($t_{V_{1}-V2}$). By definition
$\hat{t}_{V_{1}-V2}\leq t_{V_{1}-V2}$ leading to $0\leq\tau$. However, since
we do not measure $\hat{t}_{V_{1}-V2}$ directly but estimate it using the
averaged interfacial layer thickness $\hat{h}$, some particles appear as they
have $\hat{t}_{V_{1}-V2}>t_{V_{1}-V2}$. In such cases, the measured retention
time is negative, and we discard these from the dataset used for the following
stage. In summary, we have a dataset of 2,039 trajectories for which
$0\leq\tau<\infty$.

The dimensional properties of all the particles in this dataset are presented
in Table. 2. From this dataset, 18% of the analyzed particles do not
experience a minimum velocity, and 50% exhibit a minimum that is significantly
lower than $V_{2}$ (more than a 5% decrease from the terminal velocity in the
bottom layer). This leaves 32% of particles with a minimum that is smaller
than $V_{2}$ by $0.01\%-4.5\%$. Due to the large amount of uncertainty
associated with the measurements, it is difficult to determine whether these
particles truly experience a minimum or not.

Table 2: Ranges the properties of the particles and of the dimensionless numbers of the top and bottom layers varied in, across all experiments. $a$ |  $\rho_{p}$ |  $V_{1}$ |  $V_{2}$ |  $\tau$
---|---|---|---|---
(mm) |  (g cm-3) |  (cm s-1) |  (cm s-1) |  (s)
1 - 1.18 |  2,425 - 2,575 |  11.06 - 21.48 |  9 - 17.43 |  0 - 0.49

The dataset comprises of 9 dimensional properties of particles and fluids.
Prior to symbolic regression of a correlation for $\tau$, we initially create
a dataset of non-dimensional parameters.

According to the Buckingham theorem of dimensional analysis, the number of
non-dimensional parameters $\pi_{i}$ equals the number of relevant dimensional
variables minus the number of dimensional units. In our case, the number of
dimensional units is 3: either mass, length, and time, or force, length, and
time. We estimated that we need six non-dimensional parameters to formulate a
correlation for the delay time, listed in Table. 3.

Second, we had to select the form of every dimensionless parameter. This is
because every dimensionless parameter $\pi_{i}$ can be calculated in various
options based on variables with the same dimensions (between 1 and 384
different options depending on the complexity of the dimensionless parameter).
For instance, the Reynolds-like parameters can be defined using one of the two
characteristic length scales ($a,\hat{h}$) and three different velocities
($V_{1},V_{2},\overline{V}=0.5(V_{1}+V_{2})$), resulting in six plausible
options. Each of these options will lead to a different meaning of the
parameter and different relative contributions of the physical properties in
that parameter. Table. 3 includes the number of options for each parameter.

Third, we needed a non-dimensional time scale based on the delay time $\tau$
defined above. For that, we used the expected crossing time
$\hat{t}_{V_{1}-V_{2}}$ for the normalization. This value can be obtained for
the particle of a known size and density for a given interfacial layer.
Therefore, it is useful for future prediction applications. The normalized
delay time scale is then $\hat{\tau}=\tau/\hat{t}_{V_{1}-V_{2}}$.
Theoretically, it can be any value in the range $0\leq\hat{\tau}<\infty$, with
infinite values corresponding to particles entrapped in the interfacial layer.

### Selection of six non-dimensional parameters through a feature selection
process

The first stage in symbolic regression (SR) was to choose the “best” option
for each of the six $\pi_{i}$ parameters (e.g., Reynolds-like, Froude-like,
etc.) out of a total of 412 different options presented (see No. of options
column in Table 3). This is a major component of the proposed methodology,
called feature selection, and it is only available due to the large dataset of
measurements. We implemented the process in SciMED , which selects six
features (each feature being one of the options for each $\pi_{i}$) that are
“most informative” to the dimensionless delay time. “Most informative” means
that an ML model for “black box” prediction, generated the most accurate
predictions of $\hat{\tau}$ using the data of the selected six features,
compared to all other feature combinations.

During the feature selection process, SciMED incorporates the information
about the category of each feature, meaning to which $\pi_{i}$ parameter it
can be assigned. Then, it uses a genetic algorithm (GA) to generate tens of
thousands of subsets of the original dataset, where each subset includes only
six features, one from each $\pi_{i}$. The subsets are evaluated and compared
by training an AutoML model [20] that predicts $\hat{\tau}$ based on the
selected set of features. The score of each subset is determined by the
accuracy of the prediction that was achieved. Finally, the subset leading to
the most accurate result is passed as the dataset for the SR process.

### Symbolic regression for the delay time

In this step, SciMED replaces the scientist, and instead of manually searching
for the best fit, it suggests the optimal equation structure and parameters
that best predict $\hat{\tau}$ of all particles in the dataset. For that
purpose, the Las Vegas \- based SR component [20] searches for the optimal
correlation. This step results in a list of three possible equations, ordered
by their complexity (in terms of the number of operators and parameters in
it).

Each equation consists of a numerical prefactor $\alpha$ and a constant term
$\beta$ that is optimized for our dataset of samples (i.e., $\hat{\tau}=\alpha
f(\pi_{1},...,\pi_{6})+\beta$). For cases of data collected from particles or
fluid medium that differ from this experiment, such as the data by Verso  _et
al._ [1], we run another optimizer to obtain the $\alpha$ and $\beta$
numerical terms. Eventually, for the sake of generalization, we suggest one
ML-selected correlation (i.e., equation resulting from the SR) that best fits
both our data and the data from Verso  _et al._ [1].

## Results

### Feature selection

We report here the set of six $\pi_{i}$ that result from a feature selection
process, comparing 43,000 variations. This set supposedly comprises the
dimensionless numbers representing the dominant mechanisms determining the
normalized delay time $\hat{\tau}$.

In Table 3, we present a comparison between the most typical choice in the
literature and the ML-selected form of each $\pi_{i}$ parameter. In the right
column, we present the value range for each parameter according to the
selected form and our measurements.

Table 3: Comparison of the six non-dimensional parameters in the conventional form and the form selected by SciMED . For the length and velocity scales, there is only one reasonable selection. The ranges are for the selected form. $\pi_{i}$ |  No. of options |  Conventional form |  Selected form | Range
---|---|---|---|---
Length scales |  1 |  $h/a$ |  $h/a$ | 25 - 50
Characteristic velocity |  1 |  $V_{1}/V_{2}$ |  $V_{1}/V_{2}$ | 1 - 1.8
$Re$-like |  6 |  $V_{1}a/\nu$ |  $\overline{V}a/\nu$ | 103 - 213
Characteristic density |  4 |  $\rho_{2}/\rho_{1}$ |  $\rho_{p}/\overline{\rho}$ | 2.20 - 2.53
Density jump |  16 |  $\rho_{2}/(\rho_{2}-\rho_{1})$ |  $\rho_{2}/(\rho_{p}-\overline{\rho})$ | 0.67 - 0.91
$Fr$-like |  384 |  $V_{1}/a\sqrt{\frac{g}{h}\frac{\rho_{2}-\rho_{1}}{\overline{\rho}}}$ |  $V_{1}/h\sqrt{\frac{g}{a}\frac{\rho_{p}-\rho_{2}}{\overline{\rho}}}$ | 0.02 - 0.06

In Fig. 7, we compare $\hat{\tau}=f(\pi_{i})$ in the conventional (top row)
and the selected (bottom row) form of four different $\pi_{i}$ parameters. The
dimensionless parameters of length and velocity are not presented in this
comparison as, in this case, there is no other plausible form to calculate
them. As seen in the figure, the collapse of data of the $\pi_{i}$ parameters
versus the normalized delay time leads to improved correlation, although weak
for each dimensionless parameter separately. Furthermore, the parameters
leading to the most explicit trends are the selected $Re$-like and $Fr$-like
ratios, which from here on now, will be named $\hat{Re}$ and $\hat{Fr}$.

Figure 7: Normalized delay time $\hat{\tau}$ versus four different $\pi_{i}$
parameters in their conventional form (top row) and the from selected by
SciMED (bottom row).

In Fig. 8, we present our measurements (marked as A, for either particle that
experiences a significant velocity minimum or those that do not) and that of
Verso  _et al._ [1] (marked as B, varying from type P1 to P4), over the
$\hat{Fr}=f(\hat{Re})$ map, similarly to Fig. 2. We did not plot the
measurements of Abaid  _et al._ [12] and Srdić-Mitrović  _et al._ [13] that
were presented in Fig. 2, as we do not have an access to the data needed to
calculate $\hat{Fr}$. Note that particle type P3 corresponds to the case that
did not show minima (see Verso  _et al._ [1]). This new map emphasizes for the
first time that in terms of the new dimensionless parameters $\hat{Fr}$ and
$\hat{Re}$, all the particle types are different, a fact that was not
observable in the conventional $Fr_{1}=f(Re_{1})$ map [1]. Note that for our
measurements, there is also a notable though the small difference between the
particles without and with significant minima.

Figure 8: Parameter range of the experimental results of this study (A) and of
four types of particles presented in Verso  _et al._ [1] (B P1-P4), in terms
of $\hat{Fr}$ versus $\hat{Re}$ (i.e., the $Re$ and $Fr$-like parameters
selected by SciMED ).

### Symbolic regression and new correlations

Table. 4 presents the three types of equations suggested by the SR component
of SciMED together with the numerical constants optimized to our dataset. The
three typically used success metrics of $R^{2}$, mean absolute error (MAE) and
mean squared error (MSE) presented were calculated from all samples of our
study together with those of Verso  _et al._ [1]. As seen in the table, all
the suggested correlations resulted in similar scores. It is difficult to
decide which correlation is “better,” but this is also not the purpose of this
study. In order to obtain a much better correlation, we need to obtain more
measurements of very different types of particles, with different fluids,
interfacial layer thickness, and density jumps. To this end, we obtained
symbolic expressions representing the result better than the previous studies.
Furthermore, we focus on the question of why these correlations are better
than the previous ones and what underlying physical mechanisms these
selections could highlight.

Table 4: Three plausible correlations for the normalized delay time, as suggested by SciMED . The $R^{2}$, mean absolute error (MAE), and mean squared error (MSE) metric scores are reported based on measurements of this study and that of Verso  _et al._ [1]. |  Equation |  Constants |  $R^{2}$ |  MAE |  MSE
---|---|---|---|---|---
$f_{1}$ |  $\alpha_{1}\big{(}\hat{Fr}\frac{\rho_{p}}{\overline{\rho}}\big{)}+\beta_{1}$ |  $\alpha_{1}=16.19,\;\beta_{1}=-0.5$ |  $0.75$ |  0.24 |  0.23
$f_{2}$ |  $\alpha_{2}\big{(}\hat{Fr}\hat{Re}\frac{h}{a}\big{)}+\beta_{2}$ |  $\alpha_{2}=2.76,\;\beta_{2}=0.47$ |  $0.71$ |  0.23 |  0.20
$f_{3}$ |  $\alpha_{3}\big{(}\hat{Fr}^{2}\frac{\rho_{p}}{\overline{\rho}}(1+\frac{\rho_{p}}{\overline{\rho}}\frac{V_{2}}{V_{1}})\big{)}+\beta_{3}$ |  $\alpha_{3}=89.91,\;\beta_{3}=0.44$ |  $0.72$ |  0.24 |  0.23
Figure 9: Predicted normalized delay time as calculated by $f_{1}-f_{3}$
versus the actual measurements of $\hat{\tau}$. The scatter is the
experimental results of this study (A) and of four types of particles
presented in Verso  _et al._ [1] (B P1-P4).

In Fig. 9, we show the scatter of the predicted normalized delay time (using
$f_{1}-f_{3}$) versus the actual measurements of $\hat{\tau}$. In each plot, a
linear trend line (dashed line) is fitted to all predicted samples. It is
clear that the trend line fitted to the predictions of $f_{1}$ (left) is
closest to the expected $y=x$ trend (solid line). It is also noteworthy that
the results were obtained with particles measured in two separate sets of
experiments (ours and that of [1]), using two different fluid layer
combinations, five particle types, and 10 density jumps. Based on them we
suggest to use $f_{1}$ as the correlation for $\hat{\tau}$, which can be
written as:

$\displaystyle\hat{N}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{g}{a}\frac{\rho_{p}-\rho_{2}}{\overline{\rho}}}$ (3)
$\displaystyle f_{1}$ $\displaystyle=$
$\displaystyle\alpha\frac{\rho_{p}}{\overline{\rho}}\,\frac{V_{1}}{\hat{h}\hat{N}}+\beta$
(4)

and if we substitute $\hat{N}$ it reads:

$f_{1}=\alpha\frac{V_{1}}{h}\sqrt{\frac{a}{g}\frac{\rho_{p}}{\rho_{p}-\rho_{2}}\frac{\rho_{p}}{\overline{\rho}}}+\beta$
(5)

## Discussion

In this study, we took a new approach to reveal the structure of empirical
results. We utilized machine learning (ML) techniques to study the complex
behavior of inertial particles crossing an interfacial stratified layer. We
simplified the experimental setup to gain a substantial expansion of the
previously unexplored parameter range and an unprecedented number of particle
trajectories dataset. With this data, we aimed to find a correlation that
covers our and previous results. All the experimental results are open to the
community, and we hope that additional ideas and new correlations can be found
through new approaches.

We also find, consistently with previous results [13, 12, 14, 1], that not all
particles (in our case, 18% of all the trajectories) in this parameter range
experience slowdown and minimal velocity. Unfortunately, we do not have a good
explanation for this effect. We suspect that particle rotation, an
insufficient time between particle releases, non-spherical shapes, or maybe an
oscillating interface, could cause the discrepancy. New experiments would
probably be necessary to understand this problem better.

We demonstrate that the normalized delay time ($\hat{\tau}$) has a stronger
dependence on the newly selected forms of dimensionless parameters, shown in
Fig. 9. We also inspect the selected forms of the dimensionless parameters and
infer the possible meaning of the selection. Thus, $\hat{Re}$ is based on
average velocity and not on the entrance velocity as it was suggested by
previous works [13, 1]. The Reynolds number does depend on the particle size,
and it represents the so-called particle Reynolds number.

The particle density is normalized to the averaged fluid density (which is the
approximate density of the interface, assuming linear density gradient). It
seems to be more reasonable to incorporate both densities instead of the
conventional $\rho_{2}/\rho_{1}$.

We also learn that because there is a strong correlation between the Reynolds
and Froude numbers (see the new map in Fig. 8), the final correlation does not
include the Reynolds number explicitly. The most dominant effect in this
problem is stratification, and it determines the retention and the delay
times. Therefore, the most dominant dimensionless parameter here is the
Froude-like parameter. Furthermore, the newly selected Froude number differs
from the one typically suggested in the literature. It relates linearly to the
interfacial layer thickness, $h$, and only as a square root of the particle
diameter $a^{-0.5}$. In other words, we can infer that the dimensionless
frequency, $\hat{N}$ should be defined with the particle diameter, while
Froude number as $\hat{Fr}=V/\hat{h}\,\hat{N}$. The better choice of this form
for the Froude number is also supported by the better collapse of data than
the rest of the $\pi_{i}$ parameters. We also compared it to two Froude number
forms and see that it is significantly stronger related to the delay time,
i.e. ($\hat{\tau}=f(\hat{Fr})$ leads to $R^{2}=0.19$, while both $Fr_{1}$ and
$Fr_{2}$ lead to an order of magnitude weaker correlation with $R^{2}=0.01$.

Undoubtedly, our experiment has limitations. Although it enabled the sampling
of thousands of trajectories, it also led to increased uncertainty and lesser
control over particle diameter, density, and less precise release timing. The
main disadvantage is the lack of index refraction matching and the
corresponding lack of detailed trajectories inside the interfacial layer. If
such an experiment can be performed, it would be possible to implement SciMED
with respect to particle trajectories and compare the point-wise position and
velocity of each particle versus the predicted equation of motion. We expect
this approach to significantly improve our ability to predict the correct form
of the stratification force and formulate a model in a more general form of
detailed physical mechanisms.

## Acknowledgments

The authors thank Daniel Jalontzki for collecting the settling particle data,
Aviv Littman for measuring density profiles, and Dmitry Aronovsky for the
experimental setup and assistance with experiments.

## Declarations

### Funding

This research was partially supported by Israel Science Foundation (grant
number 441/2) and Gordon center for energy research at Tel Aviv University.

### Data availability

The data in this study is available in the manuscript with the relevant
sources.

### Conflicts of interest

The authors have no financial or proprietary interests in any material
discussed in this article.

## References

  * Verso  _et al._ [2019] L. Verso, M. van Reeuwijk, and A. Liberzon, Transient stratification force on particles crossing a density interface, International Journal of Multiphase Flow 121, 103109 (2019).
  * Magnaudet and Mercier [2020] J. Magnaudet and M. Mercier, Particles, drops, and bubbles moving across sharp interfaces and stratified layers, Annual Review of Fluid Mechanics 52, 61 (2020).
  * Mrokowska [2018] M. M. Mrokowska, Stratification-induced reorientation of disk settling through ambient density transition, Scientific Reports 8, 1 (2018).
  * Renggli  _et al._ [2016] C. J. Renggli, S. Wiesmaier, C. P. De Campos, K. U. Hess, and D. B. Dingwell, Magma mixing induced by particle settling, Contributions to Mineralogy and Petrology 171, 96 (2016).
  * MacIntyre  _et al._ [1995] S. MacIntyre, A. Alldredge, and C. Gotschalk, Accumulation of marines now at density discontinuities in the water column, Limnology and Oceanography 40, 449 (1995).
  * Prairie  _et al._ [2013] J. C. Prairie, K. Ziervogel, C. Arnosti, R. Camassa, C. Falcon, S. Khatri, R. M. McLaughlin, B. L. White, and S. Yu, Delayed settling of marine snow at sharp density transitions driven by fluid entrainment and diffusion-limited retention, Marine Ecology Progress Series 487, 185 (2013).
  * Turco  _et al._ [1983] R. Turco, O. Toon, T. Ackerman, J. Pollack, and C. Sagan, Nuclear winter: Global consequences of multiple nuclear explosions, Science 222, 1283 (1983).
  * Kok [2011] J. Kok, A scaling theory for the size distribution of emitted dust aerosols suggests climate models underestimate the size of the global dust cycle, Proceedings of the National Academy of Sciences 108, 1016 (2011).
  * Smith  _et al._ [1992] D. Smith, M. Simon, A. Alldredge, and F. Azam, Intense hydrolytic enzyme activity on marine aggregates and implications for rapid particle dissolution, Nature 359, 139 (1992).
  * Burd and Jackson [2009] A. Burd and G. Jackson, Particle aggregation, Annual review of marine science 1, 65 (2009).
  * Verso [2020] L. Verso,  _Particles crossing density interfaces_ , Ph.D. thesis, Tel Aviv University (2020).
  * Abaid  _et al._ [2004] N. Abaid, D. Adalsteinsson, A. Agyapong, and R. McLaughlin, An internal splash: Levitation of falling spheres in stratified fluids, Physics of Fluids 16, 1567 (2004).
  * Srdić-Mitrović  _et al._ [1999] A. Srdić-Mitrović, N. Mohamed, and H. Fernando, Gravitational settling of particles through density interfaces, Journal of Fluid Mechanics 381, 175–198 (1999).
  * Mandel  _et al._ [2020] T. Mandel, L. Waldrop, M. Theillard, D. Kleckner, S. Khatri, _et al._ , Retention of rising droplets in density stratification, Physical Review Fluids 5, 124803 (2020).
  * Wang  _et al._ [2023] S. Wang, P. Kandel, J. Deng, C. Caulfield, and S. Dalziel,  (2023).
  * Blanchette and Shapiro [2012] F. Blanchette and A. Shapiro, Drops settling in sharp stratification with and without marangoni effects, Physics of Fluids 24, 042104 (2012).
  * Doostmohammadi  _et al._ [2014] A. Doostmohammadi, S. Dabiri, and A. M. Ardekani, A numerical study of the dynamics of a particle settling at moderate reynolds numbers in a linearly stratified fluid, Journal of Fluid Mechanics 750, 5 (2014).
  * Yick  _et al._ [2009] K. Yick, C. Torres, T. Peacock, and R. Stocker, Enhanced drag of a sphere settling in a stratified fluid at small reynolds numbers, Journal of Fluid Mechanics 632, 49 (2009).
  * Camassa  _et al._ [2010] R. Camassa, C. Falcon, J. Lin, R. M. McLaughlin, and N. Mykins, A first-principle predictive theory for a sphere falling through sharply stratified fluid at low Reynolds number, Journal of Fluid Mechanics 664, 436 (2010).
  * Simon Keren  _et al._ [2023] L. Simon Keren, A. Liberzon, and T. Lazebnik, A computational framework for physics-informed symbolic regression with straightforward integration of domain knowledge, Scientific Reports 13, 1249 (2023).
  * Verso  _et al._ [2017] L. Verso, M. van Reeuwijk, and A. Liberzon, Steady state model and experiment for an oscillating grid turbulent two-layer stratified flow, Physical Review Fluids 2, 104605 (2017).
  * OpenPTV consortium [2014] OpenPTV consortium, Open source particle tracking velocimetry (2014).
  * Alahyari and Longmire [1994] A. Alahyari and E. Longmire, Particle image velocimetry in a variable density flow: application to a dynamically evolving microburst, Experiments in Fluids 17, 434 (1994).
  * White [1974] F. White, Viscous fluid flow mcgraw-hill inc, New York, New York (1974).
  * Boetti and Verso [2022] M. Boetti and L. Verso, Force on inertial particles crossing a two layer stratified turbulent/non-turbulent interface, International Journal of Multiphase Flow , 104153 (2022).