Research centers and Brazilian oil industry have reported costs of the order of U$90.000,00 per cubic meter of oil as a result of leakage . Thus, in Brazil, the environmental crime law, one of the most stringents in the world, has like objective to assess environmental disasters and applying strict fines and possible arrest orders as necessary. With such requirements, the petroleum industry must achieve high safety standards in its operations. For this reason, it is crucial compliance with technical standards related to facilities and other legal requirements required in order to improve the performance of production units.
In present day, oil transportation has occurred by pipelines. Submerged pipes are many kilometers long and are intended for the transport of oil and gas. Then, submerged ducts and floating platforms generally tend to interact with waves and water stream resulting in the pressure differences around the outer surface that can cause vibrations in the pipeline  . Pressure distribution around of the pipe is subjected to variations which may cause alteration in the force components at the pipe surface, that is, it changes the drag force parallel to flow direction. For this reason, these forces can produce vortex induced vibration in the tube, which is an important variable for oil flow and risks studies in marine environmental .
In oil transportation by pipelines, there are pressure variations that cause different fluid flow velocities and abrupt variations in area, which can cause pipeline ruptures and leaks  .
In the event of an accidental rupture, a quick and adequate actions is necessary, in order to reduce the environmental consequences. As an example, with the use of containment barriers, is possible to control the spread of oil in the sea. For this reason, it is of fundamental importance to understand the fluid flow behavior in the marine environment and the dispersion of the oil plume, with the aim to predict the leakage location   .
According to ANP (National Agency of Petroleum, Natural gas and Biofuel) which is the Brazilian institution responsible to regulate the activities related to the oil natural gas and biofuels industries in Brazil, 482 pipelines are authorized to transportation of hydrocarbons from refineries to consumption centers, and in general they are subject to leaks resulting from rupture caused by different mechanisms such as, fatigue or corrosion .
Leakage is considered extremely dangerous mainly in view of the environmental implications, especially in subsea pipelines and in remote areas .
When leaking in submerged pipes, oil tends to flow through the orifice towards the sea surface as a result of the difference in fluid densities. Thus, the water that is located in the external surface of the pipe tends to penetrate inside the pipe, occupying all the space left by the leaked oil. This process of exiting oil and entering water in the pipe is called advective migration .
For accurate leakage detection and identification, it is necessary to know the behavior of the fluids inside the pipeline, which allows determining the pressure drop along the pipeline    , and the total amount of fluid (oil, gas or even contaminated water) released, and allows to stop the pumps and the control valves avoiding or minimizing environmental damage.
The establishment of flow patterns is common in multiphase flow and it is related to the configuration of the phases present in the system. This physical situation allows that a force balance on each phase interfaces be established. This force balance depends on the parameters related to pipe geometry, for example the angle of inclination, diameter and roughness of the duct wall fluid flow (for example speed of each of the phases involved in the system) and the fluid (density, viscosity and surface tension) .
The understanding of interfacial phenomena of immiscible liquids is of great importance for the study of the advective migration process in submarine pipelines .
This information allows the prediction of flow behavior, the geometric shape of the water-oil interface and provides a great phenomenological foundation with respect to the theories of perturbation and stability criteria, and mathematical modeling.
Activities related to drilling and oil production on the high seas, oil transportation in load ships, and pipelines submerged, increase the potential for oil leakage. In this sense, several researches have focused on the study of flow behavior in the pipe with leakage. Those papers, investigate and quantifying the movements of oil leakage on the water surface. Activities related to drilling and oil production on the high seas, oil transportation in load ships, and pipelines submerged, increase the potential for oil leakage. Those papers, investigate and quantifying the movements of oil leakage on the water surface. That said, oil leakage models and numerical simulation systems were developed to predict the process behavior   .
However, few works are addressed to the spread of oil in water. Therefore, various factors affecting the spread and drift of oil leaks in the sea, such as: physical characteristics (density and viscosity), physical environment of the sea (wave, wind, current) and leakage hole parameters (leakage diameter, pressure and speed)    . In complement, this study uses computational fluid dynamics to understand important details about the oil leak in a pipeline, as the flow patterns (in the vicinity of the leakage), pressure fields, velocity, volumetric fraction of the involved phases and to evaluate the effects of the velocity of entry into the tube and into the sea and the outlet pressure in the tube on the behavior of fluids in the region of the leakage. Hence, this work aims to study numerically the oil dispersion in the seawater, originated from a subsea pipeline leakage.
2.1. Physical Domain and Geometry
The physical domain consists of a submerged pipe in the sea (two-dimensional domain with 20 m length and 10 m high) as illustrated in Figure 1. In this system, oil flows inside the pipe and water flows in the region above the pipe, which is considered the external medium. The pipe have the following dimension: 0.2 m diameter, 0.02 m wall thickness and a hole of 0.02 m diameter located 10 m from the pipe input section.
2.2. Computational Domain
The representative computational mesh of the physical domain was builded by using the ICEM-CFD software. Figure 2 shows the generated numerical mesh,
Figure 1. Physical domain analyzed.
Figure 2. Numerical mesh of the studie physical domain: (a) whole domain, (b) Central region and (c) Leakage region.
with 139,488 hexahedral elements and 93,666 nodal points, highlighting the piping and leakage regions.
2.3. Mathematical Modeling
2.3.1. Governing Equations
The free surface model was adopted in this study in order to analyze the water-oil two phase-flow, according to the following considerations:
1) Incompressible, transient and isothermal flow;
2) Fully developed laminar fluid flow inside the pipe;
3) Turbulent fluid flow at the sea domain;
4) Two-dimensional domain;
5) No mass transfer between the phases;
6) No mass generation;
7) Physico-chemical properties of fluids are constant;
8) Gravitational and drag forces were considered;
9) Smooth pipe (null roughness).
· Mass conservation equation
The following mass conservation equation has been used:
where, ρ, f and are, respectively; the specific mass, volumetric fraction and the velocity vector of the α phase. represents the mass source term, and the term of mass flow rate per unit volume of the α phase to β phase. The sub-indices α and β correspond to the phases involved in the multiphase flow.
Disregarding the mass source term and the mass flowrate term per unit volume of the phase β for phase α, , the mass conservation equation for the fluid mixture is given as follows:
· Linear momentum conservation equation
Regarding the conservation of linear momentum, the following equation has been used:
where µ is dynamic viscosity, P is the pressure, is the term that represents the external forces per volume unit; is the term that represents the momentum transfer induced by the mass transfer, and is the term that describes the total interfacial force per volume unit of the α phase due to the interaction with the β phase at the interface (drag and gravitational forces).
By using the adopted assumptions Equation (3) can be rewritten as follows:
The closure equations applied to the conservation equations are used at different manners. For the continuity equation, the sum of the volumetric fractions of the involved phases is equal to 1. Concerning to the equation of momentum transfer, the closing is done by the force acting between the phases, expressed by the drag coefficient .
The interfacial drag force is given by:
where = 0.44 represents the drag coefficient and = 0.001 m is the fluid (oil) particle diameter.
· Turbulence equation
Among the different turbulence models used for the treatment of turbulent flows, what stands out most is the k-ε model. The k-ε turbulence model is a model where the Reynolds tensors are proportional to the mean velocity gradients. The proportionality constant being identified by a turbulent viscosity . This model is defined by the transport equations applied to turbulent kinetic energy and the dissipation rate of the turbulent kinetic energy, which are solved for are turbulent length and time scale  .
For the turbulent kinetic energy, k, and the turbulent dissipation rate, ε, we can write:
where Gα is the kinetic energy generation within the α phase; C1, C2, σk and σε are constants equal to 1.44, 1.92, 1.0, and 1.3, respectively.
The term εα corresponds to the dissipation rate of turbulent kinetic energy, and kα corresponds to the turbulent energy, which are defined as follows:
where lα is the spatial scale length, qα is the velocity scale, and , is calculated using the following expression:
where Cα is an empirical constant, and μ is the dynamic viscosity of the α phase.
The turbulent viscosity is defined as follows:
In this work, the constant Cμ was considered to be equal to 0.09.
2.3.2. Initial and Boundary Conditions
For this propose the fully developed laminar single-phase flow of oil into the pipeline was considered. The initial and boundary conditions used in the simulations are shown in Figure 3.
At the pipe inlet the oil velocity profile is given as follows:
For the external domain (sea water), to it was considered water velocity according to Equation (13):
where Ue,m is the input velocity at the domain, Umax is the maximum velocity of oil in the pipe.
Further, it was considered a reference pressure equal to 1 atm.
2.3.3. Fluid Properties and Studies Cases
2.3.4. Numerical Solution Procedure
In the numerical solution of the governing equations was used the software Ansys CFX. Table 3 sumarize all numerical treatment used in each simulation.
3. Results and Discussions
3.1. Inlet Oil Velocity Effect
As reported by the relevant literature, larger pressure gradients occur proximity
Figure 3. Initial and boundary conditions adopted for pipeline and sea.
Table 1. Physical-chemical properties of fluids used in numerical simulation.
Source: *  ** , ***  .
Table 2. Analized cases in this work.
Table 3. Considerations used in the simulations.
to the leakage, that is, low variations in the pressure values are observed before and after the leakage occurs.
Figure 4 shows the behavior of the total pressure drop as a function of time in two transversal planes, located 1 m before the leakage, for the cases 1, 2, 3 and 4.
Figure 4. (a) Total pressure variation as process time and (b) parameter measurement position P11 (Cases 1 - 4, Position P11 and P13).
From the analysis of this figure, it can be seen that the total pressure increases as the oil velocity increases. Afterward, there is a greater loss of mechanical energy due to viscous friction with the pipe walls. In the initial seconds, there is a strong decline in pressure drop, being more intense increasing the inlet oil velocity in the pipe.
In agreement with the literature  , the time required for the pressure behavior to reach stability after leakage depends on the diameter of the fluid escape orifice in the pipe.
Figure 5 represents the pressure drop as a function of the position along the pipe in Y = 0.1 m with the presence of a leakage in the pipe. By analyzing this figure, in 5 s time, the total pressure drop shows a linear decreasing behavior
Figure 5. Total pressure drop as a function of axial position in the pipe (t = 5 s) and (b) measurement position indication (Cases 1 - 4, position P2.)
with a axial position before and after the leakage, according to the increase in the oil velocity at the inlet of the pipe. Higher inlet velocity the higher the pressure drop in the region of the leakage.
Figure 6 shows the behavior of the total pressure drop along the pipeline position P2, for the physical situation with and without leakage. By analyzing this figure we can see that there is a decreasing linear behavior of this parameter for the situation without leakage and a pressure discontinuous behavior for the simulation with leakage. In the six meters before the leak, a pressure drop of approximately 700 Pa is observed, clearly indicating the presence of leak due to the loss of fluid mass. After the leak, there is equilibrium state in the flow leading to a linear decrease in pressure until the pipeline outlet.
Figure 7 illustrates the surface water velocity distribution across the domain (tube and sea water), with detail of the region of the leakage, for the case 4, whore the inlet oil velocity is 2.16 m/s. By analyzing this figure we can see that, when the oil enters the marine domain, in the form of a jet resulting from the
Figure 6. (a) Total pressure drop along the pipe and (b) measurement position indication (Case 4, position P2).
Figure 7. Surface water velocity field at the domain (pipe and seawater) and detail of the oil jet in the sea region: (a) velocity field, (b) velocity vector field (Case 4).
pressure difference between the oil at the leakage and sea domain, its surface velocity reaches a value of 0.299 m/s. The vectors indicated in the figure do not represent greatness, but only the flow direction of the phases .
An experimental behavior similar to the result discussed earlier was observed in Li et al.  with different leakage rate velocities. The authors aimed to observe the process of oil droplet formation from experiments that were carried out according to the increase in the oil jet speed, so, they kept other parameters constant to increase the number of Reynolds between the phases.
The experimental images show the oil released into the water at different jet speeds from a 1.95 mm nozzle using a high resolution camera (Figure 8). In the 8a image, a stable laminar jet can be seen, once it reaches a certain height, it dissolves into individual drops. In the 8b and 8c images, it is observed that the jet has a totally turbulent behavior, being intensely presented in a large number of smaller droplets that shows a dispersed behavior.
Zhu et al.  proposed to better examine the underwater spread of the oil leak in an experimental apparatus and to assess the effects of shear flow and pressure on leakage local. A constant shear flow was generated with the average speed ranging from 8.6 cm/s to 22.0 cm/s. The schemes of the experimental configuration consist of a channel 1.5 m long and 0.1 m wide by 0.3 m high. The water depth is defined as a variable combined with the shear flow velocity. A drop of oil is initially deformed and elongated due to turbulence of the running torn in parts by the cutting force of the water flow, this is, therefore, the first way to break the jet.
Figure 9 shows the profile of the oil volume fraction, at different times (x = 10 m, case 4). After an analysis of this figure it is observed for the time equal to zero,
Figure 8. Images of oil discharged into water from a 1.95 mm nozzle. Source. Li et al. .
Figure 9. Distribution of the oil volumetric fraction in the leakage for different times (Uo = 2.16 m/s) and (b) parameter measurement positions P8 (Case 4).
that the oil volume fraction in the orifice is equal to 1, due to the initial boundary condition adopted. With the course of time, changes in the behavior of the oil volumetric fraction is observed along the leakage orifice.
Furthermore, it is verified that this parameter presents null value in some position, that suggest the presence of water in this location.
Figure 10 represents the profile of the oil volumetric fraction in the tube, at the position P14 (y = 0.198 m), located 2 mm from the leakage orifice, and, the volumetric fraction of water in the sea stream, at the position P15 (y = 0.207 m),
Figure 10. (a) Profile of the oil volumetric fraction in the tube at the position P14, (b) Axial profile of the water volumetric fraction at the sea at the position P15 and (c) measurement positions indication (position P14 and P15, t = 10 s, Case 4).
located 2 mm from the leakage orifice, for a time of 10 s, in the presence of leakage (case 4).
It can be seen in Figure 10(a) that, for higher inlet oil velocity in the tube (2.16 m/s), there is a reduction in the oil volumetric fraction from 1 to 0.8008, (20%), implying that the high losses of oil fraction in the pipeline occurs in the regions closest to the leakage.
Figure 10(b) shows that, for standing water (0 m/s), the volumetric fraction declines slightly from 1 to 0.99997575 (0.01%).
Figure 11 shows the behavior of the oil volume fraction of oil in the sea stream, 1 m before, in the center and 1 m after the leakage in a time of 10 s (case 4). By analyzing this figure, it can be observed that, for the position in which the leakage is, a peak of oil volume fraction, indicating the presence of oil in the seawater. For the positions before and after the leak, the fraction of oil is equal to zero, which indicates that there is no oil flow in the sea stream for these positions and this instant of time. This behavior occurred by there physical conditions: 1) the water velocity is zero, 2) the oil phase is less dense thar the water, and 3) pressure gradient in the vertical direction.
3.2. Pipe Output Pressure Effect
Figure 12 exhibits, graphically, the volumetric fraction profile of the oil, in the pipeline, at the position P12, located at the center of the leakage orifice, for different pressures at the pipe outlet, in t = 10 s (cases 2, 5, 6 and 7).
At this time, it can be observed a volume of oil equal to 1 at this region of the leakage, and posteriorly suffers an abrupt reduction in from 1 to 0.82 approximately, according to the conditions imposed in the duct outlet. The use of the prescribed pressure condition provides more consistent results compared to those available in the literature, as stated in . This behavior was verified in all outlet pressure of the pipe. Small variations can be observed sea the leakage.
Figure 13 shows the oil superficial velocity along the pipeline, for an oil velocity
Figure 11. Transverse profile of the oil volumetric fraction of in three positions, at the sea domain and (b) measurement positions indication (position P7, P8 and P9, case 4, t = 10 s).
at the pipe inlet of 1.08 m/s (case 6) and time equal to 5 s, at the positions P1, P2 and P3, in pipe.
The results indicate that the oil velocity at the center of the pipe reaches the maximum velocity position P2, with a value equal to 1.08 m/s.
Figure 12. (a) Transverse profile of the oil volumetric fraction the leakage region, for different prescribed out let pressures (t = 10 s) and (b) measurement indication (Position P12, Cases 2, 5, 6 and 7).
At the positions P1 and P3, in the vicinity of the pipe walls, the oil reaches a velocity of 0.65 m/s, which can be explained due to the zero-velocity condition adopted in the pipe walls, that is, the closer to the pipe wall the lower the oil superficial velocity. It can be seen that the oil has higher velocities at the center of the pipe.
Figure 14 represents the behavior of the oil superfical velocity as a function of the position along the pipe for different outlet pressures of the pipe (Pst), at the time of 5 s, and position P2, in the presence of leakage.
From these results, it can be observed, in general, that the prescribed pressure boundary conditions at the outlet have no significant influence on the behavior of the oil flowing inside the pipe.
Figure 13. (a) Oil superficial oil velocity along the pipe at different positions and (b) measurement positions indication (Position P1, P2 and P3).
It should be noted that, after the first 2.5 m of the pipe inlet, the oil superficical velocity is reduced from approximately 1.08 to 0.93, which can be explained by the inlet effects in then it remains practically constant, indicating a tendency stabilization.
Figure 15 represents a 7 meter stretch of the oil volume fraction field located 9 m from the inlet section, 4 m from the outlet and 3.5 m high from the analyzed physical domain, for different prescribed pressures at the pipe outlet, oil inlet velocity equal to 1.08 m/s and standing water for the seawater stream, at 10 seconds of process and in the presence of a leakage (case 2, 5, 6 and 7). From the analysis of this figure, in can be observed that the volumetric fraction of oil practically prevents the same behavior for the different pressures adopted at the outlet of the pipeline.
Figure 14. Oil superficial velocity along the pipe: (a) for different pressures at the pipe outlet (t = 5 s) and (b) measurement position indication (Position P2, cases 2, 5, 6 and 7).
According to Tavares  who evaluated the formation and dispersion of an oil plume in underwater from a leakage originated of in a small diameter orifice in a catenary riser, there is tendency of the oil plume to slightly incline to the right side, as observed Figure 14. This behavior can be explained by dragging the oil plume through the water stream and the difference in density between the phases (water and oil).
It was also possible to observe that the oil plume, in the established conditions is able to reach the water surface, showing spreading (upstream and downstream) in the flow direction of sea water.
Figure 16 shows the transverse profile of the volumetric fraction of the oil inside the pipeline in three axial positions (1 m before and after the orifice, and at the center of the leakage orifice in the pipeline) in the instant 5 s (case 6), in the presence of the leakage.
Figure 15. Distribution of the oil volumetric fraction for different pressures at the outlet of the pipe (t = 10 s, Case 2, 5, 6 and 7).
Figure 16. Transverse profile radial distribution of the oil volumetric fraction in three positions in the pipe (t = 5 s) and (b) measurement position indication (Position P11, P12 and P13).
By analyzing this figure, it is possible to verify that, in the region before the leak, located 9 m from the inlet section of the pipeline, the volumetric fraction of the oil phase remains constant and equal to 1. Thus, indicating that, for this situation, the disturbance of the leakage does not influence the behavior of the oil flowing inside the pipe, possibly due to the viscous resistance effect.
For the position located 10 m from the inlet of the pipe (center of the leakage orifice), the fraction decreases, indicating the entry of water into the pipe, consequently generation oil-water mixture zone. At the location 11 m from the pipe entrance, a value for this parameter of approximately 0.775 is noted, are to the difference in density between water and oil.
3.3. Inlet Seawater Velocity Effect
Figure 17 shows the profile of the water and oil superficial velocity in the sea stream for different inlet velocity of water in the sea, at the position 4, in a time instant of 10 s in the presence of leakage. In Figure 11(a), it can be seen that, in the region of leakage (x = 10 m) there is an increase in the water superficical velocity for each value of inlet water velocity.
For higher velocity (1.08 m/s), the superficical velocity reaches a value of approximately 4 m/s.
According to Zhu et al. , the higher the velocity of the water stream, the oil reachs the maximum horizontal distance from the leakage orificie when it reaches the water surface. In Figure 11(b), it is possible to observe that the oil presents low superficical velocities in the sea stream. Thus, we can states that there results are due to the buoyancy effect.
Figure 18 represents the behavior of the total pressure drop of the water in the sea stream (positions P4, P5, and P6), in the time instant of 10 s, for standing water (case 8) in the presence of a leakage. From the figure, it is observed that, for the position P4 (y = 1 m) a greater total pressure drop due to the high depth and the small distance to the leakage orificie. This is explained by the gravity effect which influences the pressure exerted by the water and the along oil volume fraction this line, close the leakage orificie.
Figure 17. Superficial velocity profile along the seawater for different inlet water velocities (a) water and (b) oil, (t = 10 s) and (c) measurement position indication (Position P4, Case 13).
Figure 18. Total pressure drop along the sea stream different axial positions (t = 10 s), and (b) measurement positions indication (Position P4, P5, and P6, Case 8).
Figure 19 shows the profile of the water volumetric fraction at different positions in the sea stream, for a time instant of 5 s in the presence of a leakage (case 8). From the analysis of this figure, it is observed that, for positions P7 and P9,
Figure 19. Transverse profile of the water volumetric fraction at different positions in the sea stream (t = 5 s), and (b) measurement position indication (Positions P7, P8, and P9, Case 8).
located 1 m before and after the leak, the water volumetric fraction in the sea remains constant and equal to 1, indicating that, for these positions, there is no oil flowing within the sea stream. For position P8, located in the pipeline center, the water volumetric fraction decreases, indicating the presence of oil into the sea, leading to water contamination with oil, and, consequently, financial and environmental damage.
In this study, it was possible to evaluate the transient and isothermal flow behavior of oil in a horizontal pipe, submerged in a sea stream, with and without leakage. From this research, it can be concluded that:
1) The proposed model showed results of pressure field, volumetric fraction and superficial velocity of the phases that represent the physics of the problem;
2) The leakage has a great influence on the total pressure drop;
3) Increasing the oil velocity in the pipeline causes a decrease in pressure drop, where the higher the oil velocity at the inlet of the pipe, the greater the difference in total pressure along the pipe;
4) The oil has a higher velocity in the center of the pipe and lower velocity in the regions close to the walls of the pipe, presenting smaller velocity differences according to the process time;
5) For the position of the leakage orifice, the volume fraction of oil increases of according with time; in other words, at the experimental level, what explains this phenomenon is the pressure difference inside and outside the leakage which is including influenced by a number of forces such as gravity, buoyancy, viscous force and surface tension.
6) The pressure boundary conditions at the pipeline outlet have no significant influence on the behavior of the oil flowing inside the pipe;
7) Have occurred in the regions around to the leakage orifice. The largest variation in the volume fraction of the phases occurs.
The authors would like to thank CAPES, CNPQ and FINEP (Brazilian Research Agencies) for the financial support to this research and the authors cited in the text which help in the improvement.
 Colombaroli, P.L. (2008) Acoustic Leakage Detection in Oil Pipelines. Federal University of Itajubá/MG.
 Biehl, J.B., Quadri, M.B., Santos, M.C. and Uliana, N.R. (2014) Simulation of the Immiscible Flow of an Oil Jet in Underwater to the Interior of a Collection Reservoir. XX Brazilian Chemical Engineering Congress, Florianópolis, 19-22 October 2014, 1-10.
 Yue, P. and Zhou, C. (2006) Phase-Field Simulations of Interfacial Dynamic in Viscoelastic Fluids Using Finite Elements with Adaptive Meshing. Journal of Computational Physics, 219, 47-67.
 Tavares, D.P.S.A. (2016) Fluid Dynamics of Oil Dispersion from Leakage in a Catenary Submarine Riser: Analysis of Multiphase Plumes. Doctoral Thesis Post Graduation in Chemical Engineering, Federal University of Campina Grande, Campina Grande.
 ANP National Agency of Petroleum (2020) Natural Gas and Biofuels.
 Baptista, R.M., Quadri, M.B., Machado, R.A.F., Bolzan, A., Nogueira, A.L. and Mariano, G.C. (2007) Effective Interfacial Tension and Geometrical Parameters Relationship for the Description of Oil Leakages from Submarine Pipelines. 8th International Conference on Chemical & Process Engineering, Vol. 11, 206-401.
 Falcao, M.V.S. (2015) Analysis of Oil Leakage in Waters at Great Depths: A Study with Computational Fluid Dynamics. Course Conclusion Depart, Graduation in Chemical Engineering, Federal University of Rio de Janeiro UFRJ, Rio de Janeiro.
 Zhu, H., You, J. and Zhao, H. (2017) An Experimental Investigation of Underwater Spread of Oil Spill in a Shear Flow. Marine Pollution Bulletin, 116, 156-166.
 Zhu, H., You, J. and Zhao, H. (2017) Underwater Spreading and Surface Drifting of Oil Spilled from Asubmarine Pipeline under the Combined Action of Wave and Current. Applied Ocean Research, 64, 217-235.
 Sun, Y., Cao, X. and Liang, F. (2019) Investigation on Underwater Spreading Characteristics and Migration Law of Oil Leakage from Damaged Submarine Pipelines. Process Safety and Environmental Protection, 127, 329-347.
 Papadimitrakis, J., Pstalki, M., Christolis, M. and Markatos, N.C. (2006) Simulating the Fate of an Oil Spill near Coastal Zones: The Case of a Spill (from a Power Plant) at the Greek Island of Lesvos. Environmental Modelling & Software, 21, 170-177.
 Neofytou, P., Venetsanos, A.G., Vlachogiannis, D., Bartzis, J.G. and Scaperdas, A. (2006) CFD Simulations of the Wind Environment around an Airport Terminal Building. Environmental Modelling & Software, 21, 520-524.
 Wols, B.A., Hofman, J.A.M.H., Uijttewaal, W.S.J., Rietved, L.C. and Van-Dijk, J.C. (2010) Evaluation of Different Disinfection Calculation Methods Using CFD. Environmental Modelling & Software, 25, 573-582.
 Maggiotto, G., Buccolieri, R.R., Antonio-Santo, M., Sandra-Leo, L. and Di-Sabatino, S. (2014) Validation of Temperature-Perturbation and CFD-Based Modelling for the Prediction of the Thermal Urban Environment: The Lecce (IT) Case Study. Environmental Modelling & Software, 60, 69-83.
 Araujo, M.V. (2014) Numerical Study of Biphasic Flow of Water and Oil in T Connections in the Presence of Leakage. Master’s Dissertation, Post-Graduation in Chemical Engineering, Federal University of Campina Grande UFCG, Campina Grande.
 Santos, W.R.G. (2014) Multiphase Flow in the Vertical Leaking Pipeline: Physical/Geometric Evaluation of the Pipeline in the Flow Hydrodynamics. Master’s Dissertation, Post-Graduation in Mechanical Engineering, Federal University of Campina Grande UFCG, Campina Grande.
 Andrade, T.H.F. (2013) Transport of Heavy and Ultravisch Oils via Core-Flow: Geometric and Thermofluidodynamic Aspects. Doctoral Thesis, Post-Graduation in Process Engineering, Federal University of Campina Grande UFCG, Campina Grande.
 Rodriguez, O.M.H. (2002) Interface Shape and Pressure Gradient in the Ascending Vertical Liquid-Liquid Annular Pattern. Doctoral Thesis, Post-Graduation in Mechanical Engineering, Federal University of Campinas, Campinas.
 Acikgoz, M., Franca, F. and Lahey, J. (1992) An Experimental Study of Three Phase Flow Regimes. International Journal of Multiphase Flow, 18, 327-336.
 Li, J., An, W., Gao, H., Zhao, Y. and Sun, Y. (2018) An Experimental Study on Oil Droplet Size Distribution in Subsurface Oil Releases. Acta Oceanologica Sinica, 37, 88-95.