A new perspective for research on the mechanism and kinetic model of aggregation between coastal spilled oil and suspended sediment

Abstract

Understanding oil spill behavior in estuarine and coastal systems requires knowledge of oil–mineral aggregate (OMA) formation, as this process significantly governs the transport of floating oil. Current understanding of OMA formation remains superficial, lacking in-depth analysis of underlying micromechanisms. Moreover, existing models rely solely on mineral concentration, resulting in limited applicability. In this study, the aggregation between moderately dispersed oil and various minerals was investigated through mesoscale simulation experiments conducted in a wave tank. Physicochemical and morphological analyses confirmed the van der Waals, electrostatic, and impact forces between oil droplets and minerals. Their relative contributions under varying conditions were qualitatively assessed. Furthermore, a new oil-attenuation equation under particle intervention was proposed, and the coupling of the oil spill dispersion model and the OMA density prediction model was achieved. Based on the experimental data, the expression of the integrative coefficient (α) in relation to the key characteristics of minerals was derived. The new model can accurately predict the time-dependent oil sedimentation at high mineral concentrations. These results can offer technical support for maritime management and marine environmental protection departments to quantitatively evaluate the settlement degree and hazard scope of coastal oil spills.

---

Environmental significance

Understanding the aggregation mechanisms between spilled oil and suspended sediment in coastal environments is paramount for effective oil spill response and ecosystem protection. This research provides a crucial new perspective on the fundamental mechanisms and kinetics governing oil–mineral aggregate (OMA) formation. By developing a refined kinetic model, we significantly enhance predictive capabilities for OMA behavior-determining whether oil is transported, sinks to smother benthic habitats, or persists as a toxic submerged mat. The model can automatically input real-time data, such as hydrological parameters and sediment properties, and output an approximate proportion of oil spill sinking or prioritize response zones. Our findings address critical knowledge gaps in near-shore spill modeling, ultimately contributing to more targeted, effective mitigation strategies that minimize ecological damage.

1. Introduction

The escalation of offshore oil exploration activities and the concomitant rise in global crude oil transportation through maritime tanker operations have precipitated a series of oil spill incidents in recent decades. Compared to floating oil, sunken oil presents significantly greater challenges due to its distinct behavioral characteristics.^1,2 It typically forms when oil and sediment are co-suspended in the water column by breaking waves or the erosion of heavily oiled sand from beaches.³ Previous studies have shown that clays can act as dispersants, contributing to the natural recovery of oil-impacted shorelines and influencing the efficacy of purification processes like surf washing.^4,5 The agglomerates formed during this process are commonly termed “oil–mineral aggregates (OMAs)”.⁶ However, while a small fraction of spilled oil associated with OMAs may be biodegraded, most OMAs persist in the marine environment. This poses a long-term threat to benthic ecosystems through contaminant retention. OMA formation is governed by a multifaceted interplay of factors, encompassing three primary categories: the properties of the oil (viscosity, droplet size and composition), the properties of minerals (wettability, concentration and size), and the environmental characteristics such as temperature, mixing energy and salinity.^7–10 OMA formation is widely attributed to Pickering emulsion mechanisms, where the emulsions show remarkable resistance to coalescence and can maintain stability for several weeks in seawater environments, as supported by extensive experimental evidence.^11,12 However, this theory was formulated solely based on surface tension while neglecting other critical factors such as hydrodynamic interactions and electrostatic forces. In contrast, alternative mechanisms have been proposed that provide comprehensive explanations for most adhesion phenomena.^13–15 Although extensive work has elucidated fundamental mechanisms underlying OMA formation and properties, current predictive capacity remains limited. Zhao et al. proposed a numerical model called A-DROP to predict the formation of droplet OMA,¹⁶ defined as individual oil droplets surrounded by micron-sized particles. Nevertheless, certain model parameters were difficult to quantify under actual hydrodynamic conditions. Payne et al. established a second-order rate equation to describe the oil–particle interactions,¹⁷ but accounted solely for particle concentration and size while neglecting other critical factors.

To synthesize the research gaps mentioned above, this investigation systematically examines the critical role of particle characteristics for OMA formation. First, the relative contributions of van der Waals forces, electrostatic force and impact force during oil–mineral interaction across varying environmental conditions were differentiated at the atomic-scale and bulk-scale. On this basis, we described an enhanced predictive model of OMA formation based on Payne et al.¹⁷ that is extended to account for the behavior of oil dispersion, OMA density and particle property modulation effects. The model coefficients were then rigorously calibrated through multivariate regression analysis of experimental datasets complemented by our previous research. The present kinetic model incorporates a parameterization scheme that accounts for surface characteristics, thereby eliminating errors arising from the use of fixed constants. This dynamic adjustment enables the model to quantify temporal variations in aggregation rates induced by seawater quality. The results will establish a mechanistic framework for predicting sunken oil formation kinetics in nearshore environments, with particular relevance to coastal oil spill mitigation strategies.

2. Materials and methods

2.1 Experimental apparatus

In this study, a stainless steel wave tank was developed (1.5 m × 0.3 m × 0.5 m) with transparent observation panels. Specific dimensions and additional equipment are shown in Fig. 1. Seawater with a salinity of 28‰ was pumped from the Bohai Sea, China, and filtered through 0.45 mm pore-size filters before use. The system maintained a 20 cm hydrostatic depth while generating breaking wave conditions through a programmable piston-type wavemaker, achieving wave characteristics analogous to whitecap formation in marine environments. Wave height and frequency were measured using a digital wave gauge. A pneumatically controlled sampling array with a peristaltic pump was used to collect and transfer the suspension into the oil in a water monitor (Turner, E09 TD-4100XDC, American). This analytical approach eliminated solvent extraction requirements inherent to conventional UV-vis and IR spectrophotometry methodologies. An industrial camera was placed outside the tank to capture the photographs of wave patterns and dispersion of floating oil.

Fig. 1 Schematic diagram of the experimental setup and process.

2.2 Physicochemical properties of typical oil and test particles

In this study, Roncador crude oil originating from Brazil was provided by China Certification & Inspection Group. The density, kinematic viscosity, asphaltene content and sulfur content at 20 °C was 0.9433 g cm⁻³, 1526 mm² s⁻¹, 4.41% and 1.76%, respectively. Five types of particles including four minerals and natural sediment were used to form the OMA. Montmorillonite, kaolin, quartz sand and illite were purchased from Andy Metal Materials Co., Ltd (Nanjing, Jiangsu Province, China). The sediment sample was collected from the tidal flats of the Xiajiehe bathing beach (39°404.82″N and 121°51′33.99″W). Following an initial rinse with deionized water and subsequent oven drying at 60 °C for 24 h, the samples underwent particle size fractionation through sequential sieving to remove gravel (>2 mm) and fine particulate impurities. A multi-technique characterization approach was employed. Specific indicators and the corresponding test methods are shown in Table 1.

Table 1 Measured indicators and reference methods for particulate properties

Name of the test	Test method
Size distribution	Particle size analyzer (LISST-100X; Sequoia, Seattle, WA)
Hydrophobicity	Optical contact angle measuring and contour analysis (Dataphysics, Germany)
Cation exchange capacity (CEC)	Barium chloride–sulphuric acid exchange method
SEM characterization	Scanning electron microscopy (Nova Nano 450, FEI)
Surface area	BET nitrogen adsorption–desorption measurements (ZEISS, Sigma 300)

---

2.3 Experimental design and procedure

First, precisely weighed particles were initially homogenized with 500 mL of seawater in a 1 L glass beaker using magnetic stirring for 5 min at 1000 rpm. The resulting slurry was subsequently transferred to the wave tank and the wave maker was deployed and continuously operated for 10 min to ensure uniform particle distribution throughout the tank. The motor of the wavemaker operated at 200 rpm to generate a breaking wave with a frequency of 2.60 Hz, wave length of 45 cm and wave height of 3.89 cm. Following stabilization, 45 g of test oil was poured onto the water surface and experimental timing commenced immediately upon formation of a continuous oil slick, with each trial maintaining a 6-hour duration. After the wave maker stopped running, suspensions were collected from 12 sampling locations: four horizontal positions (15, 45, 65, and 95 cm downstream from the oil application point) at three discrete depths (3, 9, and 15 cm below the surface level). Oil dispersion efficiency (ODE) was used to assess the effect of minerals on oil dispersion as a whole, calculated usingwhere C_sample presents the average value of C_oil from the 12 sampling sites (mg L⁻¹), V_wt is the total effective water volume of the wave tank, and m_w is the mass of stranded oil on the side wall (g).

Oil sedimentation experiments were conducted separately. The procedures for introducing oil and particles were consistent with those described above. The mixing durations were set at 10, 20, 30, 60, 120, 240, and 360 min, respectively. Triplicate experiments were performed for each condition. After the wave maker stopped running, a 2-hour quiescent phase for oil/seawater stratification was followed, and the aqueous slurry was drained from the tank. Subsequently, the sunken OMAs were quantitatively recovered using a sterile plastic syringe and aliquoted equally into separatory funnels. The extraction protocol for petroleum hydrocarbons was implemented in strict accordance with the standardized methodology described by Gao et al.¹⁸

2.4 Analytical methods

2.4.1 Morphology of OMAs. The three-dimensional architecture and mineral permeability of OMAs were characterized using a confocal laser scanning microscopy system (Andor Technologies). High-resolution imaging was performed with a 20× oil-immersion objective using standardized acquisition parameters (100 ms of exposure time). To achieve optimal differentiation between particulate matter and oil droplets, dual-wavelength excitation was employed with 488 nm and 640 nm lasers through separate detection channels.

2.4.2 Density of OMAs. In this study, the effective density of individual OMAs was indirectly determined through measurements of the settling velocity. Static sedimentation experiments were conducted in a transparent acrylic column (80 × 80 × 150 mm) filled with deionized water. To enhance optical contrast during recording, a high-intensity LED light source was positioned behind the column. 3 mL aliquot of OMA suspension was carefully injected 1 cm below the water–air interface using a calibrated micropipette, ensuring minimal disturbance from surface tension effects. Temporal measurements were initiated when OMAs crossed the upper calibration line and terminated upon reaching the lower reference line. The settling velocity was derived from the vertical distance between demarcation lines divided by the transit interval, with subsequent density calculation performed using the improved Stokes formula:¹⁹where ρ_oma is the bulk density of OMAs, V_oma is the settling velocity, D_oma is the sphere-equivalent diameter of OMAs, g is the gravitational acceleration, and ρ_s and µ_s are the density and the kinematic viscosity of water, respectively.

2.4.3 Zeta potential measurements. To validate electrostatic interactions, zeta potentials of pure mineral suspensions and oil–mineral suspensions in deionized water and seawater were compared. For pure mineral suspensions, 0.1 g of mineral powder was dispersed into 100 mL of test solution and sonicated for 15 min to ensure homogeneity. For oil–mineral suspensions, 100 mL of water-accommodated fractions (WAFs) of Roncador crude oil was prepared and incubated at 25 °C for 30 min under gentle magnetic stirring to achieve thermodynamic equilibrium. Subsequently, 1 mL of each suspension was loaded into a capillary cell, and zeta potential measurements were performed at 25 °C using a zetasizer (Malvern Panalytical, Advance). The viscosity and dielectric constant were set to 0.95 cP and ε = 80, respectively. Each measurement was repeated in triplicate, and the results were reported as arithmetic means.

2.5 Data analysis

The experimental data were statistically analyzed using Welch's t-test. The statistical analysis was done using Origin Pro 7.5 (Origin Lab, Northampton, MA).

3. Results and discussion

3.1 Analysis of the interaction force between oil droplets and particulate matter

3.1.1 van der Waals forces between oil droplets and minerals particles. In this study, the driving force intensities governing OMA formation were first evaluated. SEM analysis revealed distinct morphological variations among the four minerals (Fig. S1). van der Waals forces between the oil droplets and four mineral particles were determined using computational chemistry methods. All calculations were performed at the B3LYP-D3(BJ)/6-311+G(2d,p) theory level with the Gaussian 09-E01 software package. Since n-C₁₆–n-C₂₄ are the dominant hydrocarbons in the n-alkane fraction of Roncador crude oil and naphthalene accounts for a higher proportion in the PAH fraction, naphthalene and pristane (C₁₉H₄₀) were selected as the representative components of the oil phase. For mineral components, Al₂O₃ and SiO₂ were chosen as model systems reflecting the predominant composition of the studied particles. As depicted in Fig. 2, distinct interaction patterns emerged between different molecular pairs. The naphthalene–Al₂O₃ system revealed significant attraction between aluminum atoms and the aromatic π-system, manifesting an interaction energy of −24.14 kcal mol⁻¹. This stronger van der Waals interaction compared to the naphthalene–SiO₂ system (−19.53 kcal mol⁻¹) can be attributed to enhanced electron cloud interactions. The π-electron-rich naphthalene ring interacts synergistically with the positively charged Al centers, while in the SiO₂ system, the interaction primarily occurs between silicon atoms and individual carbon centers. For pristane and Al₂O₃, the alkyl chain of pristane showed a certain degree of curling towards Al₂O₃, indicating the existence of van der Waals forces, and the interaction energy was −20.75 kcal mol⁻¹. The pristane–SiO₂ interaction was proved weaker (−14.82 kcal mol⁻¹), consistent with the reduced polarizability of the silica surface.

Fig. 2 Diagram of van der Waals force interaction between naphthalene and pristane with silica and aluminum oxide, respectively.

3.1.2 Electrostatic force between oil droplets and minerals particles. As shown in Fig. 3a, the four mineral particles exhibited comparable zeta potential in seawater (pH 8.1, salinity 20‰), whereas they demonstrated marked increases in deionized water (pH 7.2). This phenomenon arises from the progressive adsorption of soluble seawater cations by negatively charged particle surfaces, effectively neutralizing surface charges and reducing the zeta potential magnitude. Notably, the oil emulsion system demonstrated pronounced instability, with oil droplets reflecting persistent floating and coalescence behavior. The measured zeta potential of the oil emulsion in the absence of minerals reached merely −6.6 mV. For the oil–mineral mixture system, kaolin, illite, and quartz sand displayed comparable zeta potential between their individual suspensions and oil–particle mixtures in deionized water. However, montmorillonite–oil mixtures exhibited a substantially enhanced negative zeta potential (−41.3 mV), representing a 10 mV decrease compared to pure particle suspensions. This distinctive behavior likely originates from montmorillonite's capacity for hydration-driven structural expansion, enabling persistent oil droplet interaction even in cation-depleted environments. Fig. 3b reveals that hydration swelling equilibrium was achieved within 5 days for all minerals in deionized water. While kaolin, quartz sand, and illite displayed comparable swelling capacities (≈2.3 mL g⁻¹), montmorillonite demonstrated threefold greater expansion (≈6.9 mL g⁻¹). This exceptional swelling behavior of montmorillonite can be attributed to the formation of primary and secondary hydration shells around interlayer cations.²⁰ These hydration shells can be occupied by polar organic ions in the oil, which leads to further expansion of the interlayer and promotes penetration of non-polar hydrocarbon. In seawater, the zeta potentials of montmorillonite and quartz sand in oil–mineral mixtures measured −16.3 mV and −25.5 mV, respectively, representing reductions of 15.1 mV and 8.2 mV compared to single-particle suspensions (P < 0.05). This significant decrease suggests the existence of electrostatic attraction between both montmorillonite–oil and quartz sand–oil interfaces in the marine environment, and the electrostatic force between montmorillonite and dispersed oil was stronger. In contrast, kaolin and illite showed minimal zeta potential variation between composite systems and individual suspensions, maintaining relatively stable values around −30 mV. This stability implies weaker interfacial electrostatic attractions for these clay minerals.

Fig. 3 Interfacial charge modulation in particle–oil–water systems and mineral swelling kinetics: (a) zeta potential of four particle suspensions and their mixture with dispersed oil in different water environments, and (b) the variation of the swelling volume of the four minerals with time in deionized water.

3.1.3 Impact forces between oil droplets and mineral particles. In this study, the three-dimensional morphological characteristics of various OMAs were investigated to quantitatively assess the magnitude of impact force between oil droplets and mineral particles. To isolate the effect of electrostatic force, deionized water was used for controlled environment preparation. At high wave energy, montmorillonite showed unique interfacial interaction with oil droplets, manifesting both surface adhesion and internal embedment. Additionally, the size of the encapsulated oil phase in montmorillonite systems significantly exceeded that in the other three minerals (Fig. 4). In contrast, quartz sand, kaolin, and illite showed a less pronounced coating effect and weaker penetration ability with respect to oil droplets. Previous findings have indicated that the smaller the particle size, the deeper the particles are embedded within the oil droplet.²¹ This effectively accounts for the fact that the impact force between montmorillonite and oil droplets was stronger than that between kaolin and illite. Although montmorillonite and quartz sand demonstrated comparable size and measurable electrostatic interactions with oil droplets in deionized aqueous, these intermolecular forces were insufficient to facilitate particle encapsulation within oil droplets. Comparative analysis reveals that the impact forces between quartz sand and oil droplets were substantially diminished relative to those observed in montmorillonite systems. This mechanistic disparity is primarily attributable to their contrasting surface hydrophobicity. It has been shown that hydrophobic agglomeration phenomena are fundamentally governed by non-covalent hydrophobic interactions.²² Enhanced surface hydrophobicity promotes the repulsion of water molecules and affinity for polar oil constituents. Therefore, particle size and surface hydrophobicity cooperatively govern the magnitude of impact force.

Fig. 4 Laser confocal image and three-dimensional scanning structure of OMAs formed by different minerals at high wave energy (red represents particles, while light green and yellow represent oil droplets).

3.2 Oil dispersion under the action of different mineral particles

According to the literature,²³ the average concentration of suspended sediment in the surface layer of the Bohai Sea is about 20–100 mg L⁻¹. Thus the concentration of mineral particles in this study was set at 50 mg L⁻¹ to investigate the effects of montmorillonite, quartz sand, kaolin, and illite on the oil dispersion. As shown in Fig. 5, distinct spatial distribution patterns emerged in the oil dispersion system under the influence of different mineral particles. In the control group, the dispersed oil primarily accumulated within the horizontal plane (90–130 cm) at medium vertical depth, maintaining an average concentration range of 60–65 mg L⁻¹. Notably, addition of montmorillonite induced a significant horizontal migration of the dispersed oil towards the 70–100 cm region while maintaining comparable vertical positioning, accompanied by a 15–20% concentration elevation (≈75 mg L⁻¹). Conversely, after the addition of illite, kaolin and quartz sand, horizontal distribution remained consistent with the control. However, vertical penetration of dispersed oil increased to 12–15 cm depth, and the concentration ranged from 65–85 mg L⁻¹, 60–65 mg L⁻¹ and 65–75 mg L⁻¹, respectively.

Fig. 5 Spatial distribution of Roncador crude oil after 6 h-dispersion under the action of the four mineral particles.

3.3 Oil sedimentation under the action of different mineral particles

Upon visual inspection, the sediment showed a light brown coloration within 10 min of mineral–oil interaction initiation. The oil sinking efficiency (OSE) was defined as the mass percentage of settled oil relative to the total oil, representing the proportion effectively incorporated into OMAs during vibrational mixing. This parameter was calculated using the following equation:where m_s is the mass of oil trapped in the formed sunken OMAs, m_t is the total mass of oil initially added and m_w is the mass of oil adhered to the side wall of the wave tank.

Fig. 6 demonstrates that all mineral particulates induced logarithmic progression in OSE. At a C_par0 of 50 mg L⁻¹, the addition of montmorillonite, kaolin, quartz sand, and illite significantly enhanced the OSE during the 10–360 min oscillation period, with respective increases from 0.86%, 0.33%, 0.50% and 0.26% to 8.59%, 4.92%, 7.04% and 2.61%, respectively. As revealed by the comprehensive mechanistic analysis in Section 3.1, under identical environmental conditions, montmorillonite demonstrated significantly stronger interfacial interactions with oil droplets compared to the other three minerals, including van der Waals forces, electrostatic force and impact force. This enhanced interaction theoretically leads to the highest mass of trapped oil, the prediction that aligns with experimental measurements. Quartz sand exhibited the second strongest interaction force, with its interfacial parameters being slightly lower than those of montmorillonite, resulting in correspondingly smaller OSE that was consistently reflected in experimental data. For kaolin and illite, both exhibited relatively weaker and comparable electrostatic force. However, kaolin demonstrated weaker van der Waals interaction but stronger impact force compared to illite. The observed higher mass of trapped oil bound to kaolin suggested that impact force plays a more significant role than van der Waals forces in governing OMA formation.

Fig. 6 Variation of OSE over time for different mineral particles.

3.4 Assessment of OMA driving force across varying environmental conditions

To elucidate the relative contributions of van der Waals forces, electrostatic force and impact force, a comparative analysis was conducted between oil–natural sediment and oil–montmorillonite systems. Initial observation indicated that van der Waals interaction between oil and montmorillonite was theoretically stronger than that with natural sediment, primarily attributed to montmorillonite's substantially higher specific surface area (14.31 m² g⁻¹vs. 6.73 m² g⁻¹). Electrochemical characterization in seawater revealed zeta potential values of −32.7 mV for natural sediment suspensions versus −24.1 mV for oil–natural sediment mixtures, demonstrating measurable electrostatic attraction. Notably, the potential difference magnitude (8.6 mV) approximated that observed in quartz sand systems under identical conditions (Fig. 3a), yet remained lower than the 15.1 mV differential recorded for montmorillonite. This electrostatic disparity provides indirect evidence that the electrostatic attraction in oil–natural sediment systems is comparatively weaker than that in oil–montmorillonite systems. As illustrated in Fig. S2, microscopic observations revealed that the montmorillonite coating formed a marginally more compact surface layer on oil droplets compared to natural sediment. In contrast, natural sediment presented more pronounced oil-embedding characteristics, with certain particulates penetrating inside the oil droplets. These observations suggested that the impact force between natural sediment and oil droplets was stronger than that of montmorillonite.

As depicted in Fig. 7a, OMA analysis reveals a 0.47 g lower oil content in montmorillonite–oil aggregates compared to natural sediment–oil aggregates at high wave energy. The superior formation efficiency of OMAs in natural sediment systems containing dispersed oil primarily stems from the predominance of impact force over the comparatively weaker electrostatic and van der Waals interactions. This demonstrated that impact force constitutes the predominant contributor to OMA formation at high wave energy. Notably, the oil content differential diminished by 47% as wave energy transitions from high to medium intensity, suggesting a progressive decline in impact force dominance. Interestingly, this trend reversed at low wave energy where montmorillonite–oil aggregates exhibited 0.18 g higher oil retention than that of natural sediment. This shift implies a transition in governing mechanisms, from kinetic-driven collision force at a high energy state to surface chemistry-dominated interactions under quiescent conditions. To systematically investigate the critical role of cation-mediated electrostatic interactions in OMA formation, comparative analyses of the variation of trapped oil in OMAs formed by natural sediment and montmorillonite were conducted under decreasing seawater salinity conditions (Fig. 7b). The corresponding seawater samples were prepared by diluting natural seawater with deionized water to achieve the target salinity levels. Obviously, natural sediment consistently maintained higher trapped oil than montmorillonite across all tested salinity gradients. Quantitative analysis revealed incremental disparities in trapped oil between the two kinds of OMA, with the ratio of the difference to the oil mass of natural sediment–oil aggregates progressively increasing from 0.398 to 0.507 as salinity decreased from 30‰ to 0‰. This trend suggested diminished predominance of electrostatic force and van der Waals forces accompanied by heightened significance of impact force in OMA formation under low-salinity conditions. The observed behavior may be attributed to the reduced cation-mediated double-layer compression and suppressed particle hydration expansion at lower ionic strengths.²⁴ Furthermore, enhanced particle hydrophilicity may lead to surface pore constriction and impaired physical adsorption capacity for oil components.²⁵

Fig. 7 Oil content in OMAs formed by natural sediment and montmorillonite under different (a) wave energy and (b) salinity conditions.

3.5 Simulation of OMA formation in a nearshore environment

3.5.1 Predictive model for oil dispersion. Building upon the classical oil spill models, the critical parameters governing oil dispersion processes were determined: oil slick thickness, fluid density, viscosity, oil–water interfacial tension, and seawater density. Li et al.²⁶ demonstrated that the dimensionless oil entrainment flux (Q₀) can be fundamentally characterized through two pivotal dimensionless numbers: the Weber number (We) and Ohnesorge number (Oh) as follows:

Q₀ = aWe^bOh^c (4)

where a, b and c are coefficients to be determined. We defined We as the ratio of hydrodynamically induced disruptive momentum to restoring interfacial tension, and it can be expressed as a function of multiple parameters: seawater density (ρ_w), oil density (ρ_oil), gravity, wave height (H_s), oil–water interfacial tension (IFT_ow), and the Rayleigh–Taylor (R–T) instability maximum diameter (d₀):

Oh is the ratio of viscous force to interfacial tension force and is a function of the oil dynamic viscosity (µ_oil), ρ_oil, IFT_ow, and d₀:

Taking logarithms on both sides of eqn (4) and converting to a linear equation:

ln ODE = ln a + b ln We + c ln Oh (8)

Considering that We is only related to H_s in addition to the oil properties, the mass of oil droplets dispersed into the water column at 360 min was taken as the ODE at equilibrium (ODE_e). Based on wave gauge measurements, the corresponding We and Oh were calculated by combining the results of IFT_ow, µ_oil and ρ_oil at different water temperatures (Table S1), and finally the unknown coefficients a, b and c were solved by multiple linear regression:

ln ODE = −7.442 + 1.078 ln We − 1.101 ln Oh (9)

where a, b and c are 5.861 × 10⁻⁴, 1.078 and −1.101, respectively. The coefficient of determination (R_(adj)²) of the corrected model was 0.902, indicating a high goodness of fit. The residual independence test showed that the multiple linear regression model was free of autocorrelation. The results of the F-test and t-test showed that F = 97.599 > F_α(1,20) and |t| ≥ t_α/2, where p < 0.05, suggesting that the multiple regression model and the regression coefficients were statistically significant.

3.5.2 Predictive model for OMA density. Owing to the complex interplay among multiple factors influencing OMA density, it is difficult to establish a theoretical model and impossible to circumvent the influence of interrelationships among the factors. To overcome these methodological limitations, partial least squares regression (PLS) was used to establish a prediction model for OMA density. This analytical approach can reduce dimensionality while preserving explanatory power and account for covariate interactions through latent variable extraction.²⁷ According to the results of our previous research,²⁸ the main influencing factors of OMA density were oil density, oscillation time, particle density, particle concentration, particle size and energy dissipation rate. Thus the above six factors were selected as independent variables, and regression analysis was performed on 145 sets of experimental data. The data were first normalized by deviation standardization, so that the results were mapped to the interval [0,1], eliminating the effects brought by different scales. Afterwards, the dimensionless sequence was followed by regression analysis using Simca-S (Version 13.0, Umetri AB & Erisoft AB) software, with computational conditions based on the default values:

ρ_oma = −9550.61 + 9.70ρ_oil + 1.15 × 10⁻²t + 0.43ρ_par + 1874.07C_par0 + 1.07 × 10⁶D_par + 49 148.80ε (10)

where ρ_oma, ρ_oil and ρ_par are the density of OMAs, test oil and particles (kg m⁻³); t is oscillation time (s); C_par0 is the particle concentration; D_par is particle size and ε is the energy dissipation rate (m² s⁻³).

As shown in Table 2, the cumulative cross-validation coefficients for extracted components (Q_(cum)²) were much higher than 0.5, showing good predictive ability and stability. The cumulative variance of the dependent variable (R_Y(cum)²) was 0.876, indicating that the model has a high goodness of fit. In addition, the residuals were all distributed between 2 and −2, with an overall free dispersion and a slope close to 0 (Fig. S3). This suggested that there was no systematic error in the model. The one-sample Kolmogorov–Smirnov test showed that the model residuals obeyed a normal distribution with a mean of 1.13 × 10⁻¹³ and a standard deviation of 0.98 at the 0.05 significance level (Fig. S4). It can be seen that the model residuals were mainly from chance errors and the model was reliable for estimating OMA density.

Table 2 Statistical parameters of the OMA density prediction model

n	P	A	R X(cum) ^2	R Y(cum) ^2	Q (cum) ^2	RMSE
145	6	3	0.546	0.876	0.856	142.093

---

3.5.3 Solving for the particulate characteristic coefficients. Payne et al.¹⁷ derived a kinetic equation describing the decay of free oil droplets due to the interaction with suspended particulate matter (SPM):where C_oil and C_par are the mass concentrations of the dispersed oil and the SPM (kg m⁻³), respectively, ε is the turbulent energy dissipation (m² s⁻³), ν is the kinematic viscosity of the oil phase (m² s⁻¹), t is the mixing time (s), and α is an integrative coefficient (m³ kg⁻¹), which accounts for the shape, size, and sticking properties of the SPM.

Extensive research has demonstrated that particle concentration serves as the predominant factor influencing OMA formation.^6,8 To investigate this relationship, α values for specific particle sizes across varying natural sediment concentrations were initially determined. The formula conversion process is elaborated in the SI. As shown in Fig. 8, nonlinear regression analysis revealed that at a C_par0 of 100, 150, 200, and 300 mg L⁻¹, the corresponding α values were 2.45 × 10⁻⁶, 2.57 × 10⁻⁶, 2.54 × 10⁻⁶, and 2.57 × 10⁻⁶ m³ kg⁻¹, respectively. Statistical analysis indicated no significant correlation between α and C_par0 (p > 0.05), with mean values fluctuating around 2.50 × 10⁻⁶ m³ kg⁻¹. Subsequent analysis of the particle size effect demonstrated a distinct pattern. When the average size of natural sediment was 3.24, 8.82, 10.74, 17.11 and 35.64 µm, respectively, the corresponding α decreased progressively to 2.57 × 10⁻⁶, 2.01 × 10⁻⁶, 1.82 × 10⁻⁶, 1.18 × 10⁻⁶ and 1.06 × 10⁻⁶ m³ kg⁻¹, respectively. The associated R² values confirmed robust model fitting. These results revealed a significant negative correlation between D_par and α (p < 0.05).

Fig. 8 Least squares fitting method for solving α at different initial sediment concentrations and sizes.

Following analysis of α in relation to sediment concentration and size, investigation was extended to diverse mineral types to validate the generalizability of α. As shown in Fig. 9, the α values were 1.76 × 10⁻⁶, 9.80 × 10⁻⁷, 1.54 × 10⁻⁶ and 9.87 × 10⁻⁷ m³ kg⁻¹ for montmorillonite, kaolin, quartz sand and illite, respectively, at a standardized concentration of 300 mg L⁻¹. The observed α variability arises from distinct physicochemical properties among mineral types. As established in Section 3.1, particle size, cation exchange capacity and specific surface area regulated the impact force, electrostatic force and van der Waals forces, respectively. Therefore, a semi-empirical formula for α can be expressed as

where a, b, c and d represent the undetermined coefficients. To enable analytical solution, a logarithmic transformation was applied to both sides of the equation, converting it into linear form:

ln(α) = ln a + b ln CEC + c ln SBET − d ln D_par (13)

Fig. 9 Least squares fitting method for solving α of different mineral particles.

The unknown coefficients were determined through multiple linear regression analysis performed with SPSS 26.0 software. The derived results were presented as follows:

Statistical analysis showed satisfactory model performance across multiple validation metrics. The modified model achieved a coefficient of determination (R²) of 0.851, suggesting strong explanatory power of the predictor variables. The residual independence test showed that the multiple linear regression model has no autocorrelation phenomenon. Furthermore, the results of the F-test and t-test showed that F = 16.28 > F_α(1,5) and |t| ≥ t_α/2, with p < 0.05, indicating that the regression model and regression coefficients are statistically significant.

3.5.4 Model validation. By integrating the expressions for the ODE prediction model, the OMA density prediction model and the solved α, the OSE kinetic model was obtained. To further validate the applicability of the above model, it was utilized to predict the OSE of Roncador crude oil after interaction with natural sediment at different sediment concentrations and times at an energy dissipation rate of 9.7 × 10⁻³ m² s⁻³. As shown in Table 3, there was an acceptable deviation between the predicted and experimental values, with relative error below 10%.

Table 3 Analytical accuracy of the OSE kinetic model for the verification of samples

Sediment concentration (mg L^−1)	Time (h)	Experimental results	Model calculation results	Relative error
100	4	29.87%	27.26%	−8.73%
	6	32.22%	33.71%	4.61%
150	4	38.23%	35.50%	−7.13%
	6	42.18%	41.51%	−1.59%
200	4	44.73%	44.33%	−0.89%
	6	48.30%	50.91%	5.41%
300	4	51.64%	54.25%	5.06%
	6	56.79%	61.31%	7.95%

---

As shown in Table 4, it was further compared with the empirical formula proposed by Lehr et al.:²⁹

where q_sed is the attenuated mass of dispersed oil per unit time and volume of water (kg m⁻³ s⁻¹); K_s depends on the type and size of particulate matter, with the reported value typically ranging between 10 and 20, and the present study employed a median value of 15 for subsequent calculation; ε is the energy dissipation rate (m² s⁻³); ν_w is the kinematic viscosity of seawater (m² s⁻¹); C_oil and C_sed are the mass concentrations of oil droplets and suspended particles in water, respectively (kg m⁻³).

Table 4 Comparison of prediction accuracy between the classical model and OSE kinetic model in this study

Sediment concentration (mg L^−1)	q sed (10^−6 kg m^−3 s^−1)			Relative error
	Experimental results	Lehr et al.^29 (2002)	OSE kinetic model	Lehr et al.^29 (2002)	OSE kinetic model
100	1.67	1.28	1.40	−23.35%	−16.17%
150	1.98	1.44	2.18	−27.27%	10.10%
200	2.50	3.31	2.82	32.40%	12.80%
300	2.94	3.85	3.30	30.95%	12.24%

---

The comparative analysis revealed that the empirical model developed by Lehr et al.²⁹ demonstrated limited predictive capability for the OSE kinetics of Roncador crude oil, with relative errors exceeding 23% across all test conditions. In contrast, the proposed OSE kinetic model showed improved prediction accuracy, maintaining relative errors below 17%. In general, the collision and aggregation of oil–sediment remained non-equilibrium throughout 2 h, accompanied by substantial heterogeneity in OMA characteristics. Furthermore, at low sediment concentrations (<50 mg L⁻¹), the system likely exhibited oil-dominated phase conditions, resulting in stochastic OMA formation patterns. These findings collectively demonstrated that the developed model provides reliable predictions of equilibrium OMA distribution patterns in nearshore hydrodynamic regimes. However, the predictive capability for the settlement ratio of oil spills in short-term or low sediment concentration waters still requires further improvement. The characterization of fine-grained changes in oil–sediment interactions at low concentrations can be achieved using microfluidic chips, providing more accurate data for optimizing the predictive model.

4. Conclusions

Our results provide significant insights into the mechanisms of OMA formation. Experimental data confirmed that the interaction between oil and minerals is governed by the synergistic interplay of three interparticle forces: van der Waals forces, electrostatic force, and impact force. van der Waals forces mainly depend on the specific surface area of minerals. Electrostatic force is predominantly determined by surface charge characteristics, hydration-induced swelling properties, and cation exchange capacity of mineral constituents. Comparatively, impact force demonstrates significant dependence on hydrodynamic conditions, with secondary correlations with the hydrophobicity and size of particles. In high-energy mixing regimes, impact force emerges as the predominant mechanism governing OMA formation through enhanced collision efficiency. As salinity declines, the contribution of electrostatic force and van der Waals forces to the OMA formation diminishes, suggesting weakened colloidal stabilization under low ionic strength conditions. The quantification of these three forces under standardized dynamic conditions represents a promising focus for future research. In terms of model building, following determination of the particle characteristic coefficient (α) across varied conditions via least squares optimization, two key findings emerged. First, α exhibited no dependence on the initial particle concentration but demonstrated a monotonic decrease with increasing size. Second, multivariate regression analysis revealed α to be effectively characterized by three fundamental parameters: particle size, specific surface area, and cation exchange capacity, with the constant showing correlation with oil properties. Through integration of these relationships into our enhanced dispersion oil attenuation kinetics framework, coupled with oil dispersion and OMA density prediction algorithms, we successfully established a novel model for OMA formation kinetics. This advancement represents significant progression in sunken oil modeling, particularly through its mechanistic incorporation of particle–fluid interaction thermodynamics. Subsequent studies can introduce the oil characteristic coefficient and employ the same solution strategy to enhance the generalizability of existing models across different scenarios.

Author contributions

Wanran Li: formal analysis, software, validation, supervision. Yue Yu: investigation, methodology, writing – original draft. Deqi Xiong: resources, formal analysis, project administration. Zhixin Qi: data curation, methodology. Xiaoan He: investigation, data processing.

Conflicts of 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

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: morphology details of OMA and statistical analysis that support the findings of this study. See DOI: https://doi.org/10.1039/d5em00645g.

Acknowledgements

This research was supported by the National Natural Science Foundation of China (42306181) and Dalian Science Technology Talent Innovation Support Policy Implementation Plan (2023RQ096 and 2024RY039) and in part by the Natural Science Foundation Liaoning Province, China (2025-BS-0923 and 2025-BS-0956).

References

1. M. Jacketti, C. J. Beegle-Krause and J. D. Englehardt, A review on the sinking mechanisms for oil and successful response technologies, Mar. Pollut. Bull., 2020, 160, 111626.
2. P. S. Dalyander, J. W. Long, N. G. Plant and D. M. Thompson, Assessing mobility and redistribution patterns of sand and oil agglomerates in the surf zone, Mar. Pollut. Bull., 2014, 80, 200–209.
3. J. Michel and P. Bambach, A Response Guide for Sunken Oil Mats (SOMs): Formation, Behavior, Detection, and Recovery, J. Mar. Sci. Eng., 2022, 10, 368.
4. E. H. Owens and K. Lee, Interaction of oil and mineral fines on shorelines: review and assessment, Mar. Pollut. Bull., 2003, 47, 397–405.
5. R. Jézéquel, J. Receveur, T. Nedwed and S. Le Floch, Evaluation of the ability of calcite, bentonite and barite to enhance oil dispersion under arctic conditions, Mar. Pollut. Bull., 2018, 127, 626–636.
6. Y. Gong, X. Zhao, Z. Cai, S. E. O'Reilly, X. Hao and D. Zhao, A review of oil, dispersed oil and sediment interactions in the aquatic environment: influence on the fate, transport and remediation of oil spills, Mar. Pollut. Bull., 2014, 79, 16–33.
7. X. Sang, Q. Miao, M. Bao, H. Li, D. Yan and P. Sun, The interaction between dispersed crude oil droplets and particulate matter, Environ. Sci.: Processes Impacts, 2020, 22, 1397–1407.
8. J. Sun, A. Khelifa, C. Zhao, D. Zhao and Z. Wang, Laboratory investigation of oil-suspended particulate matter aggregation under different mixing conditions, Sci. Total Environ., 2014, 473–474, 742–749.
9. S. Le Floch, J. Guyomarch, F. X. Merlin, P. Stoffyn-Egli, J. Dixon and K. Lee, The influence of salinity on oil-mineral aggregate formation, Spill Sci. Technol. Bull., 2002, 8, 65–71.
10. A. Khelifa, P. Stoffyn-Egli, P. S. Hill and K. Lee, Characteristics of oil droplets stabilized by mineral particles: Effects of oil type and temperature, Spill Sci. Technol. Bull., 2002, 8, 19–30.
11. S. M. Dieng, N. Anton, P. Bouriat, O. Thioune, P. M. Sy, N. Massaddeq, S. Enharrar, M. Diarra and T. Vandamme, Pickering nano-emulsions stabilized by solid lipid nanoparticles as a temperature sensitive drug delivery system, Soft Matter, 2019, 15, 8164–8174.
12. A. T. Tyowua, S. G. Yiase and B. P. Binks, Double oil-in-oil-in-oil emulsions stabilised solely by particles, J. Colloid Interface Sci., 2017, 488, 127–134.
13. D. Boglaienko and B. Tansel, Classification of oil-particle interactions in aqueous environments: Aggregate types depending on state of oil and particle characteristics, Mar. Pollut. Bull., 2018, 133, 693–700.
14. L. Zhao, M. C. Boufadel, J. Katz, G. Haspel, K. Lee, T. King and B. Robinson, A New Mechanism of Sediment Attachment to Oil in Turbulent Flows: Projectile Particles, Environ. Sci. Technol., 2017, 51, 11020–11028.
15. X. Zhong, H. Niu, P. Li, Y. Wu and L. Liu, An Overview of Oil-Mineral-Aggregate Formation, Settling, and Transport Processes in Marine Oil Spill Models, J. Mar. Sci. Eng., 2022, 10, 610.
16. L. Zhao, M. C. Boufadel, X. Geng, K. Lee, T. King, B. Robinson and F. Fitzpatrick, A-DROP: A predictive model for the formation of oil particle aggregates (OPAs), Mar. Pollut. Bull., 2016, 106, 245–259.
17. J. R. Payne, J. R. Clayton and B. E. Kirstein, Oil/suspended particulate material interactions and sedimentation, Spill Sci. Technol. Bull., 2003, 8, 201–221.
18. Y. Gao, X. Zhao, Z. Ju, Y. Yu, Z. Qi and D. Xiong, Effects of the suspended sediment concentration and oil type on the formation of sunken and suspended oils in the Bohai Sea, Environ. Sci.: Processes Impacts, 2018, 20, 1404–1413.
19. L. Ye, A. J. Manning and T. J. Hsu, Oil-mineral flocculation and settling velocity in saline water, Water Res., 2020, 173, 115569.
20. H. Yan and Z. Zhang, Effect and mechanism of cation species on the gel properties of montmorillonite, Colloids Surf., A, 2020, 611, 125824.
21. W. Ji, M. C. Boufadel, L. Zhao, B. Robinson, T. King and K. Lee, Formation of oil-particle aggregates: Particle penetration and impact of particle properties and particle-to-oil concentration ratios, Sci. Total Environ., 2021, 760, 144047.
22. W. Cheng, Z. Deng, X. Tong and T. Lu, Hydrophobic agglomeration of fine pyrite particles induced by flotation reagents, Minerals, 2020, 10, 801.
23. G. Zhao, W. Jiang, T. Wang, S. Chen and C. Bian, Decadal variation and regulation mechanisms of the suspended sediment concentration in the Bohai Sea, China, J. Geophys. Res.: Oceans, 2022, 127, e2021JC017699.
24. K. Yotsuji, Y. Tachi, H. Sakuma and K. Kawamura, Effect of interlayer cations on montmorillonite swelling: Comparison between molecular dynamic simulations and experiments, Appl. Clay Sci., 2021, 204, 106034.
25. C. You, C. Jia and G. Pan, Effect of salinity and sediment characteristics on the sorption and desorption of perfluorooctane sulfonate at sediment-water interface, Environ. Pollut., 2010, 158, 1343–1347.
26. Z. Li, M. L. Spaulding and D. French-McCay, An algorithm for modeling entrainment and naturally and chemically dispersed oil droplet size distribution under surface breaking wave conditions, Mar. Pollut. Bull., 2017, 119, 145–152.
27. A. R. L. Ferreira, L. F. Sanches Fernandes, R. M. V. Cortes and F. A. L. Pacheco, Assessing anthropogenic impacts on riverine ecosystems using nested partial least squares regression, Sci. Total Environ., 2017, 583, 466–477.
28. Y. Yu, Z. Qi, D. Xiong, W. Li, X. Yu and R. Sun, Experimental investigations on the vertical distribution and properties of oil-mineral aggregates (OMAs) formed by different clay minerals, J. Environ. Manage., 2022, 311, 114844.
29. W. Lehr, R. Jones, M. Evans, D. Simecek-Beatty and R. Overstreet, Revisions of the ADIOS oil spill model, Environ. Modell. Software, 2002, 17, 189–197.

---