Migration Movements of Accidentally Spilled Oil in Environmental Waters: A Review

Abstract

Accidentally spilled oil can cause great harm to the ecological balance of water once it enters the environmental waters. Clarifying its movement behavior and migration law in water has been the focus of environmental hydraulics research. This review starts from the mechanism of the oil spill migration process, and firstly reviews the kinematic characteristics of the smallest moving unit of the oil spill, the individual oil droplet, as well as focusing on several key aspects such as droplet shape, trajectory, terminal velocity and drag coefficient. Subsequently, considering the commonalities and differences between inland riverine and oceanic environments, different aspects of oil droplet collision, coalescence, breakage, particle size distribution, and vertical diffusion are discussed separately. Finally, the current status of research on the migration laws of accidental oil spills in environmental waters is summarized, and feasible future research directions are proposed to address the emerging research problems and research gaps.

1. Introduction

Oil as a transit fuel has become the world’s foremost source of energy, which has a non-negligible influence on the development of a country or even the entire global economy [1,2]. The UK, for example, requires petroleum products to meet 95% of its energy requirements in the transport sector [3]. With the growth in shipping activities and need for petroleum products such as fuel oil, the incidence of oil spills is also increasing. The level of oil in ocean, rivers, lakes and other water environments keeps rising with a serious threat to the surrounding water ecosystem [4,5,6]. Meanwhile, large oil spill accidents that involve ships or drilling rigs, etc., may cause damage to the aquatic environment for years or even decades [7]. For instance, Kuwait spilled approximately 20 million barrels of crude oil in 1991 as a result of the war, many recalcitrant contaminants (unresolved complex mixture) that require urgent degradation are still being found recently in the groundwater of the contaminated areas, and the degradation process will continue for a long time to come [8]. The Deepwater Horizon (DWH) accident occurred on 20 April 2010 in the Northern Gulf of Mexico, which spilled approximately 200 million gallons of crude oil and contaminated more than 20,000 square kilometers of sea area; the marsh soil strength is still not restored and the impact on the local biological community will last for at least 50 years [9,10]. The MV Wakashio grounding incident that occurred in the Mauritius wildlife sanctuary on 25 July 2020 resulted in the spill of 1000 metric tons of very low sulfur fuel oil into the Indian Ocean, causing varying degrees of harm to nearby coral reefs, the Île aux Aigrettes Nature Reserve, Mauritius dolphins, etc., and its ecological damage and impact on economic development will take a long time to repair [11].

Considering the enormous impact of oil spill accidents, not only oil-related projects such as wharf construction and drilling platforms need to assess the risk of oil spill in advance, but also emergency water quality prediction and early warning work must be carried out after sudden oil spill accidents [12], i.e., the emergency prediction of the spreading trend and migration paths of accidental oil spills, so as to provide technical guarantees for the emergency management of oil spills [13,14,15]. Therefore, it is crucial to grasp the physical mechanisms of the migration and transformation processes once the oil spill enters into environmental waters and to accurately determine the relevant parameters [16]. Research into the migration and transformation of oil in environmental waters has been conducted for many years and covers various aspects of ecosystems such as water bodies, sediments and aquatic organisms, which are summarized in Figure 1. Current studies have shown that, whether an oil spill occurs on the surface of the water body or as a submerged oil spill in the form of a jet stream or plume, the spilled oil will ultimately arrive at two destinations. A portion of the oil will float on the water surface as an oil film and move along with the water flow. The other portion of the oil will be retained as droplets of varying sizes for a long time in water bodies and form a certain spatial distribution (disregarding the portion that evaporates into the air and enters the sediment) [17,18,19,20,21]. Among these, the surface oil film has been focused on by most researchers due to its practical applicability that can directly show the distance and range of oil spill impacts, while there are relatively few studies on the migration process of oil completely immersed in the water body [22]. Although there are some previous reviews on the basic physical mechanisms of insoluble particles moving in quiescent water, these reviews have mainly focused on deformable fluid particles, i.e., simply considering bubbles and droplets in combination, which has many limitations in practical applications and must be verified in detail on a case-by-case basis before use [23,24]. Recent reviews have mostly focused on the construction of oil spill models, with only a few dealing with a single aspect of oil droplet aggregation, particle size distribution, and droplet–sediment interactions [5,19,25,26,27]. In conclusion, it is clear that there has been a persistent lack of reviews on the systematic study of the characteristics and fundamentals of the migration and transformation of fully submerged oils in environmental waters.

This review will commence with single oil droplet motions, gradually expanding to double oil droplet motions and the diffusion processes of oil droplet populations in different environments. This review is organized as follows: Section 2 provides a comprehensive review of the kinematic/dynamic characteristics and principles of the free-rising motion of a single oil droplet, the smallest unit of motion in an oil spill. Section 3 briefly describes the mutual motions between oil droplets of an oil spill in the inland riverine environment (droplet collision, coalescence, and breakage behaviors), as well as the study of the vertical dispersion characteristics of the oil spill under the action of water currents. Section 4 reviews the current state of research on the particle size distribution and vertical mixing characteristics of oil droplet populations in the oceanic environment. Section 5 summarizes the above contents and suggests future research directions. This review aims to summarize the present research situation on the migration process of accidentally spilled oil after submerging in environmental waters, discusses the bottlenecks in and prospects for the research, and hopes to provide some guidance for future research work.

2. Kinematics and Dynamics Characterizing of the Free-Rising Motion of Oil Droplets

The migration and transformation behavior of oil in water environments is the result of a complex combination of processes that occur simultaneously at the level of a single oil droplet. In order to clarify the macroscopic state of oil migration and transformation under different environmental conditions, it is necessary to have a clear and accurate perception of events at the microscopic scale [28]. Consequently, a single oil droplet, as the smallest microscopic scale of existence in the oil movement behavior, and its free-rising movement process in the water column become the basis of all the studies. At present, the studies on the free-rising motion of single oil droplets in water columns are generally carried out from the aspects of kinematics and dynamics, which involve the shape, trajectory, terminal velocity, and drag coefficient of the oil droplets [29,30,31]. Of these, the oil droplet terminal velocity and drag coefficient are crucial for characterizing the fluid dynamics of oil–water systems [31]. This is because both of them are the comprehensive manifestation of the joint influence of oil droplet particle size, shape, density, surface tension and other factors, reflecting the resistance to the water when the oil droplet interacts with the water body, which is one of the most vital physical quantities in the kinematic properties of oil. It can not only be directly applied to practical industrial problems such as the design of oil separator tanks in sewage treatment and oil–water separation devices in oil extraction, but can also be used immediately in research on the migration behaviors of oil spills in inland rivers or marine environments, which is a crucial technical parameter for the construction of oil spill models [32].

2.1. The Shape and Trajectory of Oil Droplets

The initial shape of oil droplets is primarily determined by their physical properties, mainly including particle size and interfacial tension. Once the oil droplets start moving, they will be deformed by the external flow field. Eventually, the shape becomes stable when the normal and shear stresses at the fluid–fluid interface balance each other [29,33]. Therefore, when compared with the infinite number of shapes possible for rigid solid particles, oil droplets at steady state are severely limited in the number of possibilities since such features as sharp corners or protuberances are precluded by the interfacial force balance [24]. Although complex shapes such as with and without skirts, crescents, biconcave discs, and annulars have been found by previous researchers in related studies, these shapes generally exist only in droplets resembling bubbles with very low surface tension and viscosity ratios, and are rarely found in other systems, especially oil–water systems [24,34,35]. The more currently accepted steady-state oil droplet shapes with increasing particle size are mainly classified into three categories: spherical, ellipsoidal, and spherical- or ellipsoidal-cap, which are illustrated in Figure 2. The specific definitions and delineation criteria basically use the generalized graphical correlation diagrams developed by Grace et al. which are based on the Eötvös number Eo, the Reynolds number Re, and the Morton number Mo [24,29,36,37]. When the oil droplet size is very large (usually the shape is characterized by a spherical-cap or an ellipsoidal-cap), the oil droplet exhibits strong self-oscillation and it is difficult to maintain a steady-state shape motion [38,39,40], but such oscillation will weaken with the decrease in the purity of the oil–water system [41].

Concurrently, related work on quantifying the degree of deformation of oil droplets and attempting to establish a relationship between them and oil properties is also underway. The parameters used to characterize the degree of oil droplet deformation vary from researcher to researcher. However, these parameters are basically benchmarked against a standard sphere, and then used to compare and establish the relationship between the horizontal (long-axis) and vertical (short-axis) axes of the oil droplet [42]. The qualitative relationship among the degree of oil droplet deformation and partial properties of oil is now largely and widely recognized. It is generally recognized that the smaller the density difference between oil and water, the smaller the particle size, the higher the interfacial tension of the oil droplet, the better the rigidity and the lower the degree of deformation [37]. However, disagreements always exist about the effect of oil viscosity on the degree of oil droplet deformation. Taylor et al. concluded that the degree of oil droplet deformation is not sensitive to the viscosity ratio κ between the two phases [43]. Myint et al. instead stated that the effect of viscosity ratio on oil droplet shape cannot be completely ignored. For oil droplets with κ > 1, the front of the droplet will be flatter than the back. Conversely, for oil droplets with κ < 1, the back will be flatter than the front [36].

Regarding the trajectory of the free-rising oil droplet, it is also essentially a reflection of the degree of deformation of the oil droplet [29,43]. While the oil droplet is approximately spherical, the oil droplet basically rises in a straight line, and when the oil droplet is deformed, the trajectory begins to show slight deviations and trajectory oscillations—even a spiral can occur [22,28,29,44]. Further analysis of the oscillation amplitude and oscillation frequency of oil droplets reveals that the equivalent diameter of oil droplets is positively correlated with the amplitude of the motion trajectory as well as the oscillation frequency [29]. In practical applications, the water flow velocity is much larger than the amplitude of transverse oscillation of small oil droplets, which generally can be approximately ignored, while large oil droplets with larger oscillations will break into small oil droplets because of water turbulence, oil droplet collision and other reasons [45]. In consequence, there are relatively few studies dedicated to the trajectory of oil droplets, which are generally considered as a verification of the degree of deformation of oil droplets, with the linear motion performed by default except in special cases.

2.2. The Terminal Velocity of Oil Droplets

Although the variable speed motion of oil droplets commonly exists, it generally only lasts for a few seconds or even a fraction of a second, so most studies focus on the uniform motion phase of oil droplets, when the speed is the terminal velocity of oil droplets [46,47]. The oil droplets, especially large oil droplets, however, are always fluctuating in instantaneous velocity due to the oscillation of their own shape, and cannot achieve velocity stability in the theoretical sense [48]. At present, it is generally considered that when the fluctuation of the vertical motion velocity of the oil droplet is much smaller than the rise velocity of the oil droplet or the regular periodic fluctuation, the oil droplet will enter the uniform motion phase, where the average value of the instantaneous velocity in a long enough uniform motion time period is the terminal velocity of the oil droplet.

Table 1. The formula for calculating the terminal velocity of oil droplets.

Time, Proposer	Expression	Applicable Situation		Eq.
1850, Stokes	${\mathrm{U}}_{\mathrm{T}}=\frac{\left({\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}\right)\mathrm{g}{\mathrm{d}}^2}{18{\mathsf{\mu }}_{\mathrm{c}}}$	$\mathrm{R}\mathrm{e}\ll 1$		(1)
1911, Hadamard and Rybczynski	${\mathrm{U}}_{\mathrm{T}}=\frac{\left({\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}\right)\mathrm{g}{\mathrm{d}}^2}{6{\mathsf{\mu }}_{\mathrm{c}}}(\frac{1+\mathsf{\kappa }}{2+3\mathsf{\kappa }})$	$\mathrm{R}\mathrm{e}\ll 1$		(2)
* 1913, Boussinesq [24]	${\mathrm{U}}_{\mathrm{T}}=\frac{\left({\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}\right)\mathrm{g}{\mathrm{d}}^2}{6{\mathsf{\mu }}_{\mathrm{c}}}(\frac{1+\mathsf{\kappa }+\raisebox{1ex}{$\mathrm{C}$}\!\left/ \!\raisebox{-1ex}{${\mathsf{\mu }}_{\mathrm{c}}$}\right.}{2+3\mathsf{\kappa }+\raisebox{1ex}{$3\mathrm{C}$}\!\left/ \!\raisebox{-1ex}{${\mathsf{\mu }}_{\mathrm{c}}$}\right.})$	$\mathrm{R}\mathrm{e}\ll 1$		(3)
1956, Klee and Treybal [49]	${\mathrm{U}}_{\mathrm{T}}=38.3{\mathsf{\rho }}_{\mathrm{c}}^{-0.45}{∆\mathsf{\rho }}^{0.58}{\mathsf{\mu }}_{\mathrm{c}}^{-0.11}{\mathrm{d}}^{0.70}$	Critical diameter d: $\mathrm{d}=0.33{\mathsf{\rho }}_{\mathrm{c}}^{-0.14}{∆\mathsf{\rho }}^{-0.43}{\mathsf{\mu }}_{\mathrm{c}}^{0.30}{\mathsf{\sigma }}^{0.248}$		(4a)
	${\mathrm{U}}_{\mathrm{T}}=17.6{\mathsf{\rho }}_{\mathrm{c}}^{-0.55}{∆\mathsf{\rho }}^{0.28}{\mathsf{\mu }}_{\mathrm{c}}^{0.1}{\mathsf{\sigma }}^{0.18}$			(4b)
* 1961, Devies and Rideal [24,50]	${\mathrm{U}}_{\mathrm{T}}=\frac{\left({\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}\right)\mathrm{g}{\mathrm{d}}^2}{18{\mathsf{\mu }}_{\mathrm{c}}}(1+\frac{\mathrm{Z}}{2})$ where: $\mathrm{Z}=\left[2/(2+3\mathsf{\kappa })\right]\left\{2(\mathrm{Y}-1)\right\}$	Presence of surface contaminants; $\mathrm{R}\mathrm{e}\ll 1$		(5)
1968, Thorsen et al. [38]	${\mathrm{U}}_{\mathrm{T}}=\frac{6.8}{1.65-\frac{∆\mathsf{\rho }}{{\mathsf{\rho }}_{\mathrm{d}}}}\sqrt{\frac{\mathsf{\sigma }}{3{\mathsf{\rho }}_{\mathrm{d}}+2{\mathsf{\rho }}_{\mathrm{c}}}}/\sqrt{\mathrm{d}}$	$400<\mathrm{R}\mathrm{e}<900$		(6)
1982, Aravamudan et al. [51]; 2018, Röhrs et al. [52]	${\mathrm{U}}_{\mathrm{T}}=\frac{\left({\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}\right)\mathrm{g}{\mathrm{d}}^2}{18{\mathsf{\mu }}_{\mathrm{c}}}$	Critical diameter $\mathrm{d}$: $\mathrm{d}=\frac{9.52{\mathsf{\mu }}_{\mathrm{c}}^{2/3}}{{(\mathrm{g}\frac{{\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}}{{\mathsf{\rho }}_{\mathrm{c}}})}^{1/3}}$	$\mathrm{R}\mathrm{e}\le 50$	(7a)
	${\mathrm{U}}_{\mathrm{T}}={(\frac{8}{3}\mathrm{g}\mathrm{d}(\frac{{\mathsf{\rho }}_{\mathrm{c}}-{\mathsf{\rho }}_{\mathrm{d}}}{{\mathsf{\rho }}_{\mathrm{c}}}))}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$		$\mathrm{R}\mathrm{e}>50$	(7b)
* 1978, Grace et al. [53]; 1978, Clift et al. [24]; 2000, Zheng and Yapa [17]	${\mathrm{U}}_{\mathrm{T}}=\frac{\mathrm{R}\mathrm{e}{\mathsf{\mu }}_{\mathrm{c}}}{\mathsf{\rho }\mathrm{d}}$ (In this formula, Re has a separate correlation, see 1978, Clift et al.)	Spherical (small size range)	$\mathrm{d}<1 \mathrm{m}\mathrm{m}$	(8a)
	${\mathrm{U}}_{\mathrm{T}}=\frac{{\mathsf{\mu }}_{\mathrm{c}}}{\mathsf{\rho }\mathrm{d}}{\mathrm{M}\mathrm{o}}^{-0.149}(\mathrm{J}-0.857)$ where $$\begin{gathered} \mathrm{J}=0.94{\mathrm{H}}^{0.757},(2<\mathrm{H}\le 59.3) \\ \mathrm{J}=3.42{\mathrm{H}}^{0.441},(\mathrm{H}>59.3) \\ \mathrm{H}=\frac{4}{3}\mathrm{E}\mathrm{o}{\mathrm{M}\mathrm{o}}^{-0.149}{(\frac{{\mathsf{\mu }}_{\mathrm{c}}}{{\mathsf{\mu }}_{\mathrm{w}}})}^{-0.14} \end{gathered}$$	Ellipsoidal (intermediate size range)	Critical diameter d $\mathrm{l}\mathrm{o}\mathrm{g}\mathrm{d}=\frac{{\mathrm{b}}_2-{\mathrm{b}}_1}{{\mathrm{a}}_1-{\mathrm{a}}_2}$ (For the values of ${\mathrm{a}}_1,{\mathrm{a}}_2,{\mathrm{b}}_1$ and ${\mathrm{b}}_2$, see Zheng and Yapa in 2000)	(8b)
	${\mathrm{U}}_{\mathrm{T}}=0.071\sqrt{\mathrm{g}\mathrm{d}∆\mathsf{\rho }/{\mathsf{\rho }}_{\mathrm{c}}}$	Spherical-cap (large size range), $\mathrm{E}\mathrm{o}>40$		(8c)
2007, Kelbaliyev and Ceylan [54]	${\mathrm{U}}_{\mathrm{T}}=\frac{22{\mathrm{d}}^2{\mathrm{M}\mathrm{o}}^{-1/4}}{1+{\mathrm{k}}_0{\mathrm{d}}^{3/4}+{\mathrm{k}}_1{\mathrm{d}}^{7/3}}$ where $$\begin{gathered} {\mathrm{k}}_0=0.2\left(0.1+{\mathrm{M}\mathrm{o}}^{-1/5}\right){\mathrm{k}}_1= \\ {(0.02+5{\mathrm{M}\mathrm{o}}^2)}^{-1} \end{gathered}$$	$9\times {10}^{-7}\le \mathrm{M}\mathrm{o}\le 78$		(9)
1976, Grace et al. [53]; 1978, Clift et al. [24]; 2010, Wegener and Kraume [31]	${\mathrm{U}}_{\mathrm{T},\mathrm{p}\mathrm{u}\mathrm{r}\mathrm{e}}={\mathrm{U}}_{\mathrm{T}}(1+\frac{\mathsf{\Gamma }}{1+\mathsf{\kappa }})$ where $\mathsf{\Gamma }=\mathsf{\alpha }·{\mathrm{e}}^{\frac{-{(\log \mathsf{\Lambda }-\mathsf{\chi })}^2}{2{\mathsf{\delta }}^2}}$ $\mathsf{\Lambda }=\mathrm{E}\mathrm{o}\frac{1+0.15\mathsf{\kappa }}{1+\mathsf{\kappa }}$	${\mathrm{U}}_{\mathrm{T}}$ used the formula proposed by Clift et al. in 1978		(10)
* 2017, Ervik and Bjørklund [55]	$\frac{{\mathrm{U}}_{\mathrm{T}}}{{\mathrm{U}}_{\mathrm{H}\mathrm{S}}}(\mathrm{x})=1$	Denotes the drop radius normalized by the critical radius as $\mathrm{x}=\mathrm{a}/{\mathrm{R}}_{\mathrm{c}}$	$\mathrm{x}\le 1$	(11a)
	$\frac{{\mathrm{U}}_{\mathrm{T}}}{{\mathrm{U}}_{\mathrm{H}\mathrm{S}}}(\mathrm{x})=\frac{3(\mathsf{\kappa }+1)}{3\mathsf{\kappa }+2}-\frac{\sqrt{2}}{3\mathsf{\kappa }+2}\frac{{\mathrm{x}}^{-2}}{\sqrt{1+{\mathrm{x}}^{-2}}}$		$\mathrm{x}>1$	(11b)

Note *: C is equal to the surface dilational viscosity divided by 1.5 times the radius; Z is degree of circulation (see specific paper for details); ${\mathsf{\mu }}_{\mathrm{w}}$ is the viscosity of water in Braida’s experiments, which may be taken as 0.0009 kg/ms (0.9 cp); ${\mathrm{R}}_{\mathrm{c}}=\mathsf{\beta }/{\mathrm{F}}_{\mathrm{\infty }}$(see specific paper for details).

summarizes some of the formulas for calculating the terminal velocity of oil droplets, with theoretically derived formulas, semi-empirical formulas and empirical formulas; the calculation methods involve the single formula method applicable to part of the Re zone and the segmental formula method is applicable to the whole subcritical zone.

A most significant theoretical analysis result when studying the terminal velocity of oil droplets is the H-R analytical solution (Equation (2)), which is obtained by considering the effect of the oil droplet’s own viscosity (i.e., internal circulation is considered to exist in the oil droplet) on the basis of Stokes solution, mainly for the motion of an ideal spherical oil droplet in the case of creeping flow (Re << 1). Although the H-R solution has only been verified in particular systems, the appearance of full circulating falling photographs of water droplets in castor oil does suggest the possibility of internal circulation [24,55,56]. However, researchers such as Bond and Newton have argued that the terminal velocity of small oil droplets should be between the Stokes solution and the H-R solution [57,58]. For this reason, Boussinesq put forward the view that the oil droplets will lack internal circulation because an interfacial monolayer which acts as a viscous membrane, constructed a constitutive equation including surface shear viscosity, surface dilational viscosity and surface tension, and finally obtained an exact solution in the case of creeping flow analogous to the H-R result (Equation (3)). But it is very difficult to obtain reliable measurements for parameter C [59,60]. In response to the discrepancy between the theoretical analysis and the experimental results, many researchers tend to believe that it is due to the impurity of the system, the H-R solution is the terminal velocity of the oil droplets in the pure system, and the Stokes formula is the terminal velocity when the surface of the oil droplets is completely covered with contaminants [61]. This view was reached subsequently, and some researchers proposed a stagnant-cap model to analyze the motion mechanism of oil droplets in this case [56]. Devies and Rideal also defined a new parameter called “a degree of circulation” Z based on the theory of surface-active pollutants, and used this parameter to revise the Stokes solution in order to provide a more reasonable explanation of the experimental results [50]. Ervik and Bjørklund later reinterpreted this problem using a continuous-interface model that takes into account the normal and tangential interfacial stresses, and solved the Stokes equation for oil droplets under different interfacial tensions, thus obtaining the critical radius between the two velocity extremes [55]. Relevant theoretical studies on the low Re range are still continuing, but it is difficult to have a completely uncontaminated pure system in the real environment [62], which leads to the fact that the surface of small oil droplets will inevitably be covered with different degrees of active contaminants; the difference between small oil droplets and rigid spheres is negligible [63]. Therefore, while constructing relevant mathematical or oil spill models, many researchers will directly choose the simpler Stokes sink velocity equation as the formula for the terminal velocity of oil droplets in the low Re range within the error tolerance [64,65].

As the particle size of oil droplets increases, Re gradually increases, the influence of oil droplets’ shape, wake flow, etc., gradually appears, and it is impossible to obtain the analytical solution via theoretical analysis, which requires an empirical or semi-empirical formula to calculate the terminal velocity of oil droplets. Thorsen et al. proposed an empirical formula for high interfacial tension systems that can be applied to calculate the terminal velocity of oil droplets in the range of 400 < Re < 900 [38]. Based on Hu and Licht’s study, Klee and Treybal combined their own experimental results to establish a new empirical equation for Re. After simplification, a set of two-stage empirical formulae was proposed for calculating the terminal velocity of oil droplets that can be applied to a wider range of Re [40,49,66]. But it is generally difficult to obtain the agreement of the majority of researchers due to the large contingency in the establishment process of purely empirical formulas. Subsequently, the U.S. Defense Technical Information Center, for its part, proposed a set of two-stage semi-empirical formulas in a report published in 1982, i.e., applying the Stokes solution at low Reynolds numbers and the Reynolds solution at high Reynolds numbers (Equation (7)), obtained when the drag coefficient is approximately constant at 0.5 [51]. This equation was later fully used by Elliot and many other researchers in constructing mathematical models of oil spills in marine environments [67,68,69]. The emergence of the multistage formula is an inevitable result of the development of research, because the increase in the particle size of oil droplets will inevitably drive the instability of the surrounding fluid flow and changes in the law of motion before and after. However, the motion of oil droplets will gradually go through multiple stages such as the emergence of front-to-back asymmetry of the steady flow, flow separation, wake stabilization, wake instability, etc. [24]. Thus the two-stage equation methodology will inevitably lead to the following problems: either the Re area that can be covered is small or the error in the transition zone between the two equations is severe. To this end, Clift and Grace et al. integrated the research results of several researchers and synthesized the findings that are consistent with liquid–liquid two-phase flow, including bubbles and raindrops, and proposed a more detailed multistage formulation using the shape of fluid particles as the classification criterion [24,70]. This formula starts from the most basic definition of Re, establishes the correlation equation between several dimensionless parameters and Re, and then derives the formula for the terminal velocity of fluid particles. It is currently the most widely used formula due to the large quantity of reference data and the wide range of Re included [71]. Later on, a further integration and analysis of the study by Clift was carried out by Zheng and Yapa which has determined the terminal velocity calculation equations and the corresponding critical diameters for three shapes of fluid particles [17]. Considering the relevance of the system purity, Grace, Clift and Wegener et al. introduced a correction function $\mathrm{“}\Gamma ”$ to the original formulation, thus obtaining a formulation that can be directly applied to pure systems [24,31,53].

The three-stage formulation proposed by Clift and Zheng et al. can basically meet the needs of fluid particle settling/floating problems that do not require high accuracy currently, but this method still has some problems that need to be considered. On the one hand, despite the refinement of the division of critical diameter by Zheng and Yapa, the division about the boundary between ellipsoid and spheroid is still not rigorous, and the two formulas cross in the scope of application. On the other hand, although bubbles and oil droplets belong to the same fluid particles’ category, the properties of $\mathsf{\gamma }\to 0$ and $\mathsf{\kappa }\to 0$ make the bubble motion belong to the extreme case of fluid particle motion, which was described by Licht et al. earlier [40]. Therefore, the use of a large number of bubble experimental results makes the accuracy of the formula doubtful when applied to oil droplets, and it may be necessary to correct or modify each paragraph of the formula and its applicability range by more oil droplet experiments. The same question applies to the formula proposed by Kelbaliyev and Ceylan [54].

2.3. The Determination of Drag Coefficient C_D in Oil Droplets Motion

Drag force is one of the most important forces controlling the motion of oil droplets in water, which is usually defined by the dimensionless drag coefficient C_D [46]. Hence, the determination of C_D is essential to obtain accurate drag force and terminal velocity [48].

A larger C_D results in greater drag force and lower terminal velocity [72]. Similar to the oil droplets’ terminal velocity, the determination of the oil droplets C_D by theoretical analysis is only applicable to the low Re stage, where C_D is estimated in most cases from experimental data and empirical correlations, especially for droplets with large Re values with corresponding flow characteristics (boundary layer separation, hydrodynamic wake, drag crisis, etc.) [37,54,73]. It is generally accepted that the empirical equation for the oil droplets’ drag coefficient is progressively dependent on the range of Re [54]. The variation process of oil droplets C_D with Re can be usually divided into the three stages. Low Re stage C_D decreases gradually with increasing Re. Medium Re stage C_D does not vary with Re and is basically stable at a fixed value (this stage sometimes does not exist). High Re stage C_D increases rapidly in a very small Re range until the oil droplet reaches a critical particle size and then breaks into small droplets [31,40,48,74]. Some researchers, however, considering the existence of deformation and the internal circulation of oil droplets, suggest that the empirical or semi-empirical equations for C_D can be established by deriving a modified model in terms of We, Mo, Eo and Oh numbers and based on these [41,48,54].

Table 2. The formula for calculating the C_D of oil droplets.

Time, Proposer	Expression		Scope of Application/Conditions		Eq.
1850, Stokes	${\mathrm{C}}_{\mathrm{D}}=\frac{24}{\mathrm{R}\mathrm{e}}$		$\mathrm{R}\mathrm{e}\ll 1$		(12)
1911, Hadamard and Rybczynski; 1962, Levich [61]	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}(\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }})$		$\mathrm{R}\mathrm{e}\ll 1$		(13)
1935, Schiller and Naumann; 2018, Cui et al. [75]	${\mathrm{C}}_{\mathrm{D}}=\frac{24(1+0.15{\mathrm{R}\mathrm{e}}^{0.687})}{\mathrm{R}\mathrm{e}}$		$\mathrm{R}\mathrm{e}\le 1000$		(14a)
	${\mathrm{C}}_{\mathrm{D}}=0.44$		$\mathrm{R}\mathrm{e}>1000$		(14b)
1964, Taylor and Acrivos [43]	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}(\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }})\left[1+\frac{\mathrm{R}\mathrm{e}}{16}\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}+\frac{1}{40}{\left(\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}\right)}^2{(\frac{\mathrm{R}\mathrm{e}}{2})}^2\mathrm{l}\mathrm{n}\frac{\mathrm{R}\mathrm{e}}{2}\right]$		$\mathrm{R}\mathrm{e}<1$		(15)
1976, Rivkind and Ryskin [76]	${\mathrm{C}}_{\mathrm{D}}=\frac{1}{1+\mathsf{\kappa }}\left[\mathsf{\kappa }(\frac{24}{\mathrm{R}\mathrm{e}}+\frac{4}{{\mathrm{R}\mathrm{e}}^{1/3}})+\frac{14.9}{{\mathrm{R}\mathrm{e}}^{0.78}}\right]$		$10\le \mathrm{R}\mathrm{e}\le 100$		(16)
1978, Clift et al. [24]; 1963, Hamielec et al. [77]	${\mathrm{C}}_{\mathrm{D}}=\frac{3.05(783{\mathsf{\kappa }}^2+2142\mathsf{\kappa }+1080)}{(60+29\mathsf{\kappa })(4+3\mathsf{\kappa })}{\mathrm{R}\mathrm{e}}^{-0.74}$		$4<\mathrm{R}\mathrm{e}<100$		(17a)
	${\mathrm{C}}_{\mathrm{D}}=\frac{48}{\mathrm{R}\mathrm{e}}\left[1+\frac{3\mathsf{\kappa }}{2}+\frac{{\left(2+3\mathsf{\kappa }\right)}^2}{{\mathrm{R}\mathrm{e}}^{1/2}}({\mathrm{B}}_1+{\mathrm{B}}_2\mathrm{l}\mathrm{n}\mathrm{R}\mathrm{e})\right]$ (${\mathrm{B}}_1$ and ${\mathrm{B}}_2$ are functions of $\mathsf{\kappa }\mathsf{\gamma }$; see 1978, Clift et al. for details)		$\mathsf{\kappa }<2,\mathrm{R}\mathrm{e}>50$		(17b)
	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{3}$		$\mathrm{R}\mathrm{e}>150,\mathrm{E}\mathrm{o}\ge 40$		(17c)
	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}\left(\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}\right)+\frac{8}{3}$		$\mathrm{M}\mathrm{o}\ge 150,\mathrm{a}\mathrm{l}\mathrm{l} \mathrm{R}\mathrm{e}$		(17d)
1998, Inaba and Sato [78]	${\mathrm{C}}_{\mathrm{D}}=2.28{\mathrm{R}\mathrm{e}}^{-0.193}{\mathrm{W}\mathrm{e}}^{-0.0326}$		$1<\mathrm{R}\mathrm{e}<80$		(18a)
	${\mathrm{C}}_{\mathrm{D}}=0.195{\mathrm{R}\mathrm{e}}^{0.344}{\mathrm{W}\mathrm{e}}^{0.622}$		$\mathrm{R}\mathrm{e}>80$		(18b)
2001, Ceylan et al. [74]	$\begin{array}{c}\hfill {\mathrm{C}}_{\mathrm{D}}=2.4{\mathrm{K}}_{13}\left(\mathrm{R}\mathrm{e},\mathrm{n}\right)+{\mathrm{K}}_{23}\left(\mathrm{R}\mathrm{e},\mathrm{n}\right)\hfill \end{array}$	of which: $$\begin{gathered} \begin{array}{c}\hfill {\mathrm{K}}_{13}\left(\mathrm{R}\mathrm{e},\mathrm{n}\right)=1-{\displaystyle \sum _{\mathrm{n}=0}^{\mathrm{\infty }}}{\mathrm{A}}_{\mathrm{n}}(\mathsf{\kappa }){\mathrm{R}\mathrm{e}}^{2\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}(-{\mathrm{m}}_{\mathrm{n}}(\mathsf{\kappa }){\mathrm{R}\mathrm{e}}^{\mathrm{n}+1})\hfill \\ \\ \hfill {\mathrm{K}}_{23}\left(\mathrm{R}\mathrm{e},\mathrm{n}\right)=1-{\displaystyle \sum _{\mathrm{n}=1}^{\mathrm{\infty }}}{(-1)}^{\mathrm{n}+1}{\mathrm{B}}_{\mathrm{n}}{\mathrm{R}\mathrm{e}}^{\frac{-2\mathrm{n}}{3}}\mathrm{e}\mathrm{x}\mathrm{p}({\mathrm{K}}_{\mathrm{n}}{\mathrm{R}\mathrm{e}}^{\frac{-\mathrm{n}}{4}})\hfill \end{array} \end{gathered}$$	$0.1<\mathrm{R}\mathrm{e}<{10}^4$		(19)
2001, Feng and Michaelides [79]	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}\left(\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}\right)\left[1+0.05\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}\mathrm{R}\mathrm{e}\right]-0.01\frac{2+3\mathsf{\kappa }}{1+\mathsf{\kappa }}\mathrm{R}\mathrm{e}\mathrm{l}\mathrm{n}\mathrm{R}\mathrm{e}$		$0\le \mathrm{R}\mathrm{e}\le 20$		(20)
2002, Saboni and Alexandrova [80]	${\mathrm{C}}_{\mathrm{D}}=\frac{\left[\mathsf{\kappa }(\frac{24}{\mathrm{R}\mathrm{e}}+\frac{4}{{\mathrm{R}\mathrm{e}}^{1/3}})+\frac{14.9}{{\mathrm{R}\mathrm{e}}^{0.78}}\right]{\mathrm{R}\mathrm{e}}^2+40\frac{3\mathsf{\kappa }+2}{\mathrm{R}\mathrm{e}}+15\mathsf{\kappa }+10}{(1+\mathsf{\kappa })(5+{\mathrm{R}\mathrm{e}}^2)}$		$\mathrm{R}\mathrm{e}\le 400, 0.01\le \mathsf{\kappa }\le 1$		(21)
2002, Rodi and Fueyo [81]; 2011, Kelbaliyev [23]	${\mathrm{C}}_{\mathrm{D}}=\frac{16}{\mathrm{R}\mathrm{e}}$		$\mathrm{R}\mathrm{e}<1.5$		(22a)
	${\mathrm{C}}_{\mathrm{D}}=\frac{16}{{\mathrm{R}\mathrm{e}}^{0.78}}$		$1.5<\mathrm{R}\mathrm{e}\le 80$		(22b)
	${\mathrm{C}}_{\mathrm{D}}=\frac{49.9}{\mathrm{R}\mathrm{e}}\left(1-\frac{2.21}{{\mathrm{R}\mathrm{e}}^{1/2}}\right)+1.17\times {10}^{-8}{\mathrm{R}\mathrm{e}}^{2.615}$		$80<\mathrm{R}\mathrm{e}\le 1530$		(22c)
	${\mathrm{C}}_{\mathrm{D}}=2.61$		$\mathrm{R}\mathrm{e}>1530$		(22d)
2006, Myint et al. [30]	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}(\frac{2+3\mathsf{\kappa }+3\mathrm{C}/\mathsf{\mu }}{1+\mathsf{\kappa }+\mathrm{C}/\mathsf{\mu }})(1+0.15{\mathrm{R}\mathrm{e}}^{0.687})$		$-11.6<{\log }_{10}\mathrm{M}\mathrm{o}<$−0.9, $1.7\times {10}^{-1}<\mathrm{R}\mathrm{e}<2\times {10}^2$, $1.7\times {10}^{-2}<\mathrm{E}\mathrm{o}<12.1$, $0.1<\mathsf{\kappa }<100$	clean systems	(23a)
	${\mathrm{C}}_{\mathrm{D}}=\frac{24}{\mathrm{R}\mathrm{e}}(1+0.15{\mathrm{R}\mathrm{e}}^{0.687})$			fully contaminated systems	(23b)
2007, Kelbaliyev and Ceylan [54]	${\mathrm{C}}_{\mathrm{D}}=\frac{8}{\mathrm{R}\mathrm{e}}\left[1+\frac{1}{1-0.5{(1+250{\mathrm{R}\mathrm{e}}^5)}^{-2}}\right]$		$0.1\le \mathrm{R}\mathrm{e}<0.5$		(24a)
	${\mathrm{C}}_{\mathrm{D}}=\frac{16}{\mathrm{R}\mathrm{e}}{\left[1+{(\frac{\mathrm{R}\mathrm{e}}{1.385})}^{12}\right]}^{1/55}+\frac{8}{3}\frac{{\mathrm{R}\mathrm{e}}^{4/3}{\mathrm{M}\mathrm{o}}^{1/3}}{24\left(1+{\mathrm{M}\mathrm{o}}^{1/3}\right)+{\mathrm{R}\mathrm{e}}^{4/3}{\mathrm{M}\mathrm{o}}^{1/3}}$		$$\begin{gathered} 0.5\le \mathrm{R}\mathrm{e}<100 \\ 9\times {10}^{-7}\le \mathrm{M}\mathrm{o}\le 7 \end{gathered}$$		(24b)
2010, Feng [42]	${\mathrm{C}}_{\mathrm{D}}=\frac{24{\mathrm{a}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{R}\mathrm{e}}\left[1+0.1935{({\mathrm{a}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\mathrm{R}\mathrm{e})}^{0.6305}\right]$		$20\le \mathrm{R}\mathrm{e}\le 260$		(25)
2017, Shao et al. [41]	the unsteady parameter A: $\mathrm{A}=(\frac{{\mathsf{\rho }}_{\mathrm{d}}}{{\mathsf{\rho }}_{\mathrm{c}}}-1)\frac{\mathsf{\sigma }\mathrm{W}\mathrm{e}}{{\mathsf{\rho }}_{\mathrm{c}}{\mathrm{U}}^4}\frac{\mathrm{d}\mathrm{U}}{\mathrm{d}\mathrm{t}}$	${\mathrm{C}}_{\mathrm{D}}={\mathrm{C}}_{\mathrm{D}\mathrm{S}}-1.3\mathrm{A}-0.4$ where:                               ${\mathrm{C}}_{\mathrm{D}\mathrm{S}}=\frac{24}{\mathrm{R}\mathrm{e}}{\left[1+18.5{\mathrm{R}\mathrm{e}}^5+{(\frac{\mathrm{R}\mathrm{e}}{2})}^{11}\right]}^{1/30}+\frac{4}{9}\frac{{\mathrm{R}\mathrm{e}}^{4/5}}{{330+\mathrm{R}\mathrm{e}}^{4/5}}$	$0<\mathrm{A}\le 2$		(26a)
		${\mathrm{C}}_{\mathrm{D}}={\mathrm{C}}_{\mathrm{D}\mathrm{S}}-0.7\mathrm{A}-1.0$ where:                               ${\log }_{10}{\mathrm{C}}_{\mathrm{D}\mathrm{S}}=-0.696 +1.259{\log }_{10}\mathrm{R}\mathrm{e}-0.465{\left({\log }_{10}\mathrm{R}\mathrm{e}\right)}^2 +0.045{\left({\log }_{10}\mathrm{R}\mathrm{e}\right)}^3$	$2<\mathrm{A}\le 6$		(26b)
2019, Joshi et al. [63]	${\mathrm{C}}_{\mathrm{D}\mathrm{\infty }}=14{\mathrm{R}\mathrm{e}}^{-0.53}$		$1<\mathrm{R}\mathrm{e}<14.5{\mathrm{M}\mathrm{o}}^{-0.154}$		(27a)
	${\mathrm{C}}_{\mathrm{D}\mathrm{\infty }}=0.115{\mathrm{R}\mathrm{e}}^{1.25}{\mathrm{M}\mathrm{o}}^{0.27}$		$14.5{\mathrm{M}\mathrm{o}}^{-0.154}\le \mathrm{R}\mathrm{e}\le 17.95{\mathrm{M}\mathrm{o}}^{-0.1985}$		(27b)
	${\mathrm{C}}_{\mathrm{D}\mathrm{\infty }}=2.66$		$\mathrm{R}\mathrm{e}>14.5{\mathrm{M}\mathrm{o}}^{-0.154}$		(27c)

summarizes some of the calculations of the C_D of oil droplets under different limiting conditions.

As can be seen from the table, the current research on the C_D relationship formula is divided into three areas. (1) Improving the applicability or accuracy of the formulas by making complex corrections to the underlying formulas. (2) Adopting the approach of relatively simple multi-segment formulas to expand the applicable range. (3) Between the first two areas, the proposed multi-segment formulas contain both complex modified equations and simple equations. But no matter which way of thinking, the research is based on the Stokes solution (Equation (12)) and H-R solution (Equation (13)). Taylor, Rivkind, Feng ([79]), Saboni, Feng ([42]) and Ceylan are typical representatives or researchers who have utilized the first approach to derive the formula. Taylor and Acrivos obtained a more accurate relation for oil droplet C_D at low Re (Re < 1) by the singular-perturbation solution of the axisymmetric equation of motion, which includes the first-order effect of deformation on the drag force [43]. Feng and Michaelides instead chose to use the finite difference method to extend the applicability of the formula to 0 ≤ Re ≤ 20 by adding a natural next order expansion in terms of Re [79]. As for the case of medium Re range, Rivkind and Ryskin obtained the C_D relation equation applicable to 10 ≤ Re ≤ 100 by numerically integrating the complete N-S equation [76]. A further extension of the equation proposed by Rivkind and Ryskin was developed by Saboni and Alexandrova to increase its applicability to Re = 400, but the equation is currently only for low viscosity oil droplets [80]. However, considering the case that the oil droplet will deform when moving in a viscous fluid and cannot maintain a fixed shape like a rigid particle, Feng used the C_D relation equation for rigid particles proposed by Clift et al. as the basis and used the radius of the cross-sectional area of the oil droplet, a_max, to describe the drag force, thus turning the C_D equation for deformable oil droplets into Equation (25) [24,42]. Ceylan et al. used droplet deformation resistance and frictional resistance as the two terms of the equation, respectively, and established the empirical relationship of oil droplet C_D by the approximate series method, which greatly broadened the scope of application of Re, but there were greater difficulties in the practical calculation of this formula [74]. Nevertheless, even for the most complex single formula, it is still inherently deficient in that it cannot be applied to the full Re range. Consequently, multi-stage formulas are also continuously improved in order to apply to larger oil droplet particle sizes and to satisfy a larger Re range. Although Schiller and Naumann corrected the Stokes solution to improve its applicability, no distinction was made between rigid particles and oil droplets, i.e., the oil droplets are invariably spherical particles by default, so the error is large and only applicable to simple calculations when accuracy is not required [75]. Kelbaliyev’s incorporation of Mo numbers to describe the degree of droplet deformation greatly improves the accuracy of the results but also limits the range of use of the formula [54]. While Clift and Zheng et al. proposed a complete set of terminal velocity calculation equations, there is no unified set of C_D relations, and each relation is proposed and applied under different applicable conditions, which are not coherent and lack consistency [17,24]. Inaba and Sato followed the idea of Klee’s earlier study and used the We number to describe the effect of interfacial tension on C_D, an important factor causing the deformation of droplets. The form of the equation is reasonable, but the parameters in it are empirical values and vary greatly between different droplet motions [49,78]. So, this equation cannot be directly applied to oil droplets’ motion, and still needs to be corrected for the parameters by experiment. Rodi and Fueyo unified the division boundaries and somewhat simplified the equations based on Clift’s study, but the conclusions, such as Equation (22b), were drawn for bubbles, which also have limitations [81]. As study in this area progressed, more details of the droplets’ motion were gradually taken into consideration. Myint et al. and Joshi et al. distinguished the difference between the droplet C_D in a completely clean system and a contaminated system in their study [30,63]. Shao et al. considered that there is no real steady motion during the droplet motion, i.e., the shape of the droplet is always in a state of change, and proposed the idea of unsteady drag coefficient and added unsteady parameters to correct the relationship equation of C_D [41]. In summary, it is necessary to establish a set of relations for the C_D of oil droplets in the full Re range with a unified classification standard for the practical application environment.

3. Migration and Transformation Characteristics of Oil in Inland Riverine Environment

The typical characteristics of an inland riverine environment are relatively narrow channel width, relatively shallow depth (the average depth of the Yangtze River is only 30–40 m) and the tendency for the water body to contain particles such as sediment, which is generally described by using the average flow velocity to characterize its movement [82,83]. Although the global surface area of rivers accounts for only 0.47% of the Earth’s land surface and inland waterway oil spills are generally small, oil spill accidents in inland rivers are more frequently occurring and often near densely populated areas, which can easily pollute the drinking water sources [25,84]. Therefore, the study and characterization of oil migration and transformation in inland waters should not be neglected. When the spilled oil enters the inland water environment, due to the turbulent movement of the water body, the spilled oil will dispersed into oil droplets and submerge in the water body [85]. Due to the comprehensive effect of the surrounding current structure, volume fraction of oil droplets and the physico-chemical properties of the oil droplets, the oil droplets will first selectively undergo zero or more behaviors such as collision, aggregation or fragmentation. After the above behaviors are basically stabilized, the large oil droplets will rise to the water surface and become an oil film, while the small oil droplets will suspend in the water column or continue to diffuse downward [47,67]. Eventually, the spilled oil will achieve dynamic stabilization of the concentration within a certain spatial distribution range.

3.1. Behavior of Oil Droplets Collision, Coalescence or Breakage

The turbulence of water flow or the increase in the volume fraction of oil droplets in the water column will lead to collision, coalescence or breakage between oil droplets, thus changing the original size and movement state of oil droplets [86]. The interrelationship between these behaviors is given in Figure 3. At present, the main causes of oil droplets collision can be summarized as four points: eddy current caused by continuous phase turbulent flow, velocity gradient or shear effect of environmental fluid, velocity difference caused by buoyancy or gravity, and wake entrainment [27,71]. Wu et al. believed that the turbulent action of water flow would induce the random collision of oil droplets, but it was limited to adjacent oil droplets [87] because the long-distance interaction was driven by large eddy currents, which transported oil droplet groups without causing significant relative movement [88]. However, Gong et al. believed that the movement of droplets can be jointly driven by small vortices and large vortices, and it is necessary to consider the contribution of vortices of different sizes to the coalescence of oil droplets of a given size [89,90]. Although they have different opinions on the effects of large eddy currents on oil droplet collision, they both indicate that the distance between oil droplets is an important factor for collision in continuous phase turbulence. Beyond that, the particle size of oil droplets has also been considered by many researchers. They have suggested that the difference in particle size causes a difference in velocity between the oil droplets, which results in mutual collisions [91]. Nevertheless, Chesters expressed that when the particle size of oil drops is far smaller than the eddy current size carrying and dissipating fluid energy, the viscous force plays a leading role in driving the movement and collision of oil droplets [92].

The collision of oil droplets is the premise of coalescence. When two oil droplets have been in contact for a long enough time during the collision process and the liquid film is drained and broken, it indicates that the oil droplets have coalesced, otherwise bouncing occurs [93,94]. Meanwhile, one of the important factors affecting the coalescence or rebound behavior of oil droplets is considered to be the Reynolds number [94]. The oil droplet coalescence frequency is determined by the collision frequency and the coalescence efficiency [27], which can generally be expressed by Equation (28):

$$\mathsf{\Gamma }\left({\mathrm{d}}_1,{\mathrm{d}}_2\right)=\mathrm{h}\left({\mathrm{d}}_1,{\mathrm{d}}_2\right)\mathsf{\lambda }({\mathrm{d}}_1,{\mathrm{d}}_2)$$

(28)

where $\mathrm{h}\left({\mathrm{d}}_1,{\mathrm{d}}_2\right)$ is the oil droplet collision frequency; $\mathsf{\lambda }\left({\mathrm{d}}_1,{\mathrm{d}}_2\right)$ is the oil droplet coalescence efficiency.

There are different ways to calculate the frequency of oil droplet collision caused by different mechanisms; the calculation of oil droplet coalescence efficiency also varies depending on the oil droplet coalescence regime. Also, both the droplet collision frequency and the coalescence efficiency rely heavily on arbitrary assumptions and empirical relationships [27,95]. Therefore, the calculation of oil droplet coalescence frequency always tends to be empirical. Oil droplets break up when the energy gained from the turbulent continuous-phase vortex is sufficient to compensate for the increased surface energy due to the expansion of the surface area of the droplet [96]. This behavior can be generally divided into two phases: the squeezing stage controlled by the squeezing pressure and the pinch-off stage controlled by viscous stresses of both phases and surface tension [97]. The construction of the oil droplet rupture model is currently based on a modified version of the Kolmogorov and Narsimhan model [98]. The collision, coalescence and breakage behaviors of oil droplets are complex behaviors, which are seldom incorporated in the construction of oil spill models at present; the particle size distribution of oil droplets in steady state is generally chosen to be applied directly.

3.2. Vertical Diffusion Characteristics of Oil

The diffusion process of oil under water is completed by the combined effect of buoyancy and turbulence [19,99]. The vertical concentration distribution of the oil is the most intuitive expression of the vertical diffusion, and the vertical turbulent diffusion coefficient is generally chosen to describe the degree of diffusion [100]. The diffusion coefficient, which first appeared in Fick’s diffusion law, is a fundamental physical quantity in the diffusion theory model and a basic parameter in transport forecasting, and has been a fundamental part of the study of diffusion laws for a long time.

The current research on the diffusion of oil to the lower water column in inland riverine environment is mainly based on sediment diffusion theory, from both the characteristics of vertical concentration distribution and vertical turbulent diffusion coefficient. However, most of them are still at the stage of qualitative analysis, i.e., the shear flow velocity magnitude is positively correlated with the vertical diffusion of oil droplets [47], lacking systematic theoretical analysis and complete experimental data. According to the oil droplet diffusion mechanism, the differential equations of oil droplet turbulent diffusion in one- and two-dimensional cases can be established, resulting in the distribution equations of oil concentration along the water depth (specifically, as follows, Equations (29) and (30)), whose results show that the oil concentration decreases exponentially along the water depth. But Imanian et al. stated that the oil concentration distribution is close to parabolic distribution in both time and space [101]. This is mainly due to the fact that the above theoretical derivation equations are realized given the corresponding initial boundary conditions and a large number of assumptions; thus, there are some difficulties and errors between them and the experimental results and practical applications:

$$\mathrm{one}-\mathrm{dimensional}: \mathrm{C}\left(\mathrm{z},\mathrm{d}\right)={\mathrm{C}}_0\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\mathrm{U}}_{\mathrm{T}}}{{\mathrm{D}}_{\mathrm{z}}}(\mathrm{z}-\mathrm{a})\right]$$

(29)

$$\mathrm{two}-\mathrm{dimensional}: \mathrm{C}\left(\mathrm{x},\mathrm{z},\mathrm{d}\right)={\mathrm{C}}_0\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\mathrm{U}}_{\mathrm{T}}+\mathrm{u}-{\mathrm{U}}_{\mathrm{w}}}{{\mathrm{D}}_{\mathrm{x}}+{\mathrm{D}}_{\mathrm{z}}}(\mathrm{z}-\mathrm{x}-\mathrm{a})\right]$$

(30)

where $\mathrm{C}\left(\mathrm{z},\mathrm{d}\right)$ and $\mathrm{C}\left(\mathrm{x},\mathrm{z},\mathrm{d}\right)$ are the oil concentration; d is the equivalent diameter of oil droplets; ${\mathrm{C}}_0$ is the oil concentration at water depth a; ${\mathrm{U}}_{\mathrm{T}}$ is the terminal velocity of oil droplets; ${\mathrm{D}}_{\mathrm{z}}$ is the is the diffusion coefficient in the z direction; z is the depth of water at a given location; u is the current velocity; ${\mathrm{U}}_{\mathrm{w}}$ is the oil film drift velocity; ${\mathrm{D}}_{\mathrm{x}}$ is the diffusion coefficient in the x direction; and x is the horizontal distance.

Two main methods were used to study the vertical dispersion coefficients for oil. One is to measure the correlation moment of oil droplet concentration fluctuation and water velocity fluctuation in water flow, which is calculated according to the gradient assumption of Reynolds average process of instantaneous equation. The other is based on the Schmidt equation, and the measured oil droplet concentration distribution is used for back calculation. According to the above methods and the current situation of sediment-related research, the vertical diffusion coefficients ${\mathrm{D}}_{\mathrm{z}}$ of petroleum pollutants in two-dimensional open channel uniform flow have different vertical distribution types, such as parabola, parabola-constant, and gradually increasing towards the water surface [102]. Graf and Cellino used Method 1 to obtain the instantaneous values of flow velocity and concentration of sediment laden flow in open channel, and calculated the vertical diffusion coefficient, which showed a typical parabolic distribution type [103]. Although Jobson believes that the vertical diffusion coefficient should be considered along the water depth direction by dividing the inner layer and outer layer, the results of his empirical formula are similar to the parabolic distribution [104]. Coleman used relevant flume experimental data to reverse the distribution of vertical diffusion coefficient by Method 2, and finds that it is distributed in sections with the increase in water depth [102], which is consistent with the relevant experimental results of Anderson [105]. Nikora and Goring believe that the vertical diffusion coefficient increases monotonously towards the water surface [106]. Currently, the vast majority of theoretical and experimental results are derived from sediment. Although petroleum pollutants have some similarities with sediment, Elliott et al. believe that buoyancy will produce resistance to the downward movement of oil droplets, in the actual vertical mixing there is a depth limit, and can not appear at the depth of the extreme value of the results [67]. Thus, there is still a lot of work to be completed on the study of the vertical diffusion coefficient of oil.

4. Migration and Transformation Characteristics of Oil in Oceanic Environment

The oceanic environment is both similar and different from inland riverine environments in that it has a wider range of space and depth. Waves are widespread as the most typical feature of water movement in the oceanic environment, which specifically manifest as the disturbance of the free surface of incompressible fluids [107]. Three types of waves, regular waves, irregular waves and breaking waves, are generally involved in the research, which is usually described by wave height, wave length, wave energy, etc. [108,109,110]. The details can be seen in Figure 4. Once an oil spill occurs in the oceanic environment, it is usually a large or medium-sized oil spill, with a spill quantity of at least 700 tons [19]. Thus, when large quantities of oil are spilled on the ocean in a short period of time, the oil film which first forms is thick so that small ocean waves passing through it tend to be damped out; meanwhile, the film itself appears to ride on any large swells present, at which point the breakup of the film into oil droplets is inhibited [111]. As the oil film continues to spread and become thinner, the high energy and strong vortex generated by the waves’ action begins to accelerate the generation of oil droplets [111], increasing the distance between each oil droplet which enters the water column, and also inhibiting the coalescence of oil droplets into large oil droplets to resurface by reducing the frequency of oil droplet collision; the spread of oil droplets and the degree of diving increases until the dynamic stability of vertical mixing is achieved [108,112,113].

4.1. Oil Droplet Size Distribution

After continuous collision, coalescence and breakage of oil droplets in wave environment, the oil droplet population gradually reaches equilibrium and maintains a stable size distribution in the water column [114].

In the early stage, based on the famous Pierson–Moskowitz spectrum, Aravamudan supposed a simple linear relationship between the number distribution of the oil-into-water droplet group and the oil droplet size and provided a method to calculate the extreme value of the oil droplet size, eventually developing a marine oil spill model to predict the oil droplet size and its distribution in the water column [51,111,115]. But the assumptions are difficult to correspond to the actual situation due to the lack of relevant data support [116]. Subsequently, some researchers gradually began to summarize the oil droplet size distribution law through experimental research and establish the corresponding empirical formula [117,118,119]. Delvigne focused on quantifying the effects of oil viscosity and wave energy on the droplet size distribution. By analyzing the experimental results, they reached the conclusion that the size distribution of dispersed oil droplets is related to the stability of dispersion and that the droplet size distribution caused by a single breaking wave is not a steady-state distribution. Finally, they proposed the first algorithm to predict the droplet size distribution based on oil properties. However, due to the limitations of the experimental equipment, it only obtained the power–law relationship between the average droplet size d₅₀ and the maximum droplet size d_max and the above two factors, and the algorithm is only applicable to oils with a viscosity below 1000 cps [18,120,121].

Afterwards, Lefebvre proposed allowing the mathematical representation of measured oil droplet size distributions based on probabilistic or purely empirical considerations in the absence of any fundamental mechanism or model for establishing a theory of oil droplet size distributions. In general, the uses include normal, log-normal, Nukiyama–Tanasawa, Rosin–Rammler, and upper limit distributions [121,122]. Based on Delvigne’s study, Reed therefore extended the experiment to the non-Newtonian region with higher viscosity and showed that the log-normal function fitted the size distribution best for different weathering conditions [121], and this conclusion was also verified in the experiments of Cui [110]. With the extensive use of the laser particle size analyzer LISST-100X in related experiments, the experimental data have progressively become more detailed and precise [123,124]. Li and Lee et al. thereby conducted a series of experiments to more accurately describe the effects of regular waves, breaking waves, and the addition of dispersants on the particle size variation in oil droplets. Their results showed that the volume average diameter of dispersed oil droplets was always high under the action of regular waves, regardless of the presence of dispersant, ranging from 400 to 450 μm throughout the experimental period. However, under the action of breaking waves, the volume average diameter decreased rapidly from 450 μm to less than 200 μm within 10 min, even without dispersant, and remained low thereafter [124,125,126,127]. Li and Miller et al. not only further proved that the size distribution of oil droplets is related to wave intensity and physical properties of oil (viscosity, surface tension), but also established a model of semi-empirical volume mean size under breaking wave conditions based on dimensionless parameters; see also We and Oh [128].

The above studies are based on surface oil spills, but the oil droplet size distribution formed during subsea spills or deepwater oil and gas well blowouts need to be considered separately due to the differences in the causes of oil droplet generation and hydrodynamic conditions [129]. Tang and Masutani were the earliest to start related small-scale experiments, eventually identifying five instability patterns and showing that the progression of instabilities leads to oil droplet generation initially moving from the jet core filament to the jet surface [130]. Johansen et al. conducted deep in situ spill experiments in the Norwegian Sea at the Helland Hansensite (65°000 N, 04°500 E) in 844 m of water roughly 125 km off the coast of central Norway in order to more accurately understand the behavior of oil during deepwater spills [131]. A major problem with in situ experiments, however, is that the environmental and experimental conditions are very difficult to precisely control, and there is a large variability in the results [132]. To this end, Brandvik et al. designed a cylindrical tower pond with a diameter of 3 m and a height of 6 m on their own to fill the gap between the two by experimental results under a mesoscale setup. Their experimental results showed that the oil droplet particle size was uniformly distributed in the range of 2.5 to 500 μm in a logarithmic scale, and the relationship between the oil droplet particle size and nozzle diameter could be described by the equivalence equation of the Weber number scaling method [133]. Johansen et al. then proposed a new method for predicting droplet size distribution in the case of subsea oil and gas release based on this experimental result, and made a distinction regarding the addition of dispersant or not [134].

In addition to the empirical equations for the particle size distribution of oil droplets under specific aqueous environment conditions obtained directly through relevant experiments, other researchers have now chosen to use population dynamics models to investigate the oil droplet dynamic. This model starts from the initial oil droplet size, and uses the mass and momentum conservation equations to simulate the process of oil droplet breakage and coalescence. Bandara and Yapa used population dynamics equations to simulate the fragmentation and aggregation of oil droplets in turbulence for the first time, and fitted the experimental data well [5,135]. Zhao et al. further considered the influence of oil surface tension and viscosity on the results [118,136]. Subsequently, Nissanka and Yapa improved the algorithm on the above model framework, so that it can be used to calculate the corresponding oil droplet size distribution under different possible release conditions of underwater oil and gas blowout (μm to mm) [119]. At the same time, the above models are all based on the empirical relationship obtained from the limited available data to set the oil droplet fragmentation and polymerization coefficients. There are no independent test data to verify the coefficients, and the initial oil droplet size needs to be provided in the calculation process, which has certain limitations.

4.2. Vertical Mixing Characteristics of Oil

The current general use of vertical mixing is to describe the characteristics of the vertical motion of oil in the wave environment [19,137]. This process involves two main components, the total amount of oil entering the water column and the distribution of oil staying in the water column. Among them, since the oil dispersion volume/rate (sometimes called wave entrainment volume/rate) is directly related to the surface oil slick volume and diffusion movement, etc., most of the current studies focus on quantifying the relationship between wave energy or oil spill characteristics and oil dispersion volume/rate [18,108,113], while relatively little research has been conducted on the spatial and temporal distribution of the concentration of oil in the water column.

Delvigne and Sweeney earlier proposed an oil entrainment equation for the broken wave case which involves parameters such as wave energy, particle size distribution, and surface oil film coverage area [18]. Although this equation is empirical and needs to rely on proportionality constants to correct for different oils, Delvigne’s experiments were by far considered the most extensive and the most complete in terms of data and empirical models obtained. Subsequently, Tkalich and Chen used a kinetic model considering only the vertical exchange of oil droplets between the oil slick and the mixed layer, and a more precise parametric description of the wave breaking phenomenon, in order to apply the above model for more diverse environmental conditions [108,138]. To address the situation that previous related studies have mainly focused on oils with viscosities less than 1000 cP, Reed et al. extended the oil diffusion model to higher viscosities and even non-Newtonian regions by using the dimensional analysis method and considering both Weber number and dimensionless viscosity groups [121]. Most of the studies showed that the oil dispersion is limited by the oil viscosity and linearly related to the oil layer thickness, while the dispersion rate is independent of the oil properties, even the dispersant addition does not change the result substantially. However, the effect of wave strength, oil viscosity and dispersant dose on the oil dispersion rate is subject to a threshold value, and when the threshold value is exceeded, the oil dispersion will be affected [12,113]. The results of Zhao’s study, for example, show that the effect of waves on the oil spill entrainment rate is much greater than that of dispersants, but the coupling effect is only apparent when both exceed a threshold value and the effect can be cumulative over time [113].

Oil droplets entrained by waves into the water column will continue to vertically migrate due to turbulence or eddy effects, resulting in changes in temporal and spatial distribution. Li et.al first conducted a series of studies on this aspect and obtained the vertical oil concentration distribution in different sections, different time periods and different wave conditions. Their results show that the oil concentration increases with time with certain fluctuations, but eventually gradually levels off. The greater the wave energy, the greater the oil concentration in the water column, of which the most significant increase in oil concentration is near the surface. The penetration depth of dispersed oil in water is shallow, and the oil concentration at the bottom is very low, of which the most significant is in regular waves [126,139,140]. In Li’s study, wave conditions were mainly described by wave intensity, which could not reflect the relationship between wave characteristics such as wave period, wave height, etc., and experimental results. Parsa et al. therefore considered the effect of wave characteristics on the vertical concentration distribution of oil. Their study found that wave steepness has a direct effect on the oil dispersion in the water column, the oil concentration tends to vary linearly with wave steepness, and the oil concentration increases with the increase in wave height or the shortening of the wave period [109,141]. However, the existing research results obtained are basically qualitative conclusions, lacking theoretically rigorous quantitative results, and also cannot properly reflect the vertical diffusion characteristics of oil, while the characterizing parameters such as vertical diffusion coefficient are still mainly referred to the relevant forms of sediment in the construction of oil spill models.

4.3. Diffusion Characteristics of Oil under the Combined Action of Waves and Current

In coastal and estuarine areas, waves and water currents/tidal currents co-exist, the hydrodynamic structure becomes complex via the interaction between them and the turbulent kinetic energy generated by shear and wave breaking activities, and the movement pattern of petroleum pollutants is exceptionally complex and more stochastic [142,143,144]. Although 3D oil spill models containing complex conditions such as wind, currents, and waves have been developed, the vertical motion of oil droplets is generally determined by the Langeven equation, using stochastic perturbations from Markov’s chain for deceleration [143,145]. However, this method’s final result only reflects whether the oil particles invading into the water body will resurface to form an oil film, it cannot reflect the vertical concentration distribution of oil particles that are completely submerged into the water body or may continue to move to the lower water body. There are very few studies on the vertical dispersion characteristics of oil under the joint action of waves and currents, and only relevant simple qualitative conclusions are available [146], so the law of sediment movement proposed by Van Rijn can generally be directly referred to [99,147]. That is, the vertical diffusion coefficient of oil under the joint action of wave and current can be expressed as a linear superposition of the diffusion coefficient under the action of pure wave and pure current, and the diffusion coefficient needs to include both the vertical oil diffusion caused by wave current turbulence and the vertical oil diffusion caused by wave trajectory flow velocity [99,148]. The specific equation is as follows:

$${\mathsf{\epsilon }}_{\mathrm{w}\mathrm{c}}^{(1)}={({({\mathsf{\epsilon }}_{\mathrm{w}}^{(1)})}^2+{({\mathsf{\epsilon }}_{\mathrm{c}}^{(1)})}^2)}^{0.5}$$

(31)

where $({\mathsf{\epsilon }}_{\mathrm{w}}^{(1)})$ and $({\mathsf{\epsilon }}_{\mathrm{c}}^{(1)})$ are diffusion coefficients under the action of pure wave and pure current, respectively.

5. Conclusions and Perspectives

The impact on the surrounding ecological environment after an oil spill accident is self-evident. Therefore, the accurate prediction of the migration path and dispersion amount of oil spill in environmental waters is the critical key to rationalize the emergency response to accidents, which has received extensive attention. This review starts from the mechanism of the oil spill migration and transformation process, and firstly reviews the kinematic characteristics of a single oil droplet, which is the smallest moving unit of the oil spill, and focuses on several key aspects, such as the droplet shape, the trajectory, the terminal velocity, and the drag coefficient. Subsequently, considering the commonalities and differences between inland riverine and oceanic environments, different aspects of oil droplet collision, coalescence, breakage, size distribution, and vertical diffusion are discussed, respectively. Up until now, some reliable research conclusions have been obtained on the free-rising motion of single oil droplets, the particle size distribution of oil droplets in the oceanic environment and the vertical diffusion characteristics, but it is still necessary to continue to provide more accurate and applicable models or methods. There is still a big gap in the related research on the movement process of oil spills in inland rivers [149].

In conclusion, the current status of the existing research as well as future research directions and outlook are as follows:

• The related research on the free rising motion of a single oil droplet began earlier, and the research results are also more numerous. Early researchers, due to research conditions and other reasons for the limitations, have obtained results that are generally lower in precision. With the continuous refinement of related theories and equipment, more researchers choose to start from the subtle factors affecting the movement of oil droplets (e.g., the concentration and type of surface-active pollutants, different droplet shapes produced by different nozzle types, etc.), in order to obtain more accurate theoretical results and quantitative relationships. However, theoretical or experimental studies ultimately serve for engineering applications. Future research should focus more on the optimization of relevant parameters in engineering applications, mainly considering the following two points. (1) In capturing the range of relevant physical properties according to the possible types of oil spills, through theoretical analysis and experimental research, a more comprehensive set of models applicable to oil droplets can be integrated by adopting the method of re-proposing or amending the original correlation equation. (2) Surface-active pollutants are bound to exist in environmental waters but are difficult to quantify, and optimization results of relevant parameters within a certain error range or fluctuable range can be obtained based on the average situation under different water conditions.
• At present, the research on the migration and transformation characteristics of oil spills in the ocean environment is mainly divided into two categories: experimental research and numerical simulation. Although experimental research can precisely control the influence factors, the final results are more idealized. Numerical simulation reveals the actual situation more realistically, but its accuracy depends on the hydrodynamic conditions of the provided waters and the accuracy of the relevant parameters. Therefore, on the one hand, future experimental studies should be based on the existing research, through the gradual coupling of more hydrodynamic conditions (e.g., the combination of wave, current, etc.) in order to obtain results closer to the actual water environment. On the other hand, numerical simulation should be combined with experimental research, and the results of the two can be the prerequisite for and verify each other.
• The amount of accidental oil spills occurring in inland rivers is relatively small but more frequent, many studies on the migration and transformation process of oil spills in inland rivers are currently in a state of blankness, and a large number of relevant studies are urgently needed in the future. Initially, the research can be carried out with reference to the relevant research in the ocean environment. Firstly, the particle size distribution characteristics of oil droplets under different flow conditions are the basis for further research. Secondly, the vertical water depth of the river is much lower than that of the ocean, which makes it more important to determine the vertical dispersion characteristics of the oil spill compared to the marine environment. Finally, the interaction between oil spills and sediment, although not considered in this review, is also a focus for future research.