Numerical Simulation on Oil Spilling of Submarine Pipeline and Its Evolution
on Sea Surface

Due to the interaction and corrosion of the seawater, submarine pipelines are easy to be broken to spill oil. The special environment of subsea restricts the technical development of pipeline maintenance. Therefore, the study on the oil spilling model of submarine pipeline is very important for predicting the movement and diffusion of spilled oil, so that oil spilling traces and relating strategies can be determined. This paper aims to establish an oil spilling model of a submarine pipeline, study the movement characteristics of spilled oil in seawater by numerical simulation, and determine the traces, diffusion range, time to sea surface, etc. Then, the maximum horizontal migration distance (MHMD) with corresponding time are analyzed under different oil densities, spilling speeds and seawater velocities. Results show that the MHMD decreases first and then increases while the time to achieve the MHMD increases along with increasing oil density. The MHMD increases while the time to achieve the MHMD decreases, along with increasing spilling speed. Both the MHMD and corresponding time increase along with increasing seawater velocity. Based on numerical results, a correlation of spilling distance and spilling time is proposed to give fast and accurate predictions. After the oil reaches sea surface, oil expansion and transport are simulated. Euler-Lagrange method is used in the simulation. Dynamic and non-dynamic factors are considered. Results show that wind velocity and water velocity are dominant in dynamic factors. When they are large, spilled oil moves very fast with variable directions in complex flow field. Nondynamic factors such as evaporation, emulsion and solution mainly reduce the volume of oil film. They almost do not affect the direction and displacement of spilled oil. Quick response should be made for large wind and water velocities when the placement of oil boom is given. With the correlation and simulation, emergency responses can be guided effectively to reduce the impact of submarine oil pollution. The computational results benefit pollution control and environmental protection in marine petroleum engineering.


Introduction
With the development of submarine oil drilling and production [1], submarine pipelines become the main transportation manner of oil and gas in sea oil fields. During the operation of the submarine pipelines, some kinds of damages are inevitable, including internal corrosion by petroleum, external corrosion by seawater, stress corrosion by pressure and other forces. These damages lead to crack and leakage of the pipeline so that crude oil will spill, which causes large wastes of petroleum resources and severe pollutions of environment [2,3]. To minimize the economic losses and social risks induced by the oil spilling of submarine pipelines, fast and sufficient responses should be taken after the spilling. This relies on accurate prediction of oil diffusion in the entire spilling process.
At the initial stage of spilling, crude oil usually has a certain velocity. Moreover, the density of oil is smaller than seawater so that oil is also driven by buoyancy force. Thus, the initial stage of spilling can be considered as buoyant jet. List et al. [4] compared the experimental results with the dimensional analysis. They indicated that the entrainment coefficient of rising-up turbulent buoyant jet was varying continuously and had close relations with the Froud number and diffusivity of the buoyant jet. Fischer et al. [5] explained basic mechanism of buoyant jet underwater and solved the relation of the diffusion height with spilling speed and concentration. McDougall [6], Fannelop et al. [7] established oil spilling models in static water and simulated the vertical spilling diffusion, but the influence of surrounding water on spilled oil was not considered. Milgram [8] introduced the momentum coefficient into his plume model by considering the effect of turbulent fluctuation. Papanicolaou et al. [9] obtained the oil spilling trace in turbulent shear flow in experiments. Rye [10] used Fannelop et al.'s model [7] to simulate the oil spilling diffusion under different situations. Wright [11] and Wood [12] studied the asymptote solution for simple cases and established a simple buoyant jet flow model which can simulate limited cases.
Further developments of oil spilling models refer to a k-ε model [13] and the integral method, which includes Euler method and Lagrangian method. Euler method fixes the control volume to simulate the buoyant jet flow [14][15][16][17][18][19]. Lagrangian method tracks the traces of all fluid particles to simulate the flow field [20,21]. Yapa et al. [22] and Zheng et al. [23] established a three-dimensional (3D) oil spilling model for a submarine pipeline using the Lagrangian method and introduced effects of diffusion, solution and entrainment to simulate the buoyant jet flow induced by submarine oil spilling. Wang et al. [24] put the work of [22] and [23] forward to consider the emulsion of the spilled oil and improve the buoyant jet theory. Zheng et al. [25] further compared the simulation results with experimental data and found they agree well with each other. Based on this, Reed et al. [26] used a combination model to simulate the migration and volume of spilled oil. Dasanayaka et al. [27] applied a convection-diffusion model and a buoyant jet model to do the simulation and compared the results of the two models. Lardner et al. [28] improved Malačič's model [29] and tested parameters in plume model.
Gao et al. [30] established a prediction model of oil diffusion through an ostiole on submarine pipeline using a velocity-pressure coupling algorithm. They obtained the rising transport process of the spilled oil droplet by the buoyancy of seawater and established 2D and 3D models for predicting the oil spilling. Wang [31] and Xiao [32] established a micropore oil spilling model and simulated the formation, rising speed and scale of the oil droplet. The influence of the initial velocity of the spilled oil on the oil trace was not considered. Liao et al. [33] studied the migration process of spilled oil in deep sea of Norway. Simulation results agree well with field test data. Chen et al. [34] simulated the undersea oil spilling in static and wavy water environment but did not consider other factors. Liu [35] simulated the diffusion process of spilled oil in water using the volume of fluid (VOF) method and analyzed the effects of spilling speed, diameter of spilling hole, etc.
Spilled oil can expand on sea surface very quickly to pollute the environment. It is quite necessary to study oil spilling behaviors to reduce the environmental hazards. Numerical simulations can obtain the traces of spilled oil more easily than experiments and field tests. Prediction data from numerical simulations play more and more important roles in dealing with the oil spilling accidents. Blokker [36] first proposed an oil expansion model which assumes the oil film was circle and expanded everywhere. His model only considered gravity. Fay's expanded formula of oil spilling law was applied widely [37]. His model was divided by three stages with different equations. Sea surface was assumed to be static. Turbulence of seawater and viscosity of spilled oil were not considered. Properties of oil film were assumed to be constant in the process of oil expansion. Thus, this model has some limitations in applications. Based on Fay's theory, more researchers made improvements and developments. Mackay [38] first proposed the transitional relation between thick and thin oil films by improving Fay's model. Wang et al. [39] added the influence of the Coriolis force on oil spilling. Lehr et al. [40] found the asymmetric expansion phenomenon of oil film and introduced the effects of ocean currents and winds. Elliott et al. [41] modified Fay's formula via experiments. Diffusion coefficients in vertical and horizontal directions are determined to obtain results closer to reality. Venkatesh et al. [42] included variations of oil viscosity in oil expansion process and developed Fay's model. After 1980s, particle random diffusion theory was applied to the study of oil spilling models. Johansen [43] and Elliott [44] proposed an "oil particles" spilling model. They considered the continuous oil film as the combination of amount of oil particles and tracked the random movements of each particle. Zhao et al. [45] considered wind effects. Yang et al. [46,47] added solutions and emulsifications in the oil spilling model. Ji et al. [48] simulated the transport diffusion dominated by ocean currents, waves and winds. Zhang et al. [49] established a 3D model for oil spilling on sea surface using the concept of "oil particles". Lou et al. [50] simulated the diffusion traces of oil membrane on sea surface using Euler-Lagrange method, considering winds and evaporations. They considered the effects of turbulent diffusions, vertical diffusions, surface expansions, convections, evaporations and emulsifications to establish a 3D comprehensive model including spilling transports and destinations [51]. Yang [52] studied the oil transport diffusion of Beibu Gulf using the oil particle model but did not consider laying phenomenon. Li [53] established a three-dimensional dynamic model and applied it to the oil spilling simulation in the South China Sea. The distribution of flow field was simplified. Fannelop et al. [7], Hirst [14] and Doneker [54] also obtained experimental results of oil spilling traces. Their results have been cited broadly so that they can be used as benchmark solutions of oil spilling. Ghanmi et al. [55] studied the electromagnetic signature from sea surface covered by oil spills in bistatic case using the numerical forward-backward method. Some recent progresses can be found in [56][57][58][59].
As stated previously, few studies combined the oil diffusion on sea surface and rising process from the submarine pipeline. The simulation and prediction of the whole process of oil spilling has not been fulfilled. In this paper, the whole process of oil spilling from submarine pipeline to sea surface is studied under the combined effect of spilling speed, oil density and seawater velocity. The spilled oil can hardly be monitored by devices before it reaches to the free surface of the sea. Moreover, the oil boom floats on the free surface. Thus, the maximum horizontal migration distance (MHMD) is the most important quantity for preventing oil pollution. The correlations to quantitatively determine MHMD with corresponding time are proposed in this paper. The correlations can give reference values for the positions of the spilled oil droplets at any time and provide the basis for the emergency response strategy. Then, the transport of spilled oil on sea surface and corresponding pollution control method are studied. The systematic research connecting oil spilling from submarine pipeline and pollution control on sea surface is completed for the first time by numerical simulation. The results are meaningful for guiding the applications in petroleum engineering using the computer modeling approach.

Computational Models
Computational models to simulate the whole oil spilling process from submarine pipeline to sea surface are given below.

Governing Equations of Underwater Oil Transportation
Oil and seawater are considered as incompressible fluid. No phase change and slip are assumed on phase interfaces. Fluid flow is governed by the continuity equation and the Reynolds-Averaged-Navier-Stokes (RANS) equation as follows: where u i is the component of the mean velocity in the i (i = x, y, z) direction, u 0 i u 0 j is the Reynolds stress, x i is the coordinate, g i is the gravitational acceleration, t is the time, p is the mean pressure, q and t are the density and the kinetic viscosity respectively. The multiphase flow of spilled oil is described by the VOF. Fluid volume functions of water and oil (F w and F o ) are defined as follows: where v c , v o and v w are the cell volume, oil fraction and water fraction. The transport equations of F w and F o are: Density and viscosity can be evaluated by Eqs. (7) and (8): where the subscripts a, o, w represent air, oil and water. In this study, Reynolds number is in the range of 181.25~4531.25. At high Reynolds number, realizable k-ε model is used to describe turbulent flow: where where k and ε are turbulent kinetic energy and its dissipation, S ij is the strain rate, G k and G b are the turbulent kinetic energy caused by the mean velocity gradient and the buoyancy respectively, l is the dynamic viscosity of fluid, l t is the turbulent viscosity, Pr t is the Prandtl number (here takes the value 0.85). C l , C 1e , C 2 and C 3e are empirical parameters taking the values 0.09, 1.44, 1.9 and 0.9. r t and r e are turbulent Prandtl numbers with the values of 1.0 and 1.2 respectively.

Expansion of Spilled Oil on Sea Surface
After arriving sea surface, spilled oil forms a film. The oil film is thicker around the spilling area and acted by gravity, surface tension and viscous force. It is expanded to nearby area by these forces. At the initial stage of spilling, expansion is dominant, determining the area of the film and affecting evaporation and solution with increasing film area. Oil expansion is mainly affected by shear of seawater and turbulence. The former causes deformation of oil film. The latter affects the size of oil film.
Area of oil film in the process of expansion can be calculated as follows: where A o is the area of oil film, V o is the volume of spilled oil, K a is coefficient, R o is the radius of oil film, h o is the initial thickness of oil film (0.1 m in this paper).

Transport of Spilled Oil on Sea Surface
Oil film is driven by wind and water in the transport process. Thus, wind velocity and water velocity are two important parameters. Transport process includes convection and diffusion. Convection velocity can be calculated as follows: where U p and V p are convection velocity components in the x and y directions, U s and V s are seawater velocity components on sea surface in the x and y directions (U s ¼ V water sin h water , V s ¼ V water cos h water , h water is the direction of seawater velocity), C wind is the transport coefficient of wind ( = 0.03), U wind is the wind velocity at 10 m above sea surface, h wind is the angle of wind deflection, h is the angle of wind direction.
Turbulent diffusion can be expressed by the random movements of oil particles: where Da is the distance of turbulent diffusion in the a direction, K a is the turbulent diffusion coefficient in the a direction, R d is a random number in the range of [-1, 1], Dt is time step.

Evaporation of Spilled Oil on Sea Surface
Evaporation is the process that light components of oil turn from liquid to gas. When the spilled oil expands on sea surface, the area of oil film increases to promote evaporation. The mass of oil film losses in the evaporation, especially for light oil. Influential factors of evaporation include oil film area, seawater temperature, oil film thickness, oil components, wind velocity, etc. The evaporation rate can be calculated as follows: where N e i is the evaporation rate of the ith component, P i is the vapor pressure of the ith component, M i is the molecular weight of the ith component, q i is the density of the ith component, T is temperature, R is gas coefficient, K ei is calculated by the following expression: (24) where K is the evaporation coefficient, here takes 0.29, Sc i is the Schmidt number of the ith component.

Emulsion of Spilled Oil on Sea Surface
Spilled oil is in the action of winds and waves. The continuous oil film is broken by and mixed with seawater to form oil-water emulsion, which has much larger volume and viscosity than spilled oil. It can float on sea surface for long time. The emulsion usually appears several hours after spilling. At the initial stage, oil film is thick so that it is hard to be broken. When the spilled oil is expanded, the oil film becomes thin so that it is easy to be broken to form emulsion. The mass loss of oil film induced by emulsion can be calculated as follows: where D a and D b are the masses entering water and remaining in water respectively, l o is oil viscosity, c is surface tension. The rate for spilled oil returning from water to oil film is: Water holdup can be expressed by: where R a and R b are absorption rate and release rate of water in emulsification process, Y w is water holdup of emulsion, Y max w is the maximum water holdup of emulsion (usually 0.85), K a and K b are absorb coefficient and release coefficient of water, W a and A s are contents of wax and asphalt.

Solution of Spilled Oil on Sea Surface
Solution of spilled oil has weaker impact than other processes. It also causes mass loss of oil film. Solution rate can be calculated using following equations: where V di is the solution volume of the ith component, K di is the mass transfer coefficient of the ith component, X i is the molar fraction of the ith component, C s i is the solubility of the ith component, M i is the molar weight of the ith component, e i is relating with components of spilled oil (2.2 for aromatic hydrocarbon, 1.4 for alkane).

Numerical Methods
The above equations can be numerically solved using proper methods. Numerical methods for all the processes are introduced as follows.

Computational Domain and Boundary Conditions
Computational domain is shown in Fig. 1, including a submarine pipeline and seawater. The diffusion of the spilled oil on sea surface is also considered in this paper. The pipeline is on the seabed. A hole is on the top of the pipeline. Oil is spilling to the sea through the hole at a certain velocity and is diffusing when it reaches sea surface.
As shown in Fig. 1, boundary conditions are as follows. Logarithmic profile is used for seawater velocity as the inlet velocity distribution to fit the real conditions of the offshore flow. The equation of the profile is: where V wmax is the maximum velocity of seawater on free surface (usually V wmax = 0.1 m/s), H is the distance from the spilling hole to sea surface, D is the diameter of the pipe, 0 y H þ D.
The outlet boundary condition is set as static pressure distribution: The bottom of the domain is no-slip boundary condition (u = v = 0). This no-slip boundary condition is also used for pipe wall. The top of the domain is set as free surface boundary condition as it contacts air (p = 101325 Pa). A constant spilling speed v o is assumed vertical to the hole.

Mesh Generation and Numerical Methods
Unstructured triangular mesh is used for domain discretization. To obtain high quality mesh and minimize computational time, local refinement is taken near the hole. The mesh is gradually coarser for the area farther to the hole, as shown in Fig. 2. The computational time is reduced from 30 min to 10 min by using this local refinement compared with non-refinement unstructured mesh, showing the high efficiency of the computation.
Finite volume method (FVM) is used for the equation discretization. SIMIPLE algorithm is applied to deal with the pressure-velocity coupling. Second-order upwind and central difference schemes are used for the discretization of the convection and diffusion terms to ensure the stability and accuracy. The convergence criterion is that the residual of each equation is less than 10 -4 .

Numerical Methods for Expansion Process
The expansion process of oil film driven by surface tension and viscous stress is slow, so that explicit scheme is used to discretize Eq. (18): o is the volume of the spilled oil at the nth time step. Eq. (18) also has the analytical solution by using the following integration: The analytical solution can be obtained as:

Numerical Methods for Transport Process
Euler-Lagrange method is used to simulate the transport process of spilled oil. Flow simulation on sea surface is on the fixed Euler grid. Lagrange method is used to track oil film transport. Lagrange velocity of oil film are: where U n ð Þ p;i and V n p;i are the velocity components of oil film in the x and y directions, U n ð Þ s;i and V n s;i are the velocity components of seawater in the x and y directions respectively, U n wind;i is the wind velocity in the x direction, the subscript i represent the ith point, the superscript n represent the nth time step. Based on the velocities of oil film, positions of oil film at any time can be calculated as: where x nþ1 ð Þ p;iþ1 and y nþ1 ð Þ p;iþ1 are positions of oil film in the x and y directions respectively. When the velocity on sea surface is large, the transport process also needs consider turbulent diffusion. Eq. (22) can be calculated as: where Dx nþ1 ð Þ turbulent and Dy nþ1 ð Þ turbulent are displacements of oil film in the x and y directions respectively, K x and K y are turbulent diffusion coefficients in the x and y directions. Final positions of oil film considering turbulence are: where x nþ1 ð Þ pt;iþ1 and y nþ1 ð Þ pt;iþ1 are final positions at the (n+1)th time step in the x and y directions respectively.

Numerical Methods for Evaporation Process
The time step in the process of transport is short so that evaporation can be assumed as flash process, i.e., phase change of oil film only occurs at every moment. Evaporation inside the time steps is neglected. Thus, from Eqs. (23) and (24), instantaneous evaporation rate can be calculated as: where N n ð Þ et is the total evaporation of oil film, P n ð Þ i is the vapor pressure of the ith component, T n ð Þ is temperature at the time n.

Numerical Methods for Emulsion Process
Similarly, emulsion is also assumed to happen only on snapshots: where D n ð Þ is the total mass loss of oil film by emulsion at the time n, h n s is the height of oil film.

Numerical Methods for Solution Process
Solution is also assumed to flash, so that the instantaneous solution volume can be calculated as: where V n ð Þ dt is the total solution volume at the time n.

Model Verification
To verify the above models and numerical methods, our numerical results are compared with experimental data of previous researchers as benchmark solutions [6,13,53]. Fannelop et al. [7] and Hirst [14] all did experiments of oil spilling traces at R 0 ¼ 8. The difference is the former did it at F r 0 ¼ 10 compared with the latter at F r 0 ¼ 20. R 0 is the velocity ratio between buoyant jet and surrounding fluid velocity: where V j ! is the initial velocity of the spilling hole, V j ! is the velocity of surrounding seawater. F r 0 is the Froude number: where V j ! represents the momentum of the jet, g 0 represents the buoyancy effect: First, our numerical results of spilling traces are compared with their experimental data under different Froude number at the same velocity ratio. The comparison results are shown in Fig. 3. The agreement of numerical results and experimental data of oil spilling traces is quite well. The relative errors are listed in Tabs. 1 and 2. Compared with Fannelop et al.'s results [7], the maximum error is 11% while the minimum error is only 0.7%. Average error is about 5.5%. Compared with Hirst's results [14], the maximum error is 10.9% while the minimum error is only 0.9%. Average error is about 5.2%. To further verify the numerical results with experimental data at different R 0 , we compare them in Fig. 4. Experimental data are Doneker's results [54]. Oil spilling traces agree well at R 0 ¼ 5 and R 0 ¼ 10.      The numerical results of oil expansion on sea surface are compared with the analytical results from Eq. (38) to verify the numerical methods. Fig. 5 shows that they agree well with each other, indicating the numerical methods in this paper has high precision.
In this paper, a typical submarine pipeline in the largest marine oil field in China-Bohai oil field is selected as the research object, because this pipeline is longer with higher accident risk. Moreover, the population near the pipeline is high so that the influence of oil spilling is severe. In this region, average depth of water is 15 m. Thus, the size of the computational domain is selected as 20 m-length and 15 mheight. Larger length and closer position of the pipe to the inlet are helpful for studying the development of the spilled oil on downstream. The diameter of the pipe is 0.6 m. The distance of the pipe center to the domain inlet is 1.5 m. The diameter of the hole on the top of the pipe is 0.06 m. The density of the seawater is 1025 kg/m 3 . The dynamic viscosities of the crude oil and seawater are 0.048 Pa·s and 1.8 × 10 -5 Pa·s. 8 cases are designed and simulated under different oil density, spilling speed and water velocity, which are shown in Tab. 5.

Standard Case
Case 2 is selected as the standard case (oil density 900 kg/m 3 , spilling speed 5 m/s, water velocity 0.1 m/s). Fig. 6 is the distribution of the two phases (oil and water) with time. It represents the oil spilling process. Different colors represent volume fractions of oil and water. Continuous oil stream flows out of the small  leakage hole. It becomes discontinuous in the action of turbulence of seawater. Thus, the spilled oil exists in the form of droplets, which are acted by the gravity, inertial force, buoyancy force and shear stress. The droplets become more scattered when they are rising up. With the decreasing depth of the droplets, the migration distance increases apparently. High speed water current offers large shear stress on oil stream so that it is accelerated in the horizontal direction. On the other hand, gravity provides the resistance of rising movement so that the oil stream is decelerated in the vertical direction. Therefore, the migration distance in the horizontal direction is much larger than in the vertical direction in the same time scope. Further analyzing the results, we can find that the first oil droplet takes about 19.8 s moving to the free surface. After that, oil droplets continuously move to the downstream by the shear stress of water current. At about 24.6 s after spilling, the migration distance in the horizontal direction reaches to its maximum value (MHMD) 8.3 m. At this moment, the ratio of the horizontal distance and the vertical distance is 0.59.
In the following sections, oil density, spilling speed and water velocity are changed based on the standard case to study their effects on the trace of spilling oil.

Effect of Oil Density on MHMD
In this section, all parameters are the same as the standard case except oil density to study its influence on the MHMD with corresponding time. Fig. 7 shows the process of oil spilling from the submarine pipeline to free surface at three different oil densities. At each moment (4 s, 12 s, 20 s), lower-density oil rises up higher than higher-density oil. Oil droplet is acted by buoyancy and gravity in the vertical direction at the same time. For the droplets with the same diameter, their buoyancies are the same, but the droplet with higher density has larger gravity so that it rises slowly. The droplet with lower density rises faster and enters high-speed water layer earlier, so that it is driven by larger velocity of seawater. Thus, the lower oil density corresponds to the shorter time to free surface. The MHMD of the oil with density 950 kg/m 3 is 8.4 m while the MHMD of the oil with density 850 kg/m 3 is 8.6 m. However, the time to achieve the MHMD for the former oil is 1.34 times of the latter oil. Thus, oil density can be considered to influence the MHMD very slightly but influence the time to achieve MHMD significantly. Lower-density oil can achieve free surface much faster, but the positions of the oil boom for different oils are almost the same.

Effect of Spilling Speed on MHMD
Spilling speed is an important factor to affect the diffusion of the spilled oil. The process of oil spilling under different spilling velocities is shown in Fig. 8. When the water velocity is 0.1 m/s, the shortest time of oil arriving sea surface is 24.0 s while the MHMD is 7.8 m at the lower spilling speed (v o = 3 m/s in Fig. 8a), but the shortest time becomes 19.8 s while the MHMD becomes 8.3 m at higher spilling speed (v o = 5 m/s in Fig. 8b). The reason is that the higher spilling speed leads to larger upper momentum and more scattered diffusion of oil. To reduce the environmental impact, fast responses should be made for high speed spilling.

Effect of Seawater Velocity on MHMD
Seawater is the carrier of spilled oil so that it is important for oil diffusion. In this section, we change the seawater velocity to evaluate its impact on oil spilling. The process of oil spilling under different seawater velocities is shown in Fig. 9. Larger velocity provides more shear stress and kinetic energy to spilled oil so that the oil trace shifts to the downstream further. The MHMD at lower seawater velocity (v wmax = 0.1 m/s in Fig. 9a) is 8.3 m, which is 1.1 m shorter than that at higher seawater velocity (v wmax = 0.2 m/s in Fig. 9b). Corresponding time is 19.8 s and 20.8 s respectively. The time depends on the trace of the spilled oil, i.e., d ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi h 2 þ l 2 p (d, h, l are total displacement, vertical displacement, horizontal displacement respectively). The vertical displacement is constant since the pipeline and the hole are not moving, i.e., h = H. The total displacement only depends on the horizontal displacement. For the case v wmax = 0.1 m/s, the horizontal displacement is shorter than the case v wmax = 0.2 m/s, so that the total displacement is shorter leading to shorter rising time. Thus, lower velocity of seawater needs shorter response time while higher velocity of seawater corresponds to longer oil boom.

Correlation of MHMD and Spilling Time
Effects of oil density, spilling speed and seawater velocity on spilling trace are studied above. However, for engineering application, a general correlation is needed to fast predict more situations so that quick responses can be taken by engineers. Thus, it is necessary to propose the correlations of MHMD and corresponding time.
Through the computations of oil spilling using the 5 densities in Tab. 1 (750 kg/m 3 , 800 kg/m 3 , 850 kg/ m 3 , 900 kg/m 3 , 950 kg/m 3 ), one can obtain the dimensionless distance (L max /H) under different oil densities and the correlation of the dimensionless distance with oil density (Eq. (53)). The data and correlation are  Similarly, through the computations of oil spilling using the 5 densities in Tab. 1, one can obtain the dimensionless time (v o t/H) to achieve the MHMD under different oil densities and the correlation of the dimensionless time with oil density (Eq. (54)), as shown in Fig. 11. The coefficient of the correlation is very high (R 2 = 0.9994), so that the correlating precision is enough.

Dynamics of Spilled Oil on Sea Surface
Once the spilled oil reaches sea surface, it will transport and expand to downstream. The dynamics of spilled oil on sea surface directly affect the pollution control strategy. Thus, it is important to study oil spilling process on sea surface. It is simulated under different water velocities and wind velocities in 10 km × 10 km region. The initial distribution of spilled oil is shown in Fig. 12, where the blue region represents spilled oil on sea surface.
Wind velocity at 10 m above sea surface is assumed to be 10 m/s. Seawater velocity is 0.4 m/s. Their directions are both Northeast 45°. Through simulation, oil distributions at different time (1 h, 3 h, 5 h, 7 h) are shown in Fig. 6. Initial position of spilling center is about (1.3 km, 1.3 km). Spilled oil is moving to the Northeast and arrives (8.6 km, 8.6 km) after 7 h. At the time of the 7 h transport, range of the spilled oil film expands apparently. Oil edge tends to disperse instead of continuous at the initial time. The reasons are randomly moving oil particles in turbulence and thinner oil film, which cause oil film is not as tight as initial time. To study effects of water velocity and wind velocity on the transport of spilled oil, simulations are also made at other velocities. Fig. 14 shows the instantaneous distribution of spilled oil at wind velocity 10 m/s and water velocity 0.8 m/s. Comparing Fig. 14 with Fig. 13, we can find that the center position of spilled oil is (6.6 km, 6.6 km) at 5 h when the water velocity is 0.4 m/s. It moves to (7.8 km, 7.8 km) at the same time for water velocity 0.8 m/s, which is close to the position (8.5 km, 8.5 km) at 7 h for water velocity 0.4 m/s. However, the spilling center is out of the computational domain at 7 h for water velocity 0.8 m/s. Thus, variation of water velocity changes the transporting speed of spilled oil apparently at the same wind velocity. center is out of the domain at 5 h for wind velocity 15 m/s. Therefore, wind velocity also apparently affects the transport speed of spilled oil at the same water velocity. Larger wind velocity corresponds to larger transport speed so that the displacement of oil film is larger at the same period.
Through the above simulation, spilling range extends continuously in the transport of oil. The edge of oil film becomes thinner and scattered. Evaporation, emulsion and solution have small impact for transport at the initial stage. The transport of spilled oil is dominated by water velocity and wind velocity. Larger velocities of water and wind leads to faster transport so that shorter response time is needed to dealing with spilling accidents.

Spilling Traces in Non-Uniform Flow Field
For the oil spilling on sea surface, the flow field is usually complex. Spilling traces in different flow fields are simulated in this section. Fig. 14 shows two different flow fields in the computational domain. In the non-uniform flow field I (Fig. 16a), water velocity has different directions with the same magnitude for different positions. In the non-uniform flow field II (Fig. 16b), directions and magnitude of water velocity are all different for different positions.  From above analyses, we know that water velocity is a key factor for oil spilling. When water flushes the oil boom and its velocity exceeds 0.386 m/s, spilled oil will escape from the oil boom. To prevent this problem, oil boom should be placed not vertically to water flow direction. Thus, the maximum angle of oil boom can be calculated as follows (Fig. 19): where α is the angle between water velocity and oil boom, V lim is the maximum velocity component vertical to oil boom. The ultimate angles under different water velocities can be calculated according to Eq. (55). The results are shown in Fig. 20. With increasing water velocity, the angle of the oil boom decreases dramatically but the decreasing is more and more slowly.
Oil boom should be restricted not only by the angle but also by the pulling force because too large pulling force may cause the oil boom broken. The maximum pulling force of oil boom acted by water can be calculated as follows: where F is the maximum pulling force of oil boom, L is the length of oil boom, H is the height of the oil boom (usually 0.3 m), g is gravitational acceleration. Oil boom is usually semi-circle as shown in Fig. 21. Its diameter should be at least larger than the widest dimension of the spilled oil. Thus, the length of the oil boom can be calculated as: where l o is the width of the spilled oil.
From Figs. 13d, 14c and 15b, the width of the spilled oil is about 0.5 km. The length of the oil boom calculated by Eq. (56) is about 785 m. When the water velocity is 0.4 m/s (Fig. 13d), the angle α calculated from Eq. (55) is 75°while the maximum pulling force of oil boom in Eq. (56) is 36926 N. When the water velocity is 0.8 m/s (Figs. 14c and 15b), α is 29°while the maximum pulling force is 147706 N. Therefore, larger water velocity leads to much larger pulling force, requiring higher strength of oil boom.

Conclusions
The whole process of oil spilling is numerically studied, including underwater spilling from the pipeline and movement of spilled oil on sea surface, which can give valuable reference for marine oil pollution control in engineering. Oil spilling and diffusion process is simulated in this paper. The effects of oil density, spilling speed and seawater velocity on spilling trace of oil are analyzed. The maximum horizontal migration distance (MHMD) with corresponding time are obtained via the simulation. Correlations of MHMD and its time are proposed. Dynamic factors (expansion, convection, turbulent diffusion) and non-dynamic factors (evaporation, emulsion, solution) are considered in the numerical simulation. Governing equations based on oil particle dynamics are solved using Euler-Lagrange method. Based on these numerical results and discussions, conclusions can be made as follows: 1. Oil is continuously flowing out of the spilling hole, but continuous oil flow is interrupted by the turbulence of seawater to form discontinuous oil droplets. The first droplet arriving sea surface does not achieve the MHMD. It needs further movement. For the standard case, the MHMD is 8.3 m with corresponding time 24.6 s. 2. Higher density of oil leads to longer time to achieve MHMD, but the MHMDs for different densities have much smaller difference. Thus, oil density mainly affects the time to achieve MHMD. Higher spilling speed provides larger rising energy of oil, so that the MHMD is larger with shorter corresponding time. Higher velocity of seawater puts larger shear stress on oil droplets, so that the MHMD is larger. 3. Based on numerical results, correlations for MHMD and corresponding time are fitted. For parameters different from the simulation, the MHMD and time can be calculated quickly according to the two correlations to provide reference information for emergent responses in engineering. 4. The transport speed and direction are also related with the local velocity in flow field. The impact of flow field on the prediction of spilling position is very important. For large wind velocity, large water velocity and complex flow field, transport of spilled oil may very fast with variable directions. Quick response is needed to prevent larger pollution area. Evaporation, emulsion and solution have limited impact on oil spilling. They mainly reduce the volume of oil film and have few influences on direction and displacement of spilling transport. During the initial serval hours of transport, oil particles move randomly by turbulence. Oil film becomes thinner and larger with scattered edge.
Funding Statement: This research was funded by National Natural Science Foundation of China (NSFC), grant number No. 51576210.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.