Improved prediction of settling behavior of solid particles through machine learning analysis of experimental retention time data

1. Introduction

Settling of inertial particles across layers of fluids of different densities appears in various engineering and environmental fluid mechanics problems (Magnaudet and Mercier, 2020; Mrokowska, 2018; Renggli et al., 2016). The settling velocity of particles can be estimated from the first principles only for cases of particles settling with low Reynolds (𝑅𝑒_𝑝 = 𝑎𝑉_𝑝∕𝜈, 𝑎 is particle diameter, 𝜈 fluid kinematic velocity) and low Stokes number (𝑆𝑡 = 𝑇_𝑝∕𝑇_𝑓, 𝑇_𝑝 particle response time scale, 𝑇_𝑓 flow response time scale). Furthermore, to be able to estimate the settling velocity, the particle must also settle through a fluid of homogeneous density (𝜌 = const) or weakly linearly stratified fluids (𝜕𝜌(𝑧)∕𝜕𝑧 = const, where 𝑧 is fluid depth). When inertial particles cross interfacial layers, which are fluid layers with sharp density changes, there may be a noticeable amount of lighter fluid that follows the particle into interfacial layer with different densities. This lighter fluid is sometimes referred to as a ‘‘caudal wake’’. The coupled dynamics of the particle motion and of its caudal wake with respect to the surrounding fluid, together with the flow due to particle motion, lead to additional resistance on the particle. We will denote this additional resistance force by 𝐹_𝑠 as a single quantity, although there is an active discussion in the community about the origin, nature, and magnitude of various sources of the resistance force (Magnaudet and Mercier, 2020; Verso et al., 2019; Verso, 2020; Abaid et al., 2004; Srdić-Mitrović et al., 1999; Mandel et al., 2020). Modeling accurately the time it takes for an inertial particle at a higher Reynolds and Stokes number to move through a stratified fluid with sharp density changes (hereinafter called interfacial layers) would be beneficial, for example, in problems of marine snow aggregations (MacIntyre et al., 1995; Prairie et al., 2013), airborne or waterborne pollutant dispersion (Turco et al., 1983; Kok, 2011; MacIntyre et al., 1995), and oxygen levels regulation in the ocean (Smith et al., 1992; Burd and Jackson, 2009).

Estimating the force components for different parameter regimes is extremely challenging due to the complexity of the dynamics of all the components, including the particle, the fluid layers, and the wake. Instead of integrating the force model in time, we suggest modeling the settling time directly, combining the experimental and machine learning methods.

1.1. Problem definition

Let us consider the problem at hand schematically in Fig. 1 . In panel (a), we show the scheme of the physical process: a particle heavier than 0301-9322/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

the fluid is settling from the top, lighter fluid layer, through the interfacial layer of thickness ℎ, and into the heavier fluid layer at the bottom (i.e., 𝜌₁ < 𝜌₂ < 𝜌_𝑝). In panel (b), we plot the fluid density profile of the two homogeneous fluid layers at densities 𝜌₁, 𝜌₂, and a continuous smooth transition 𝜌(𝑧) between the two homogeneous layers, called the stratified interfacial layer of thickness ℎ (green curve). We also plot the curve corresponding to the typical particle velocity expected to vary from 𝑉₁ to 𝑉₂ as a smooth monotonic function (cyan curve). This curve describes the case of a sphere settling without any additional resistance force stemming from the stratification, 𝐹_𝑠 = 0 (Magnaudet and Mercier, 2020). In some parameter regimes, the additional resistance is non-negligible, and there is possibly a non-monotonic change of velocity with a local minimum near the edge of the interfacial layer (Srdić- Mitrović et al., 1999; Verso et al., 2019; Verso, 2020; Mandel et al., 2020) (blue curve).

In this study, we focus on the settling time estimate. Suppose the additional resistance force in the stratified interfacial layer is negligible. In that case, we will measure the theoretically predictable settling time of the particle,̂ 𝑡_{𝑉1−𝑉2}, defined as the time the particle moves from one homogeneous-density layer to another. Instead of settling time, we can determine the so-called retention time; the interval during which particle velocity changes from one terminal velocity value, 𝑉₁ to another terminal velocity value, 𝑉₂, marked as 𝑡_{𝑉1−𝑉2}. This definition is more beneficial for several cases, such as cases when the particle is heavier than both fluids 𝜌_𝑝 > 𝜌₂, when the particle is lighter 𝜌_𝑝 < 𝜌₁ and rises through the interfacial layer, as well as for the cases when particle temporarily changes its direction of motion and levitates (Abaid et al., 2004; Verso et al., 2019).

For the monotonic, theoretically predictable case, the theoretically predictable settling time through the interfacial layer is equivalent to the retention time i.e.,̂ 𝑡_{𝑉1−𝑉2} = 𝑡_{𝑉1−𝑉2}. Suppose there is an unknown additional resistance force. In that case, the retention time is longer because it also contains the interval during which particle settling velocity is lower than both the steady-state values or it levitates and changes the velocity sign. We define the difference between the expected settling time and observed retention time as the ‘‘delay time’’, marked as 𝜏 (see Fig. 1). In this study, we develop the method to predict 𝜏 based on particle and fluid properties, using a machine learning model trained on experimental data.

1.2. Existing data

Only a few studies address the problem of inertial particles crossing sharp interfacial stratified layers (ℎ∕𝑎 ∼ 𝑂(10), ℎ is the interfacial layer thickness, and 𝑎 is the particle diameter) between two miscible fluids of densities 𝜌₁ and 𝜌₂. The key parameters are the ‘‘entrance’’ Reynolds number, defined with the particle size and the velocity and viscosity of the layer from which the particle enters the interfacial layer: 𝑅𝑒₁ = 𝑉₁𝑎∕𝜈> 10, and the corresponding entrance Froude number, 𝐹 𝑟₁ = 𝑉₁∕𝑁𝑎< 100, where 𝑁 is calculated as defined below in Eq. (1). In Fig. 2, we summarize the existing results on the map of 𝑅𝑒₁, 𝐹 𝑟₁. The figure also demonstrates the limited number of studies in the literature that referred to this parameter range, in which 𝑅𝑒₁ > 10, i.e., the particles have significant inertia, and at the same time, the stratification is relatively strong, i.e., 𝐹 𝑟₁ < 100. Outside of this parameter range, there are many more studies and comprehensive results. Most of these results are for cases of linear stratification and small particles at creeping flow regime, where 𝑅𝑒₁ < 1

In the suggested parameter regime, the first experimental study of inertial particles settling through an interfacial layer of finite thickness is by Srdić-Mitrović et al. (1999). The authors used a water–alcohol– brine system and particles in the ranges 1.5 < 𝑅𝑒₁ < 15, 3 < 𝐹 𝑟₁ < 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 (Srdić-Mitrović et al., 1999). 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. Therefore, we do not have their estimate of the retention time (as it requires the tracking of the particle within the bottom layer), and thus cannot infer the particles’ delay time 𝜏. This prevents us from using their measurements to train or validate our model.

Abaid et al. (2004) performed similar experiments with sharper interfaces in the ranges of 1 ≤ ℎ∕𝑎 ≤ 5, 20 < 𝑅𝑒₁ < 450, and 5 < 𝐹 𝑟₁ < 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 Abaid et al. (2004). Despite the similar range of Reynolds numbers to those of Srdić-Mitrović et al. (1999), 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. (2004) 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. (2019) estimated the retention time based on the results of an experiment similar to Srdić-Mitrović et al. (1999). The authors used a water–alcohol–brine system in ranges of 2 < 𝑅𝑒₁ < 14, 0.6 < 𝐹 𝑟₁ < 4, and ℎ∕𝑎 ∼10. Verso et al. (2019) 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 𝑅𝑒, 𝐹 𝑟, and ℎ∕𝑎 in which the levitation is possible. In addition, Verso et al.

(2019) observed very prolonged 𝑡_{𝑉1−𝑉2} timescales and developed their parametric model for the stratification force 𝐹_𝑠. The authors also show that the crossing time 𝑡_{𝑉1−𝑉2} is inversely proportional to the particle Reynolds number in the bottom layer, 𝑅𝑒₂ = 𝑉₂𝑎∕𝜈₂. Their parametric model also predicted the data from the literature (Srdić-Mitrović et al., 1999). This is the only data that we could use for training, along with the data measured in our experiments.

Recently, Wang et al. (2023) reported in their preprint a detailed experiment on the bouncing effect of particles, similar to the levitation reported by Abaid et al. (2004). The authors measured and numerically simulated velocity fields and suggested a model of the resistance force. Unfortunately, their data is not yet available for comparison.

The studies mentioned above are based on detailed measurements of a relatively low number of particles, 5−50, mainly due to the technically challenging experiments. The major challenges relate to control of the thickness and location of the interface without mixing the fluid layers, handling small spheres (𝑎 ∼10 ÷ 10³ μ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 delay times.

We also suspect that there are more useful forms of Reynolds and Froude numbers using different velocity and length scale combinations. It appears that the present typical distinction of parameter regimes used in the literature in terms of 𝑅𝑒₁ and 𝐹 𝑟₁ (Magnaudet and Mercier, 2020; Blanchette and Shapiro, 2012; Doostmohammadi et al., 2014; Yick et al., 2009; Camassa et al., 2010) 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, 𝐹_𝑠, particularly the need for index refraction matching. Note that the correlation between Fs and index refraction stems from the need to precisely track the particles within the interfacial layer to measure Fs. This entails meticulously matching the refractive indexes of the top and bottom layers, ensuring minimal refractive differences throughout the particle’s trajectory. By avoiding the need to match the refraction index of the top and bottom fluid layers, we can simplify the experimental design and significantly extend the parameter space to include previously unexplored regimes.

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 (Simon Keren et al., 2023). We developed this tool to find symbolic regression correlations of 𝜏, 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 𝜏, using a ‘‘data-driven’’ methodology (Simon Keren et al., 2023).

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. (2019). It is clear that we arrived at a different parameter regime for which previous correlations do not match.

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 (Simon Keren et al., 2023). Our main interest at the end of the process is to find a new correlation that fits the available data of 𝜏 as a function of particle and fluid parameters.

Experimental setup

We conducted experiments in a glass tank with a cross-section of 200 × 200 mm² 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. (2017, 2019). 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 ℎ = 10𝑎 ÷ 100𝑎. We varied the ratio of densities of the fluid layers 𝜌₂∕𝜌₁ 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 (Srdić- Mitrović et al., 1999), we assume equal kinematic viscosity of both layers, 𝜈₁ = 𝜈₂ = 𝜈 = 10⁻⁶ m² s⁻¹. The values of fluid density (𝜌₁, 𝜌₂) and Brunt–Vaisala frequency (𝑁) are shown in Table 1, where in this study, 𝑁 is calculated as follows:

𝜌1	𝜌2	T̂		ℎ		𝑁	𝜈
(g cm−3)	(g cm−3)	(◦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 × 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 (Co-spheric Inc., Santa Barbara, CA). These smooth spherical polystyrene particles have a density of 2.5 g cm⁻³ and a diameter range of 1000– 1180 μ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 h (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 × 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 (OpenPTV consortium, 2014; Verso et al., 2019). 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 (Alahyari and Longmire, 1994), 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.

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 𝑉 (𝑧) = 0.99𝑉₁, and the exit point where either the particle reaches its minimum velocity, or if it does not experience a minimum, the terminal velocity 𝑉₂, similar to Refs. Abaid et al. (2004), Srdić-Mitrović et al. (1999), Verso et al. (2017, 2019). We average these position values for groups of particles that were released in the same run and use this estimate as the approximate interface thicknesŝ ℎ. In other words,̂ ℎ is an average of individual interface thickness measurements, collected from trajectories of particles that were released in proximity to one another. Fig. 6 demonstrates the individual entrance and exit points (dashed lines) of three archetypes of particle trajectories that we measured, alongside the interface thicknesŝ ℎ calculated based on 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 (Srdić-Mitrović et al., 1999; Wang et al., 2023). Our method avoids mixing and significantly prolongs the experimental run time, measuring thousands of trajectories.

As seen in Fig. 6, the three archetypes of particle trajectories are particles with a clear local minimum velocity (𝐼), particles without a minimum (𝐼𝐼), and particles that increase their velocity before entering the interfacial layer (𝐼𝐼𝐼). 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 𝑉₁, 𝑉₂, the position of the minimum velocity, and the actual retention time 𝑡_{𝑉1−𝑉2}. Specifically, to determine 𝑉₁, 𝑉₂, 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 𝐼𝐼𝐼, we adjusted the number of averaged measurement points to avoid the time with increased velocities. For trajectories of archetype 𝐼, where the trajectory fails to stabilize within the limits of the frame, we determined 𝑉₂ 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 ≤ 𝑉 (𝑡)∕𝑉₂ ≤ 1.01 and the time instant of departing from 𝑉₁ before entering the interfacial layer.

Next, similar to previous work (Verso et al., 2019), 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, 𝐹_{𝑊 𝐵} = (𝜌_𝑝−𝜌_𝑓) –𝑉 𝑔, and the drag 𝐹_𝐷 = 1 𝐶_{𝐷 2}𝜌_𝑓|𝑉 |𝑉 𝐴. Here 𝜌_𝑓 is the density of the surrounding fluid (either 𝜌₁ or 𝜌₂ in the homogeneous layers), –𝑉 = 𝜋𝑎³∕6 is the particle’s volume, 𝐴 = 𝜋𝑎²∕4 is the projected surface area, and 𝐶_𝐷 is the drag coefficient.

The drag coefficient is estimated using the correlation (White, 1974; Verso et al., 2019):

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

In analogy to previous studies looking for correlations (Srdić-Mitrović et al., 1999; Abaid et al., 2004; Verso et al., 2019; Boetti and Verso, 2022), we also used two average properties: the average density of the fluid medium 𝜌, and the average of terminal velocities 𝑉 = 0.5(𝑉₁ + 𝑉₂). From these properties, we are able to calculate the expected ‘‘crossing time’’̂ 𝑡_{𝑉1−𝑉2} =̂ ℎ∕𝑉 that represents the time it takes the particle to crosŝ ℎ assuming that 𝐹_𝑠 is negligible. Finally, we calculated the delay times 𝜏 as the delta between the actual retention time 𝑡_{𝑉1−𝑉2} and the expected crossing timê 𝑡_{𝑉1−𝑉2} (0 ≤ 𝜏< ∞).

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 (̂𝑡_{𝑉1−𝑉 2}) that is greater than the measured retention time (𝑡_{𝑉1−𝑉 2}). By definition̂ 𝑡_{𝑉1−𝑉 2} ≤ 𝑡_{𝑉1−𝑉 2} leading to 0 ≤ 𝜏. However, since we do not measurê 𝑡_{𝑉1−𝑉 2} directly but estimate it using the averaged interfacial layer thicknesŝ ℎ, some particles appear as they havê 𝑡_{𝑉1−𝑉 2} > 𝑡_{𝑉1−𝑉 2}. 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 ≤ 𝜏< ∞.

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 𝑉₂ (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 𝑉₂ 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.

The dataset comprises of 9 dimensional properties of particles and fluids, i.e., particle diameter and density 𝜌_𝑝, 𝑎, the terminal velocities 𝑉₁, 𝑉₂, the delay time 𝜏, the fluids density and viscosity 𝜌₁, 𝜌₂, 𝜈, and the averaged interfacial layer thicknesŝ ℎ. Prior to symbolic regression of a correlation for 𝜏, we initially create a dataset of non-dimensional parameters.

According to the Buckingham theorem of dimensional analysis, the number of non-dimensional parameters 𝜋_𝑖 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 𝜋_𝑖 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 (𝑎,̂ ℎ) and three different velocities (𝑉₁, 𝑉₂, 𝑉 = 0.5(𝑉₁ + 𝑉₂)), 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 𝜏 defined above. For that, we used the expected crossing time 𝑡_{𝑉1−𝑉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̂ 𝜏 = 𝜏∕̂𝑡_{𝑉1−𝑉2}. Theoretically, it can be any value in the range 0 ≤̂ 𝜏< ∞, 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 𝜋_𝑖 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 𝜋_𝑖) 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̂ 𝜏 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 𝜋_𝑖 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 𝜋_𝑖. The subsets are evaluated and compared by training an AutoML model (Simon Keren et al., 2023) that predictŝ 𝜏 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̂ 𝜏 of all particles in the dataset. For that purpose, the Las Vegas - based SR component (Simon Keren et al., 2023) 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 𝛼 and a constant term 𝛽 that is optimized for our dataset of samples (i.e.,̂ 𝜏 = 𝛼𝑓 (𝜋₁, … , 𝜋₆) + 𝛽). For cases of data collected from particles or fluid medium that differ from this experiment, such as the data by Verso et al. (2019), we run another optimizer to obtain the 𝛼 and 𝛽 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. (2019).

Results

Feature selection

We report here the set of six 𝜋_𝑖 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 timê 𝜏.

In Table 3, we present a comparison between the most typical choice in the literature and the ML-selected form of each 𝜋_𝑖 parameter. In the right column, we present the value range for each parameter according to the selected form and our measurements.

In Fig. 7, we comparê 𝜏 = 𝑓 (𝜋_𝑖) in the conventional (top row) and the selected (bottom row) form of four different 𝜋_𝑖 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 𝜋_𝑖 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 𝑅𝑒-like and 𝐹 𝑟-like ratios, which from here on now, will be named̂ 𝑅𝑒 and̂ 𝐹 𝑟.

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. (2019) (marked as B, varying from type P1 to P4), over thê 𝐹 𝑟 = 𝑓 (̂𝑅𝑒) map, similarly to Fig. 2. We did not plot the measurements of Abaid et al. (2004) and Srdić-Mitrović et al. (1999) that were presented in Fig. 2, as we do not have access to the data needed to calculatê 𝐹 𝑟. Note that particle type P3 corresponds to the case that did not show minima (see Verso et al. (2019)). This new map emphasizes for the first time that in terms of the new dimensionless parameterŝ 𝐹 𝑟 and̂ 𝑅𝑒, all the particle types are different, a fact that was not observable in the conventional 𝐹 𝑟₁ = 𝑓(𝑅𝑒₁) map (Verso et al., 2019). Note that for our measurements, there is also a notable though the small difference between the particles without and with significant minima.

𝜋𝑖	No. of options	Conventional form	Selected form	Range
Length scales	1̂	ℎ∕𝑎̂	ℎ∕𝑎	25–50
Characteristic velocity	1	𝑉1∕𝑉2	𝑉1∕𝑉2	1–1.8
𝑅𝑒-like	6	𝑉1𝑎∕𝜈	𝑉𝑎∕𝜈	103–213
Characteristic density	4	𝜌2∕𝜌1	𝜌𝑝∕𝜌	2.20–2.53
Density jump	16	𝜌2∕(𝜌2 √ −𝜌1)	𝜌2∕(𝜌𝑝−𝜌) √	0.67–0.91
𝐹𝑟-like	384	𝑔̂ 𝜌2−𝜌1 𝑉1∕𝑎 ℎ 𝜌	𝑔 𝜌𝑝−𝜌2 𝑉1∕ℎ 𝑎 𝜌	0.02–0.06

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 coefficient of determination 𝑅², 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. (2019). 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.

In Fig. 9, we show the scatter of the predicted normalized delay time (using 𝑓₁ −𝑓₃) versus the actual measurements of̂ 𝜏. 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 𝑓₁ (left) is closest to the expected 𝑦 = 𝑥 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 Verso et al. (2019)), using two different fluid layer combinations, five particle types, and 10 density jumps. Based on them we suggest to use 𝑓₁ as the correlation for̂ 𝜏, which can be written as:

	Equation	Constants	𝑅2	MAE	MSE
𝑓1	𝛼1 (̂𝐹𝑟 𝜌𝑝 ) + 𝛽1 𝜌	𝛼1 = 16.19, 𝛽1 = −0.5	0.75	0.24	0.23
𝑓2	(̂𝐹𝑟̂𝑅𝑒ℎ ) 𝛼2 𝑎 + 𝛽2	𝛼2 = 2.76, 𝛽2 = 0.47	0.71	0.23	0.20
𝑓3	(̂𝐹𝑟 2 𝜌𝑝 𝜌𝑝 𝑉2 )) 𝛼3 𝜌(1 + 𝜌 𝑉1 + 𝛽3	𝛼3 = 89.91, 𝛽3 = 0.44	0.72	0.24	0.23

and if we substitutê 𝑁 it reads:

Discussion

In this study, we took a new approach to reveal the structure of empirical results. We utilized machine learning 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 trajectory datasets. 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.

Consistently with previous results (Srdić-Mitrović et al., 1999; Abaid et al., 2004; Mandel et al., 2020; Verso et al., 2019), we also find that not all particles in this parameter range experience slowdown and minimal velocity. In our case, 18% of all the trajectories did not experience a minima. More research is needed to understand why this effect occurs, but possible factors include particle rotation, the timing of particle release, non-spherical shapes, or interface oscillation. New experiments would be necessary to better understand this effect.

We demonstrate that the normalized delay time (̂𝜏) has a stronger dependence on the newly selected forms of dimensionless parameters rather than the conventional form, as shown in Fig. 7. We also inspect the selected forms of the dimensionless parameters and infer the possible meaning of the selection. Thus,̂ 𝑅𝑒 is based on average velocity and not on the entrance velocity as it was suggested by previous works (Srdić-Mitrović et al., 1999; Verso et al., 2019). 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 𝜌₂∕𝜌₁.

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, ℎ, and only as a square root of the particle diameter 𝑎^−0.5. In other words, we can infer that the dimensionless frequency,̂ 𝑁 should be defined with the particle diameter, while Froude number aŝ 𝐹 𝑟 = 𝑉 ∕̂ℎ̂ 𝑁. The better choice of this form for the Froude number is also supported by the better collapse of data than the rest of the 𝜋_𝑖 parameters. We also compared it to two Froude number forms and see that it is significantly stronger related to the delay time, i.e.̂ 𝜏 = 𝑓(̂𝐹 𝑟) leads to 𝑅² = 0.19, while both 𝐹 𝑟₁ and 𝐹 𝑟₂ lead to an order of magnitude weaker correlation with 𝑅² = 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 for 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. Additionally, we appreciate the idea of an anonymous reviewer to try to construct additional dimensionless features using a new length scale that instead of ℎ will be of the type of 𝑙_𝜌 = 2(𝜕𝜌∕𝜕𝑧)(𝜕²𝜌∕𝜕𝑧²). We encourage the readers to use SciMED’s open-source code and the openly shared data to verify whether this analysis leads to a better prediction of the retention time scale.

CRediT authorship contribution statement

Liron Simon Keren: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft.

Teddy Lazebnik: Software, Supervision, Visualization, Writing – review & editing. Alex Liberzon: Conceptualization, Data curation, Funding acquisition, Investigation, Project administration, Supervision, Validation, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

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.

Funding

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

Article notes

Publication history

Received 17 September 2023 · Accepted 22 December 2023 · Published 2 January 2024

References

• Abaid, N., Adalsteinsson, D., Agyapong, A., McLaughlin, R., 2004. An internal splash: Levitation of falling spheres in stratified fluids. Phys. Fluids 16 (5), 1567–1580. link
• Alahyari, A., Longmire, E., 1994. Particle image velocimetry in a variable density flow: application to a dynamically evolving microburst. Exp. Fluids 17 (6), 434–440. link
• Blanchette, F., Shapiro, A., 2012. Drops settling in sharp stratification with and without marangoni effects. Phys. Fluids 24 (4), 042104. link
• Boetti, M., Verso, L., 2022. Force on inertial particles crossing a two layer stratified turbulent/non-turbulent interface. Int. J. Multiph. Flow 104153. link
• Burd, A., Jackson, G., 2009. Particle aggregation. Ann. Rev. Mar. Sci. 1, 65–90. link
• 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 (2010), 436–465. link
• Doostmohammadi, A., Dabiri, S., Ardekani, A.M., 2014. A numerical study of the dynamics of a particle settling at moderate Reynolds numbers in a linearly stratified fluid. J. Fluid Mech. 750, 5–32. link
• Kok, J., 2011. A scaling theory for the size distribution of emitted dust aerosols suggests climate models underestimate the size of the global dust cycle. Proc. Natl. Acad. Sci. 108 (3), 1016–1021. link
• MacIntyre, S., Alldredge, A., Gotschalk, C., 1995. Accumulation of marine snow at density discontinuities in the water column. Limnol. Oceanogr. 40 (3), 449–468. link
• Magnaudet, J., Mercier, M., 2020. Particles, drops, and bubbles moving across sharp interfaces and stratified layers. Annu. Rev. Fluid Mech. 52, 61–91. link
• Mandel, T., Waldrop, L., Theillard, M., Kleckner, D., Khatri, S., et al., 2020. Retention of rising droplets in density stratification. Phys. Rev. Fluids 5 (12), 124803. link
• Mrokowska, M.M., 2018. Stratification-induced reorientation of disk settling through ambient density transition. Sci. Rep. 8 (1), 1–12. link
• OpenPTV consortium, 2014. Open source particle tracking velocimetry. link
• Prairie, J.C., Ziervogel, K., Arnosti, C., Camassa, R., Falcon, C., Khatri, S., McLaughlin, R.M., White, B.L., Yu, S., 2013. Delayed settling of marine snow at sharp density transitions driven by fluid entrainment and diffusion-limited retention. Mar. Ecol. Prog. Ser. 487, 185–200. link
• Renggli, C.J., Wiesmaier, S., De Campos, C.P., Hess, K.U., Dingwell, D.B., 2016. Magma mixing induced by particle settling. Contributions Mineral. Petrol. 171 (11), 96. link
• Simon Keren, L., Liberzon, A., Lazebnik, T., 2023. A computational framework for physics-informed symbolic regression with straightforward integration of domain knowledge. Sci. Rep. 13, 1249. link
• Smith, D., Simon, M., Alldredge, A., Azam, F., 1992. Intense hydrolytic enzyme activity on marine aggregates and implications for rapid particle dissolution. Nature 359 (6391), 139–142. link
• Srdić-Mitrović, A., Mohamed, N., Fernando, H., 1999. Gravitational settling of particles through density interfaces. J. Fluid Mech. 381, 175–198. link
• Turco, R., Toon, O., Ackerman, T., Pollack, J., Sagan, C., 1983. Nuclear winter: Global consequences of multiple nuclear explosions. Science 222 (4630), 1283–1292. link
• Verso, L., 2020. Particles Crossing Density Interfaces (Ph.D. thesis). Tel Aviv University. link
• Verso, L., van Reeuwijk, M., Liberzon, A., 2017. Steady state model and experiment for an oscillating grid turbulent two-layer stratified flow. Phys. Rev. Fluids 2, 104605. link
• Verso, L., van Reeuwijk, M., Liberzon, A., 2019. Transient stratification force on particles crossing a density interface. Int. J. Multiph. Flow 121, 103109. link
• Wang, S., Kandel, P., Deng, J., Caulfield, C., Dalziel, S., 2023, Bouncing behaviour of a particle settling through a density transition layer, https://dx/doi/10.48550/arXiv.2301.01484. link
• White, F., 1974. Viscous fluid flow mcgraw-hill inc. New York, New York. link
• Yick, K., Torres, C., Peacock, T., Stocker, R., 2009. Enhanced drag of a sphere settling in a stratified fluid at small Reynolds numbers. J. Fluid Mech. 632, 49–68. link